WO2013097413A1

WO2013097413A1 - 一种二倍体单体构建方法和系统

Info

Publication number: WO2013097413A1
Application number: PCT/CN2012/076324
Authority: WO
Inventors: 黄树嘉; 孙鹏; 吴红龙; 汪建; 王俊; 杨焕明
Original assignee: 深圳华大基因科技服务有限公司
Priority date: 2011-12-31
Filing date: 2012-05-31
Publication date: 2013-07-04
Also published as: US20150120256A1

Abstract

本发明公开一种二倍体单体构建方法和系统，涉及生物信息学领域。该方法包括：将根据所有至少包含有一个共同位点的序列片段集合在一起构建由三元字符{A，B，C}组成的m×n的序列片段矩阵M；根据序列片段矩阵M初始化两个片段集合S和T，S∪T=M，并且S∩T=Φ，Φ表示空集；确定目标函数ζ（S，T）=∑∑ε（M，i，j）和初始参考温度T₀，ε（M，i，j）表示片段i和j之间碱基型相同的所有数目和碱基型不同的所有数目之差值；基于目标函数和初始参考温度T₀进行模拟退火过程，达到收敛时输出最终的集合S和T；根据最终的集合S和T通过最小错误纠正模型推断出单体型h。本发明的二倍体单体构建方法和系统，能够获得全局最优解的单体型，准确性高、速度快。

Description

一种二倍体单体构建方法和系统技术领域

本发明涉及生物信息学领域，特别涉及一种二倍体单体构建方法和系统。背景技术

在基因组中不同个体 DNA序列上的同一个位点上单个碱基的差异称为单核苷酸多态（ single nucleotide polymorphism, 简称 SNP ) 。 SNP 是基因组中最常见的遗传变异，据估计人类群体中大约存在一千万左右的 SNP 位点，其中大约 90%是人类群体间共有的。单体型（haplotype )指的是位于一条染色体上或某一区域的一组相关联的 SNP 等位位点。单体型是描述人类基因组遗传差异的一种主要方式，同时也广泛用于基因组关联研究，群体遗传学研究等。

2002 年美国、英国、中国等国家发起了国际单体型图 ( Haplotype Map , HapMap ) 计划，目的是创建一个公共的，基因组范围的常见人类序列变异的数据库，为临床表型的基因研究提供有价值的信息。随着 HapMap I期、 Π期、 m期计划的逐步完成，各国科学家已经对非裔、亚裔、欧裔等十一个不同群体一千多个样本的数百万 SNP 进行了成功分型，构建了不同群体的单体型图谱。除此之外，对人类基因组结构了解的不断深入（包括基因组中的 SNP, 纯合 \杂合缺失，染色体倒位，拷贝数变异等等），使得对单体型的结构和在群体中的分布也有了越来越清晰的认识。同时，随着二代测序技术的曰趋成熟和广泛应用，获得了大量人类基因组中染色体的片段信息。这使可以不用再完全依赖统计学方法基于群体的 SNP 基因型信息重构单个个体的单体型，而可以通过拼接单个个体含有杂合 SNP 位点的染色体片段来直接完成重构。有多篇文献已经证明，基于片段拼接方式构建的单体型具有更高的准确性和完整性。

为了基于染色体片段信息构建单体型，很多研究提出了不同的目标函数用于获得最佳的单体型重构结果，包括最少片段移除 ( minimum fragment removal, MFR ) ，最少错误校正 ( minimum error correction， MEC ) ，最少 SNP 移除 ( minimum SNP removal, MSR )等等。目前，有许多方法以 MEC 为目标函数对单体型进行重构（即使得重构的单体型与已知的染色体片段间差异最小），这些方法有：

1. Levy etal.等提出的贪婪启发式算法：其核心思想是基于贪婪启发式算法使已知的染色体片段与重构的单体型间有最少的差异。当片段中杂合 SNP位点不存在测序错误时，该方法可以很快获得最优的单体型。当杂合 SNP 位点存在测序错误时，该方法耗时较长且结果的准确性较低。

2. HapCUT: 其核心思想是通过初始化单体型和建立染色体片段矩阵计算 SNP位点间的权值（基于 MEC ) ，根据权值大小构建二分图将 SNP 分为两类，多次迭代后取最优，并按照最优的 SNP 分类结果对单体型进行重构。当数据包含较多的染色体片段和杂合 SNP 位点时，该方法运行时间较长且获得的结果通常只是局部最优解而不是全局最优解。

3. ReFHap: 该方法与 HapCUT 类似都是构建二分图，但不是对其中的 SNP 进行分类，而是将所有的染色体片段按照彼此相似程度的高低分为两类，多次迭代取其中差异最高的两类染色体片段集合作为最终结果并根据其对单体型进行重构。该方法虽然耗时较短且具有较高的准确性，但是仍无法摆脱结果容易陷入局部最优解的缺点。

综上所述，开发一种新的准确性高、耗时短、可以得到全局最优解的单体型重构方法是待解决的问题发明内容

^开的一个方面要解决的一个技术问题是提供一种准确性高、速度快的二倍体单体构建方法和系统。

根据本发明的一个方面，提供一种二倍体单体构建方法，包括：根据所有至少包含有一个共同位点的序列片段构建由三元字符 {A，B，C}组成的 m x n 的序列片段矩阵 M，其中，在所述序列片段矩阵 M中将染色体片段中 SNP位点的两个等位碱基分别用 A和 B标记， m为矩阵的行数，表示染色体片段的数目， n为矩阵的列数，表示杂合 SNP位点的数目；根据所述序列片段矩阵 M初始化两个片段集合 S和 T， S U T=M，并且 S fl T= >， φ表示空集；确定目标函数 (S，T) = ∑∑£(M，i，j) 和初始参考温度 T。，其中 S， j 6 T, £(M，i，j)表示片段 i和 j之间碱基型相同的所有数目和碱基型不同的所有数目之差值;基于所述目标函数和所述初始参考温度 T。进行模拟退火过程，达到收敛时输出最终的集合 S和 T; 根据所述最终的集合 S和 T通过最小错误纠正模型推断出单体型 h。

可选地，初始参考温度 T。为 T。=-|厶 _max|/ln(p_r)，其中， I厶 _max|= _max- ζπήη, p_r为初始接受概率，所述 _max和 _min表示根据矩阵 M随机生成 K组由 S和 T两个片段子集所构成的集合分别计算每一组 S和 T的 ζ值中的最大值和最小值， Κ为大于等于 2的自然数。

可选地，基于所述目标函数和所述初始参考温度 Τ。进行模拟退火过程包括:通过高温退火过程锁定最优解范围，高温退火函数为： T=T_Qexp(-_a(j-l)^1/2); 当所述高温退火过程稳定时，转入低温退火过程，低温退火函数为： T=T。exp(-(x(j-k。/p)^1/2); 其中 T。为初始参考温度， j表示退火次数， k。表示高温状态下退火的次数， =0.98, β=1.2。

可选地，通过 Metropolis抽样稳定准则判断是否停止迭代进入退火。

可选地，模拟退火过程收敛准则为：当连续退火预定数次目标函数的值 ζ都保持不变时则认为整个算法已达到最终收敛目。

可选地，三元字符 {A，B，C}为 {0，1，-}。

可选地，根据所述最终的集合 S和 T通过最小错误纠正模型推断出单体型 h 包括：对于每一列 j， j 6 [1 , n] , mj,。表示该列中 0的数目，表示该列中 1 的数目， hj表示该列被推断出来的碱基型；若！！^,。〉!!^，则 hj=0; 若 my >_mj,。，则 hj=l ; 若非以上两种情形，则 hj= - 。

可选地，该方法还包括：对染色体片段中 SNP 位点进行过滤，去掉纯合及包含两个以上等位碱基的 SNP位点。

根据本发明的另一方面，提供一种二倍体单体构建系统，包括：序列矩阵构建模块，用于根据所有至少包含有一个共同位点的序列片段构建由三元字符 {A，B，C}组成的 m x n的序列片段矩阵 M，其中，在所述序列片段矩阵 M 中将染色体片段中 SNP位点的两个等位碱基分别用 A和 B标记， m为矩阵的行数，表示染色体片段的数目， n为矩阵的列数，表示杂合 SNP 位点的数目；初始条件确定模块，用于根据所述序列片段矩阵

M初始化两个片段集合 s和 τ， S U T=M，并且 s n T= >， Φ 表示空集；确定目标函数 (S，T) = ∑∑£(M，i，j) 和初始参考温度 T₀，其中 S， j 6 T, £(M，i，j)表示片段 i和 j之间型相同的所有数目和碱基型不同的所有数目之差值；模拟退火迭代模块，用于基于所述目标函数和所述初始参考温度 T。进行模拟退火过程，达到收敛时输出最终的集合 S和 T; 单体型确定模块，用于根据所述最终的集合 S和 T通过最小错误纠正模型推断出单体型 h。

可选地，初始参考温度 T。为 T。=-|厶 _max|/ln(p_r)，其中， I厶 _max|= _max- _min， Pr为初始接受概率，所述 _max和 ζπήη分别表示根据矩阵 Μ随机生成 Κ组由 S和 Τ两个片段子集所构成的集合分别计算每一组 S和 Τ的 ζ值中的最大值和最小值， Κ为大于等于 2的自然数。

可选地，模拟退火迭代模块包括：高温退火执行单元，用于通过高温退火过程锁定最优解范围，其中，高温退火函数为： T=T。exp(-a(j-l)^1/2); 退火稳定判断单元，用于判断高温退火过程是否稳定，当所述高温退火过程稳定时，转入低温退火过程；低温退火执行单元，用于执行低温退火过程，其中，低温退火函数为： T=T。exp(-a(j-k。/p)^1/2); 其中 T。为初始参考温度， j 表示退火次数， k。表示高温状态下退火的次数， =0.98, β=1·2。

可选地，模拟退火迭代模块通过 Metropolis抽样稳定准则判断是否停止迭代进入退火。

可选地，模拟退火迭代模块采用的模拟退火过程收敛准则为：当连续退火预定数次目标函数的值 ζ都保持不变时则认为整个算法已达到最终收敛目。

可选地，单体型确定模块对于每一列 j， j [1 , n] , _mj,₀ 表示该列中 0的数目， ιη_μ表示该列中 1 的数目， hj表示该列被推断出来的碱基型；若 m ！！!^ ，则判断 hj=0 ; 若 _{mj l} >m_{j 0} , 则判断 hj=l; 若非以上两种情形，则判断 hj= - 。

可选地，该系统还包括： SNP位点过滤模块，用于对染色体片段中 SNP 位点进行过滤，去掉纯合及包含两个以上等位碱基的 SNP位点。

本发明的二倍体单体构建方法和系统，基于目标函数和初始参考温度进行模拟退火过程，收敛时输出最终的集合 S和 T，通过最小错误纠正模型推断出单体型，从而获得全局最优解的单体型，准确性高、速度快。附图说明

图 1示出本发明的二倍体单体构建方法的一个实施例的流程图；

图 2示出本发明的二倍体单体构建方法的另一个实施例的流程图；

图 3示出本发明的二倍体单体构建方法的又一个实施例的流程图；

图 4 示出高温过程退火函数不同 α取值对降温速度的影响；

图 5 示出低温过程退火函数不同 β 取值对降温速度的影响；

图 6示出染色体片段间重叠程度对单体型重构结果准确性的影响；

图 7示出重构出来的单体型的错误率在不同深度下与 SNP 由测序带来的错误率之间的关系；

图 8 示出本发明的二倍体单体构建系统的一个实施例的结构图；

图 9 示出本发明的二倍体单体构建系统的另一个实施例的结构图；

图 10示出二分图示例。具体实施方式

下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性实施例。

本发明的基本思路在于基于模拟退火算法将通过测序所得的序列片段按照差异度构建出二分图从而实现单体型重构的方法，旨在海量测序数据背景下，快速准确的完成（例如人类）基因组中单体型的重构。

一、模拟退火（ simulate annealing ) 算法的简要说明：模拟退火算法来源于物理系统上固体退火的原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而通过徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据 Metropolis准则，粒子在温度 T时趋于平衡的概率为 e-A E/(kT)，其中 E为温度 T时的内能， Δ Ε为其改变量， k为 Boltzmaim常数（数值为： K=l .3806505 x 10^-23 J/K ) 。该算法是一种用于求解大规模优化问题的随机搜索算法，它以优化问题求解过程与物理系统退火过程之间的相似性为基础；优化的目标函数相当于金属的内能；优化问题的自变量组合状态空间相当于金属的内能状态空间；问题的求解过程就是找一个组合状态，使目标函数值最小（或最大）。利用 Metropolis准则并适当地控制温度的下降过程实现模拟退火，从而达到在多项式时间内求解全局优化问题的目标。

重构单体型的问题，是一个求解组合优化的问题。用固体退火模拟组合优化问题，将内能 E模拟为目标函数值 f，温度 T 演化成控制参数 t，即得到解组合优化问题的模拟退火算法：由初始解 i和控制参数初值 t开始，对当前解重复 "产生新解→计算目标函数差→接受或舍弃" 的迭代，并逐步衰减 t 值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。

二、二分图（ Bipartite graph ) 的基本概念：

二分图又称作二部图，是图论中的一种特殊模型。设 G=(V，E)是一个无向图，如果顶点 V可分割为两个互不相交的子集 (V1，V2)，并且图中的每条边（i， j ) 所关联的两个顶点 i 和 j分别属于这两个不同的顶点集 ( VI； j ^ V2 ), 则称图 G为一个二分图，如图 10所示。

图 1示出本发明的二倍体单体构建方法的一个实施例的流程图。

如图 1 所示，步骤 102，根据所有至少包含有一个共同位点的序列片段构建由三元字符 {A，B，C}组成的 m x n的序列片段矩阵 M，其中，在序列片段矩阵 M中将染色体片段中 SNP位点的两个等位 ^分别用 A和 B标记， m为矩阵的行数，表示染色体片段的数目， n 为矩阵的列数，表示杂合 SNP位点的数目。例如，将染色体片段中 SNP位点的两个等位碱基按照 ASCII码的顺序改用 0和 1标记，即 ASCII码较小的用 0表示，较大的则用 1表示，并将所有至少包含有一个共同位点的片段集合在一起构建成一个 m x n 的二维矩阵，记为 M，其中 m为矩阵的行数，表示染色体片段的数目， n 为矩阵的列数，表示杂合 SNP位点的数目，若某染色体片段不包含某一 SNP位点则在矩阵中该点记为 "-" 。从而构建出一个由 {1，0，-}三元字符组成的 m x n的序列片段矩阵。

步骤 104，根据序列片段矩阵 M初始化两个片段集合 S和 τ， S U T=M，并且 s n T= >， Φ表示空集。片段集合 s和 τ可以随^ 取。

步碌 106，确定目标函数 (S，T) = ∑∑ ε (M,i,j)和初始参考温度 T0，其中 S， j Τ， ε (M，i，j)表示片段 i和 j之间型相同的所有数目和 ^型不同的所有数目之差值。假设片段集 S 中所有的序列片段都来自于同一个单体型，而片段集 T中所有的序列片段则完全来自于另外一个单体型，那么此时这两个集合 S 和 T的差异将达到最大。当 (S，T)的值达到最大时，则可以认为 S和 T中所有的片段皆分别来自于不同的两个单体型。稍后将通过一个例子介绍初始参考温度的选取。

步骤 108，基于目标函数和初始参考温度 TO进行模拟退火过程，收敛时输出最终的集合 S和 T。收敛^ Γ以由本领域的技术人员根据需要进行选调整。

步骤 110，根据最终的集合 S 和 T通过最小错误纠正模型 ( MEC )推断出单体型 h。例如，通过最小^纠正模型推断出单体型 h: 对于每一列 j ( j G [l， n] ), mj,。表示该列中 0 的数目， my表示该列中 1 的数目， hj表示该列被推断出来的>?½ 型。若 m >_{mj l}，则

若 my >m_{j 0}，则 h_rl; 若非以上两种情形则 hj= -表示无法确定。

上述实施例中，通过序列片段构建片段矩阵，基于目标函数和初始参考温度进行模拟退火过程，收敛时输出最终的集合 S和 T，通过最小错误纠正模型推断出单体型，从而获得全局最优解的单体型，准确性高。目标函数同时考虑了片段间碱基型的相同与不同，因此还能够有效避免只考虑碱基型不相同位点的分值这种由于信息量利用不足所带来的缺陷。该缺陷可以通过以下例子进行说明：

片段 1011111011111

片段 f₂: 1011001

片段 f₃: 1011010111111 若只是考虑碱基型不相同的位点将由于 f₂信息的不全而致使认为 f₂和 f₃的差异最小，而在这情形下，事实则是和 f₃的差异才是最小的。

下面介绍一个实施例中初始状态（初始参考温度）的确定：初温函数为 T。叫厶 _max|/ln(p_r)，其中 T。表示初始参考温度。利用矩阵 M随机生成 K (例如， 30~200 )组由 S和 T两个片段子集所构成的集合，然后分别计算每一组 S和 T的 ζ值，这 Κ个 ζ值间最大的差值即为 |A_max |，即 |A_{max max}-_min； p_r为初始接受概率（例如， p_r=0.9 )。

图 2示出本发明的二倍体单体构建方法的另一个实施例的流程图。

如图 2 所示，步骤 202，将根据所有至少包含有一个共同位点的序列片段集合在一起构建由三元字符 {0，1，-}组成的 mxn 的序列片段矩阵 M。

步骤 204，根据序列片段矩阵 M初始化两个片段集合 S和

T。

步骤 206，确定目标函数 ζ ,Τ) = ∑∑ ε (M,i,j)o 例如，确定打分矩阵： s{a_x, a₂) ― )

其中， _ai， a₂分别为矩阵 M中第 i条序列片段和第 j条序列片段在同一位点坐标上的^ ^型。片段 i和片段 j 的差异分值记为 c(M,i,j )=∑ε(Μ [i,k] ,M [j ,k] ) , 其中 k G [1，η】，也即是计算片段 i和 j之间^ ^型相同的所有数目与^ ^型不同的所有数目之差值。该打分矩阵的意义在于通过将片段间所有碱基型相同的数目与碱基型不同的数目做差，从而给出片段差异程度的分值，分值越大表示两个片段间的差异就越大。

步骤 208，根据初温函数 -|Amax|/ln(pr)确定初始参考温度 T₀=-|A_max|/ln(p_r)。

退火过程是算法的一个重要过程，影响着状态转换时的接受概率，为了更准确地获得全局最优解避免在退火的后期陷入局部最优解，将退火过程分为两个：高温退火和低温退火过程。

步骤 210，通过高温退火过程锁定最优解范围。高温退火过程的目的在于快速锁定最优解（也即 ζ最大值）的范围，缩小解区间。高温退火函数例如为： T=T。exp(-a(j-l)^1/2)，其中 T₀ 为初始温度， j 是迭代次数， a=0.98。相应的模型扰动方式为：

μ是 [0，1】内均匀分布的随机数， [A^为波动范围，并且 [Αί,Βί] 例如，假设有 50 个序列片段，那么波动范围就为 [1, 50】。

步骤 212，当高温退火过程稳定时转入低温退火过程，达到收敛时输出最终的集合 S和 Τ。当高温退火过程稳定时（判断稳定与否的条件参见下文实施例中的说明），则转入低温退火过程。低温退火函数例如为： T=T。exp(-(x(j-k。/p)^1/2)，其中 T₀为初始温度， j 表示退火次数， ko表示高温状态下退火的次数， α=0.98， β=1.2。相应的模型扰动方式为： 111 =1^+ (8「 Ai)， μ 是 [0,1】内均匀分布的随机数， [AbBi]为波动范围，并且 mi 6 [Ai,Bi] 从该退火函数可以看出，在刚进入低温阶段时，系统进行了一定程度的回火升温，这一升温有利于跳出在高温退火过程可能进入的局部最优解。

步骤 214，根据最终的集合 S和 T通过最小错误纠正模型推断出单体型。对于每一列 j ( j 6 [1， n] ) ， _mj,。表示该列中 0的数目， ιη_μ表示该列中 1 的数目， hj表示该列被推断出来的碱基型。若 m >_{mj l}，则 hj=0; 若则 hj=l; 若非以上两种情形则 hj= - 表示无法确定。当 h确定之后，利用互补的关系就能准确的获得与其相对应的另一个单体型。

上述实施例中，通过高温和低温两次退火过程，设定双阈值的方法来构建二分图，最终就能重构出二倍体基因组的两条单体型 H= ( hl， h2 ) ，既保证了单体型重构的准确性，又能提高收敛速度，耗时短。

退火过程（包括高温和低温）中分为两个紧密相连的步 ( 1 ) Metropolis抽样稳定准则：在同一退火温度 t下目标函数 ζ值达到稳定，则停止迭代进入退火。在同一温度 t 下，采用 Metropolis抽样准则每次都随机从 S或 T 中抽取一条序列片段 v，若 vGS，则 TUv，若 vGT，则 SUv，并计算经此变换后 ζ的值。若 ζ值变大则接受变换，若 ζ值变小或不变则计算此时接受概率函数 ρ(-Δ ζ/Τ)的值，将其与 0~1 之间的随机数比较，判断是否接受变换。具体的 Metropolis抽样准则如下：假设第 i-1次迭代得到的目标函数的值为 ζ。_ld，第 i次经过变换后目标函数的值为 ζ _newJ 二者的差值为 Δ ζ = ζ _new - ζ old. 其中 new ^ 。_ld。生成位于 [0， 1】间的随机数 ρ，若 p<exp(-A ζ/Τ)则接受变换，并令。_ld= _new; 若 ρ>εχρ(-Δ ζ /Τ)则取消变换，并回滚到上一状态。重复上述抽样步骤（m 为序列片段的数目 ) ，当目标函数值连续 W 次（例如，若 m<10， W= 15 m; ^ m>10, W=150)没有变化，认为达到 Metropolis抽样稳定准则，此时停止 Metropolis抽样，返回当前所得最优解以及相对应的状态，接着利用退火函数进行退火降温，然后再次进入 Metropolis抽样过程进行判断。

( 2 ) 算法收敛准则：当连续退火 N 次，（例如，若 m<10， N=5xm次，若 m>10， N=50次， m为矩阵中序列片段的数目），目标函数的值 ζ都保持不变时则认为整个算法已达到最终收敛，将此时得到的两个集合 S和 Τ作为最终最优解输出，并停止算法。

通常来说，在构建序列片段矩阵前可以对染色体片段中

SNP位点进行过滤，去掉纯合及包含两个以上等位碱基的 SNP 位点。

图 3示出本发明的二倍体单体构建方法的又一个实施例的流程图。

如图 3所示，步骤 300，对 SNP位点进行过滤，去掉纯合及包含两个以上等位碱基的 SNP位点。

步骤 302，构列片段矩阵 M。

步骤 304，根据序列片段矩阵 M随;^始化两个片段集合 S 和 T，确定目标函数和初始参考温度。

步骤 306，判断片段集合 S和 Τ是否满足算法收敛准则？如果是，输出结果（步骤 326 )，否则，继续步骤 308。

步骤 308，判断是否满足抽样稳定准则？如果是，则继续步骤 310，否则，继续步骤 314。

步骤 310，判断当前是否属于高温退火过程？如果是，则继续步骤 312a，采用高温退火函数进行退火；否则，继续步骤 312b, 采用低温退火函数进行退火。

步骤 314，由当前状态产生新状态。

步骤 316，判断接受函数是否成立，如果是，则接受新状态 (步骤 320 )，否则，保持当前状态不变（步骤 318 )。返回步骤 308。

需要指出，在图 3中步骤的序号并不表示该步骤的执行先后顺序。在步骤 308判断满足了 "抽样稳定准则" 之后进入退温过程，而退温的时候需要选择高温退火过程还是低温退火过程，这两个过程用到的是两个不同的退温函数，步骤 310判断这个时候是高温退火过程（步骤 312a )还是低温退火过程（步骤 312b )，如果为高温过程就用高温退火函数，如果是低温过程就用低温退火函数。不过需要特别说明的一点就是：算法在重构单体型的过程中高温和低温过程是不会交叉出现的，一定是严格分开的两个过程，即在经过高温退火过程之后，接下来的就是低温退火过程，不会出现从低温退火过程进入高温退火过程的现象，回温步骤也算低温过程。

图 3 示出了以 ζ为目标函数基于模拟退火的单体型重构方法流程，下面将以实例具体阐述退火过程中相应参数的设定及适合的数据类型。

实施例参数 α和 β的取值：退火函数中参数 α和 β的不同取值将影响降温的速度。申请人模拟了在初温 Τ。=100，降温次数 j 6 [1,100] , 取值分别为 0.98、 0.95、 0.9 时，退火温度 Τ 下降的情况。从图 4 中可以看到当（X值减小时，退火温度下降速度无明显差异。但为了加快高温退火的速度以减少迭代次数，在本发明中选取 α=0.98。申请人模拟了在初温 Τ。=100，降温次数 j [51,150】，高温退火次数 k。=50， β取值分别为 1、 1.2、 1.5时，退火温度 Τ下降的情况。从图 5可以看到当 β值减小时，退火温度下降速度差异明显。回火升温有助于结果跳出局部最优解，为了不使低温退火速度过快在本发明中 β=1·2。

实施例一：结果评估指标 SE : 申请人使用交变错误率 ( switch error, SE , 也可称为重构错误率）来评估基于本发明的单体型重构结果的准确性。计算 SE的公式为：

SE⁼min{d(h_rec0n_Struct， h_rean), d(h_reconstruct， h_real2)}/n，其中 hreconstruct表示重构结果中的一条单体型， d(h reconstruct? h_reail)和 d(h_reconstruct, h_real2)表示一条重构单体型与模拟产生的两条标准单体型间 SNP 不匹配的数目， n表示为重构结果中 SNP的数目，那么 SE就表示为重构结果与模拟的真实结果间 SNP位点不一致最小数目的百分比。 SE 值越小则表明根据模拟数据重构的单体型与真实结果越相似，准确性就越高。

先评价 SE 与序列片段重叠程度的关系：染色体片段重叠程度 ( overlap level )定义为：

Poverlap_level=∑ N。_V x )，

其中， m表示序列片段的数目（即上述序列矩阵 M 的行数）， n表示 SNP位点的数目（即上述序列矩阵的列数）， ∑ N。_verlap表示所有存在于两条以及两条以上序列片段的 SNP位点深度之总和。如下表 1中的例子:

表 1

∑ N_overia_P = 4+4+2 = 10

m=3

n=13

P。_verlap__level反映了整个序列矩阵 M 中那些存在重叠关系的序列片段以及 SNP 位点所占到的比重及可利用信息的充裕程度，这个比重既能说明每个 SNP 位点平均被重用到的次数还能说明序列片段间的紧凑程度，片段间相互重叠的 SNP 位点越多（即 SNP 位点平均被重用到的次数越大）或序列片段间约紧凑，那么 P。_verlap__level就越大，可利用的信息也越充裕。

申请人随机生成了 39 套模拟数据，用于评估染色体片段间重叠程度与 SE的关系。每套数据由随机产生的 50对标准单体型及根据单体型生成的染色体片段组成。每套数据间，标准单体型数据包含的 SNP数目相同为 200，相应生成的染色体片段数相同为 20，但染色体片段包含的 SNP数目则是由 10开始每套增加 5。为了使模拟数据更接近真实数据，还考虑了染色体片段中 SNP 的缺失及颠换。在生成模拟数据过程中染色体片段每个杂合 SNP位点缺失的概率为 0.9，这个概率是远高于实际情况的，在这个条件下假如同样能得到好的结果那么就说明本发明的方法是很有实际应用意义的，同时，设置发生颠换的概率为 0.05。最后，通过统计生成的 39 套模拟数据 Poverlapjevel (即重叠程度）的范围为 0·18~0·90，相应的 SE (即重构错误率）范围为 0~0.03。见图 6，横坐标为重叠程度（下称 Ρ overlap level ) ，纵坐标为重构错误率（下称 SE ) ，可以看出，随着 P overlap— level的升高 SE 迅速下降，要注意的是，在这种数据之下重构错误率最高也不到 0.03，相当于准确性至少是 97% , 这个准确性已经相当高了，而当 P。_verlap__level达到 20% (也就是图中的 0.2 ) 以上时， SE已经在 0.01以下。这表明本发明中实施例的方法在进行单体型重构时可以获得极高的准确度。

实施例二： SE与 SNP颠换率的关系：为了评估 SNP颠换率与 SE间的关系，申请人生成了 SNP在不同覆盖深度下颠换率由低到高的模拟数据。

SNP覆盖深度包括 10X, 20X, 30X, 40X, 50X, 每个深度又包含 SNP颠换率为 0.01， 0.05, 0.1 , 0.2, 0.3, 0.4, 0.5 的 7套模拟数据，每套数据由随机产生的 50对标准单体型及根据单体型生成的染色体片段组成。 SNP覆盖深度的计算公式为：

C=m L (L-d)/N,

其中 m表示染色体片段数， L 表示染色体片段的平均长度， d表示 SNP 缺失率， N表示标准单体型含有的 SNP 数目。对于不同深度的模拟数据 N=200， 1=20， m分别为 110， 220， 330, 440， 550。通过计算不同深度下 SE和 SNP错误率的值，如图 7 所示，从图中可以很明显得看出随着 SNP 由于测序导致的错误率的增加单体型重构的错误率也在不断增加，特别是当 SNP错误率达到 0.3以上之后重构的错误率就直接增大到 0.25以上，这是不可接受的，但是一般情形下，测序并不会带来如此大的错误率，一般都在 0.001 以下，特殊请况下个别最大到达 0.05，因此本发明实施例的方法同样具备实际意义。

图 8 示出本发明的二倍体单体构建系统的一个实施例的结构图。如图 8所示，该实施例中二倍体单体构建系统包括：序列矩阵构建模块 81，将根据所有至少包含有一个共同位点的序列片段集合在一起构建由三元字符 {A，B，C}组成的 m X n的序列片段矩阵 M，其中，在序列片段矩阵 M 中将染色体片段中 SNP位点的两个等位碱基分别用 A和 B标记， m为矩阵的行数，表示染色体片段的数目， n为矩阵的列数，表示杂合 SNP 位点的数目；初始条件确定模块 82，根据序列片段矩阵 M初始化两个片段集合 s和 τ， S U T=M，并且 s fi T= >， Φ表示空集；确定目标函数 (S，T) =∑∑£(M，i，j) 和初始参考温度 T₀，其中 S， j 6 T, £(M，i，j)表示片段 i和 j之间型相同的所有数目和碱基型不同的所有数目之差值；模拟退火迭代模块 83, 基于目标函数和初始参考温度 T。进行模拟退火过程，达到收敛时输出最终的集合 S和 T; 单体型确定 84，根据最终的集合 S和 T通过最小错误纠正模型推断出单体型 h。在一个实施例中，单体型确定模块 84对于每一列 j， j [1， n] , _mj,₀表示该列中 0的数目， my表示该列中 1 的数目， hj表示该列被推断出来的碱基型；若 m >_{mj l} ，则判断 hj=0；若 _{mj l} >m_{j 0} , 则判断 hj=l; 若非以上两种情形，则判断 hj= - 。

在一个实施例中，初始参考温度 T。通过如下方式确定： T₀=-| A _max|/ln(p_r)，其中， | A _max|= _max- _min， p_r为初始接受概率， max和 _min分别表示根据矩阵 Μ随机生成 Κ组由 S和 Τ 两个片段子集所构成的集合分别计算每一组 S和 Τ的 ζ值中的最大值和最小值， Κ为大于等于 2的自然数。

在一个实施例中，模拟退火迭代模块通过 Metropolis抽样稳定准则判断是否停止迭代进入退火；模拟退火迭代模块采用的模拟退火过程收敛准则为：当连续退火预定数次目标函数的值 ζ都保持不变时则认为整个算法已达到最终收敛目。

图 9 示出本发明的二倍体单体构建系统的另一个实施例的结构图。如图 9 所示，该实施例中，还可以包括 SNP位点过滤模块 90，用于对染色体片段中 SNP位点进行过滤，去掉纯合及包含两个以上等位碱基的 SNP位点。在一个实施例中，模拟退火迭代模块 93 包括:高温退火执行单元 931，通过高温退火过程锁定最优解范围，其中，高温退火函数为： T=T。exp (- a(j-l)^1/2); 退火稳定判断单元 932，判断高温退火过程是否稳定，当所述高温退火过程稳定时，转入低温退火过程；低温退火执行单元 933，执行低温退火过程，其中，低温退火函数为： T=T₀exp(-a(j-k₀/p)^1/2); 其中 T。为初始参考温度， j表示退火次数， k。表示高温状态下退火的次数， α=0.98, β=1.2。

对于图 8和图 9中名斗装置或单元的功能，可以参考上文中关于本发明方法的实施例中对应部分的说明，为简洁起见，在此不再伴述。

本领域的技术人员应当理解，对于图 8 至图 9 中的各个装置，可以通过单独的计算处理设备实现，或者将其集成为一个独立的设备实现。在图 8至图 9中用框示出以说明它们的功能。这些功能块可以用硬件、软件、固件、中间件、微代码、硬件描述语音或者它们的任意组合来实现。举例来说，一个或者两个功能块都可以利用运行在微处理器、数字信号处理器（DSP )或任何其他适当计算设备上的代码实现。代码可以表示过程、功能、子程序、程序、例行程序、子例行程序、模块或者指令、数据结构或程序语句的任意组合。代码可以位于计算机可读介质中。计算机可读介质可以包括一个或者多个存储设备，例如，包括 RAM 存储器、闪存存储器、 ROM 存储器、 EPROM 存储器、 EEPROM存储器、寄存器、 ^ i、移动硬盘、 CD-ROM或本领域公知的其他任何形式的存储介质。计算机可读介质还可以包括编码数据信号的载波。

本领域技术人员将意识到硬件、固件和软件配置在这些情况下的可替换性，以及如何最好地实现每个特定应用地该功能。

本专利在已有的成熟测序技术及单体型重构方法的基础上提出一种新的单体型重构方法和系统，可以在较短时间内确定目标函数的全局最优解，进而完成单体性的重构。

本发明的描述是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显然的。选# ^描述实施例是为了更好说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。

Claims

权利要求

1. 一种二倍体单体构建方法，其特征在于，包括：根据所有至少包含有一个共同位点的序列片段构建由三元字符 {A，B，C}组成的 m x n 的序列片段矩阵 M，其中，在所述序列片段矩阵 M中将染色体片段中 SNP位点的两个等位分别用 A和 B标记， m为矩阵的行数，表示染色体片段的数目， n为矩阵的列数，表示杂合 SNP位点的数目；

根据所述序列片段矩阵 M初始化两个片段集合 S和 T， S U T=M，并且 S fl T= >， φ表示空集；

确定目标函数 (S，T) =∑∑£(M，i，j) 和初始参考温度 To, 其中 S， j 6 T, £(M，i，j)表示片段 i和 j之间型相同的所有数目和^ ^型不同的所有数目之差值；

基于所述目标函数和所述初始参考温度 T。进行模拟退火过程，达到收敛条件时输出最终的集合 S和 T;

根据所述最终的集合 S和 T通过最小错误纠正模型推断出单体型 h。

2. 根据权利要求 1 所述的方法，其特征在于，所述初始参考温度 Τ。为 T。=-|厶 _max|/ln(p_r)，其中， | A _max|= _max— ζπΰη ' p_r 为初始接受概率，所述 _max和 _min表示根据矩阵 M随机生成 K 组由 S和 T两个片段子集所构成的集合分别计算每一组 S和 T 的 ζ值中的最大值和最小值， Κ为大于等于 2的自然数。

3. 根据权利要求 1 所述的方法，基于所述目标函数和所述初始参考温度 Τ。进行模拟退火过程包括：

通过高温退火过程锁定最优解范围，高温退火函数为： T=T₀exp(-a(j-l)^{1 2});

当所述高温退火过程稳定时，转入低温退火过程，低温退火函数为： Τ=Τ χρ(-αα- β)¹²); 其中 T。为初始参考温度， j表示退火次数， k。表示高温状态下退火的次数， =0.98, β=1.2。

4. 根据权利要求 3 所述的方法，其特征在于，通过 Metropolis抽样稳定准则判断是否停止迭代进入退火。

5. 根据权利要求 1 所述的方法，其特征在于，所述模拟退火过程收敛准则为：

当连续退火预定数次目标函数的值 ζ都保持不变时则认为整个算法已达到最终收敛目。

6. 根据权利要求 1 所述的方法，其特征在于，所述三元字符 {A，B，C}为 {0，1，-}。

7. 根据权利要求 1 所述的方法，其特征在于，所述根据所述最终的集合 S和 T通过最小错误纠正模型推断出单体型 h 包括：

对于每一列 j， j ^ [1 , n] , mj,。表示该列中 0的数目， my 表示该列中 1的数目， hj表示该列被推断出来的型；

若 m_{j 0} > _hl，则 hj=0;

若 mj ! >m_{j 0}，则 hj=l;

若非以上两种情形，则 hj= - 。

8. 根据权利要求 1所述的方法，其特征在于，还包括：对染色体片段中 SNP位点进行过滤，去掉纯合及包含两个以上等位威基的 SNP位点。

9. 一种二倍体单体构建系统，其特征在于，包括：序列矩阵构建模块，用于根据所有至少包含有一个共同位点的序列片段构建由三元字符 {A，B，C}组成的 m x n的序列片段矩阵 M，其中，在所述序列片段矩阵 M 中将染色体片段中 SNP位点的两个等位碱基分别用 A和 B标记， m为矩阵的行数，表示染色体片段的数目， n为矩阵的列数，表示杂合 SNP 位点的数目；初始条件确定模块，用于根据所述序列片段矩阵 M 初始化两个片段集合 s和 τ， S U T=M，并且 s fi T= >， Φ表示空集；确定目标函数 (S，T) =∑∑£(M，i，j) 和初始参考温度 T。，其中 S， j 6 T, £(M，i，j)表示片段 i和 j之间型相同的所有数目和^ ^型不同的所有数目之差值；

模拟退火迭代模块，用于基于所述目标函数和所述初始参考温度 T。进行模拟退火过程，达到收敛时输出最终的集合 S 和 τ;

单体型确定模块，用于根据所述最终的集合 S和 T通过最小错误纠正模型推断出单体型 h。

10. 根据权利要求 9所述的系统，其特征在于，所述初始参考温度 Τ。为 T。=-|厶 _max|/ln(p_r)，其中， | A _max|= _max— ζπΰη ' p_r 为初始接受概率，所述 _max和 _min分别表示根据矩阵 M随机生成 K组由 S和 T两个片段子集所构成的集合分别计算每一组 S 和 T的 ζ值中的最大值和最小值， Κ为大于等于 2的自然数。

11. 根据权利要求 9所述的系统，所述模拟退火迭代模块包括：

高温退火执行单元，用于通过高温退火过程锁定最优解范围，其中，高温退火函数为： T=T。exp(-a(j-l)^1/2);

退火稳定判断单元，用于判断高温退火过程是否稳定，当所述高温退火过程稳定时，转入低温退火过程；

低温退火执行单元，用于执行低温退火过程，其中，低温退火函数为： T=T。exp(-a(j-k。/p)^1/2);

其中 T。为初始参考温度， j表示退火次数， k。表示高温状态下退火的次数， =0.98, β=1.2。

12. 根据权利要求 9所述的系统，其特征在于，所述模拟退火迭代模块通过 Metropolis抽样稳定准则判断是否停止迭代进入退火。

13. 根据权利要求 9所述的系统，其特征在于，所述模拟退火迭代模块采用的模拟退火过程收敛准则为：

14. 根据权利要求 9所述的系统，其特征在于，所述单体型确定模块对于每一列 j， j ^ [1, n] , _mj,。表示该列中 0 的数目，表示该列中 1 的数目， hj表示该列被推断出来的>?½ 型；

若！^,。〉!!!^，则判断 hj=0;

^ m _i >_mj,₀，则判断 hj=l ;

若非以上两种情形，则判断 hj= - 。

15. 根据权利要求 9所述的系统，其特征在于，还包括： SNP位点过滤模块，用于对染色体片段中 SNP位点进行过滤，去掉纯合及包含两个以上等位碱基的 SNP位点。