2. 广西师范大学广西多源信息挖掘与安全重点实验室, 广西 桂林 541004;
3. 广西区域多源信息集成与智能处理协同创新中心, 广西 桂林 541004
2. Guangxi Key Lab of Multi-source Information Mining & Security, Guangxi Normal University, Guilin 541004, Guangxi, China ;
3. Guangxi Collaborative Innovation Center of Multi-source Information Integration and Intelligent Processing, Guilin 541004, Guangxi, China
单核苷酸多态性(single nucleotide polymorphism, SNP)作为基因组上单个核苷酸碱基变异频率大于1％的DNA序列多态性, 广泛存在于人类基因组中, 其对解释个体表型差异和疾病易感性等方面起重要作用。单体型由染色体上某区域一组连锁遗传的SNP位点组成, 研究发现其比单个SNP位点携带更多的遗传信息, 为个性化诊断及药物设计等的研究提供重要信息。然而, 采用生物实验方法直接获取单体型成本很高, 通常采用计算方法将测序片段组装成单体型, 单体型组装问题应运而生, 也称为单体型重建问题。
针对二倍体单体型的重建问题, 已提出一些理论和算法[1-6]。LANCIA G等人2001年首次提出个体单体型重建问题, 给出求解该问题的最少片段删除模型(the minimum fragment removal, MFR)、最少SNP删除模型(the minimum SNP removal, MSR)和最长单体型重建问题模型(the longest haplotype reconstruction, LHR)[1]。2002年, LIPPERT R等提出求解该问题的最少错误更正模型(the minimum error correction, MEC)[2]。针对这些模型, 研究者们提出一系列精确算法和启发式求解方法, 取得较好的重建效果[7-16], 其中包括著名的迭代算法DGS[4]及启发式算法Fast Hare[12]。2012年, AGUIAR D等利用SNP测序片段构造带权无向图Compass, 提出基于Compass图重建单体型的MWER模型, 并给出求解模型的启发式算法HapCompass[5]。该算法在Compass图中构造最大生成树, 通过迭代删除权值绝对值最小的边以消除图中的冲突环基, 以得到唯一的单体型分型方式。然而, 当1个冲突环基中同时存在多条权值绝对值最小的边时, 随机选择会保留测序错误, 从而降低重建效果。本研究针对这种情况, 提出新的去边规则, 进一步通过各分型的片段支持差异数与总片段数之间的比值来确定删除边。此外, 针对HapCompass算法未对Compass图中孤立SNP位点赋值的问题, 本研究给出有效的解决方案。试验结果显示, IHapCompass算法能有效求解MWER模型, 获得较算法HapCompass[5]、DGS[4]和Fast Hare[12]更高的重建精度。
1 问题模型及符号定义给定1组SNP位点, 可以确定来自1对同源染色体上的单体型。由于SNP有二态性, 常采用0和1编码其野生型和变异型取值, 则单体型可用定义在字符集｛0, 1｝上的二进制序列表示, 而无需用｛A, T, C, G｝ 4个碱基符。若1对单体型在某SNP位点上的碱基取值相同, 则该位点称为纯合位点, 否则称为杂合位点[17]。
单体型重建问题的输入为1组来自1对同源染色体、联配好的SNP片段所构成的二维矩阵Fm×n, 其每行表示1条SNP片段fi, 每列对应1个SNP位点sj, 每个元素fij∈{0, 1, -}(i=1, 2, …, m, j=1, 2, …, n) (-表示片段fi未覆盖SNP位点sj, 或者片段fi在位点sj的取值未知)。令Nxj为第j (j=1, …, n)列中取值为x(x∈{0, 1})的元素数量, Pxj为第j列中取值为x的元素在该列所有非空元素中所占比例[18]:
| $ P_x^j = \frac{{N_x^j}}{{N_x^j + N_{1 - x}^j}},\left( {x \in \left\{ {0,1} \right\},j = 1,\cdots ,n} \right). $ | (1) |
给定两条字符串X=〈x1, …, xn〉和Y=〈y1, …, yn〉, 其中xj, yj∈｛0, 1, -｝ (j=1, …, n), Def(X, Y)定义为X和Y在对应位点取值相异的位数之和:
| $ {\rm{Def}}\left( {X,Y} \right) = \sum\limits_{j = 1}^n {d\left( {{x_j},{y_j}} \right)} , $ | (2) |
式中,
| $ d\left( {{x_j},{y_j}} \right) = \left\{ \begin{array}{l} 1\quad 若{x_j} \ne - ,{y_j} \ne - ,且{x_j} \ne {y_j},\\ 0\quad 若则. \end{array} \right. $ | (3) |
若将X和Y看成两条SNP片段, Def(X, Y)等于0表示片段X和Y是相容的, Def(X, Y)大于0表示片段X和Y是冲突的, 即它们不是源于同1条染色单体或片段数据中存在测序错误。
将矩阵Fm×n中所有杂合位点及覆盖该位点的片段抽取出来组成新矩阵, 仍记为Fm×n, 这里假设每条片段至少覆盖两个杂合SNP位点。由新矩阵F生成带权无向图GC(F)=(VC, EC, w), 其中顶点表示SN位点, 即sj∈VC, 边(si, sj)∈EC表示至少有1条片段同时覆盖SNP位点si和sj(1≤i≤n, 1≤j≤n), 每条边(si, sj)赋权值w(si, sj)。该图称为Compass图[5]。对任意给定的杂合SNP位点si和sj, 二倍体生物单体型存在两种确定的分型, 即
| $ \begin{array}{c} w\left( {{s_i},{s_j}} \right) = \left( {{n_e}\left( {00} \right) + {n_e}\left( {11} \right)} \right) - \\ \left( {{n_e}\left( {01} \right) + {n_e}\left( {10} \right)} \right), \end{array} $ | (4) |
式中ne(x, y)表示矩阵F中, 在位点si和sj上取值分别为x和y的片段数目, 见图 1。
|
图 1 SNP矩阵及Compass图 Figure 1 SNP Matrix and Compass graph |
图 1(a)表示SNP矩阵F, 其中包含5条片段和6个SNPs。图 1 (b)表示GC(F),
给定GC(F)中的1个简单环, 若其中包含奇数个负权值边或至少包含1条0权值边, 则称该简单环为冲突环(conflicting cycle)[5]。当图GC(F)中不存在冲突环时, 可为图GC(F)中的每个顶点唯一地确定其取值, 此时称图GC(F)可唯一分型(unique phasing)[5]。
基于上述定义, 2012年, AGUIAR D等人提出求解重建二倍体个体单体型的MWER模型[5]:给定Compass图GC(F), 试图删除边集L⊂EC, 以满足以下3个条件:
(1)被删除边的权值绝对值之和最小, 即
(2)新图GC′(F)=(VC, EC′, w)(EC′=EC-L)中的每条边都是决定性的;
(3)按多数决定原则确定EC′中每条边的分型可使图GC′(F)唯一分型。
2 IHapCompass算法针对MWER模型提出求解二倍体个体单体型重建问题的改进环基算法IHapCompass (Improved HapCompass)。算法输入为1个m×n的SNP片段矩阵F, 输出为1对长度为n的单体型h=(h1, h2)。IHapCompass首先预处理SNP矩阵F, 删除矩阵中的纯合SNP位点, 并建立Compass图; 然后删除权值绝对值之和最小的边集以使得图中每条边都是决定性的, 且图是可唯一分型的, 以此得到1对仅含杂合位点的单体型h′=(h1′, h2′); 最后将单体型h′扩展成为h=(h1, h2)。为方便描述, 下面给出生成树环基的概念, 然后介绍算法IHapCompass。
给定图GC(F)的最大生成树T, 由1条非树边与T组成的简单环称为图GC(F)的1个生成树环基, 则图GC(F)有|EC|-(|VC|-1)个生成树环基[5]。
2.1 预处理利用文献[12]的列删除法处理矩阵F, 删除所有满足条件Pxj≤e(x∈ ｛0, 1｝, j=1, …, n)的列, 并称该列为(1-x)列, e设为0.2。列删除后, 某些行会变为空行, 也将其删除。预处理后保留的SNP位点均看作杂合位点。统计被删除的1列和0列的数目, 分别记为C1和C0。新矩阵仍用Fm×n表示, 根据新矩阵F构建Compass图GC(F)。
2.2 求解MWER模型给定Compass图GC(F), 首先, 从中删除所有权值为0的边, 该操作不会影响目标值
(1)随机选择1个冲突的生成树环基, 从中删除1条权值绝对值最小的边e。权值绝对值最小表示支持
| $ {\rm{pw}}\left( e \right) = \frac{{\left| {w\left( e \right)} \right|}}{{\sum\limits_{x,y \in \left\{ {0,1} \right\}} {{n_e}\left( {xy} \right)} }}, $ | (5) |
pw(·)体现在某对SNP位点上支持
(2)若删除边e为非生成树边, 则生成树仍有效; 否则将1条非生成树边ent加入其中, 形成新的生成树T-e∪ent。计算图GC(F)-e中相对于树T-e∪ent的所有生成树环基, 标记其是否为冲突环。令GC(F)=GC(F)-e, T=T-e∪ent, 若仍存在冲突的生成树环基, 则返回(1)。
(3)根据图GC(F)边上权值的正负情况确定非孤立点的SNP取值, 得到1条仅含杂合位点的单体型, 并根据杂合位点的取值互补性导出1对仅含杂合位点的单体型。对图中所有孤立点, 根据C0和C1的取值情况为其赋值, 即若C0>C1, 孤立点赋值为0, 否则赋值为1。由此得到单体型h′=(h1′, h2′)。
2.3 扩展结果将预处理时删除的纯合位点重新加回[12]。对于单体型h′=(h1′, h2′), 若某个已删除位点记为1列(或0列), 则将1(或0)插回单体型h1′和h2′的相应位点处, 由此扩展成最终的单体h=(h1, h2)。根据上述算法设计, 算法IHapCompass的详细描述如下:
算法 IHapCompass
输入 SNP矩阵Fm×n,
输出 单体型h=(h1, h2)。
算法步骤:
(1)预处理矩阵F, 得到新矩阵Fm×n, 构造图GC(F);
(2)删除GC(F)中所有权值为0的边;
(3)构建图GC(F)的最大生成树T, 确定图GC(F)中相对T的所有冲突生成树环基;
(4)随机选择1个冲突生成树环基, 从中删除1条|w(·)|值最小的边e。若同时存在多条权值绝对值最小的边, 则删除pw(·)值最小的边e, GC(F)=GC(F)-e;
(5)若被删除的边e∈ T, 则T=T-e∪en(en为非生成树边);
(6)若图GC(F)中仍存在相对T的冲突的生成树环基, 返回步骤4;
(7)确定图中各SNP位点的取值, 得到单体型h′=(h1′, h2′);
(8)扩展单体型h′=(h1′, h2′), 输出h=(h1, h2)
2.4 算法复杂性算法IHapCompass在预处理的时间复杂度为O(mn), 构造图GC(F)的时间复杂度为O(mn2); 删除权值为0的边用时O(n2); 构造生成树在O(n2lbn)时间内完成; 确定生成树环基的时间复杂度为O(n4); 处理生成树的所有简单环花费时间O(mn4); 第6步时间的复杂度为O(n4); 第7步图中的每条边花费至多O(n2)时间, 处理孤立点花费至多O(n)时间; 扩展单体型的时间复杂度为O(n)。算法总时间复杂度为O(mn4)。
3 试验结果试验在联想工作站(英特尔酷睿i5 3.10 GHz, 内存为4.0 GB)上进行, 操作系统为Windows 7, 程序编译器为Microsoft Visual C＃ 2013。
3.1 试验数据原始单体型来自于国际人类基因组单体型图计划[19] (The International HapMap Project)发布的genotypes_chr1_CEU_r22_nr.b36_fwd_phase.gz数据文件, 其包含CEPH样本中60个个体在1号染色体上的单体型。随机选取某个体长度为n的1对单体型作为试验单体型数据。
利用CELSIM[20]和MetaSim[21]两种模拟器生成测序片段数据, 分别称为CELSIM实例和MetaSim实例。生成CELSIM实例的方法如下:对每组测试数据, 生成m1条单片段和m2条mate-pair片段, 单片段的长度范围为[fmin, fmax], mate-pair片段长度为n／10。两种片段的覆盖率均为c／2, 总片段覆盖率为c。根据错误概率Ps植入测序错误。fmin和fmax分别为3和7, c取值范围为5～10, Ps取值范围为2％～5％[12]。
生成MetaSim实例的方法如下:对每组测试数据, 生成m1=(1-Pm)×m条单片段和m2=Pm×m条mate-pair片段, m为片段总数, Pm为mate-pair片段占总片段数目的比例, 试验中取值0.25。单片段的期望长度为l, mate-pair片段的期望长度为3l。由于1条mate-pair片段由同条单体型的1对单片段构成, 则片段覆盖率c=[(m1+2m2)×l]／2n。
3.2 评价指标采用平均重建率和运行时间评价算法性能。重建率(reconstruction rate, RR)用于衡量重建的准确度, 给定真实单体型h=(h1, h2)和重建单体型
| $ {\rm{RR}}\left( {h,\hat h} \right) = 1 - \frac{{\min \left\{ {{r_{11}} + {r_{22}},{r_{12}} + {r_{21}}} \right\}}}{{2n}}, $ | (6) |
式中
测序过程中, 染色单体被随机打断, 某些片段仅重叠于纯合位点。由于1对单体型之间大约仅有20％的杂合位点, 且测序片段受到长度限制, 片段重叠于纯合位点的几率很大[12]。这使得无法根据片段的重叠情况构造全基因组完整的单体型, 取而代之的是构造单体型块(haplotype block)。基于上述原因, 试验中首先根据片段数据对位点的覆盖情况将SNP位点分成nb块, 并分别重建各块单体型, 其重建率为RRi(i=1, …, nb), 则平均重建率
| $ \overline {{\rm{RR}}} = \frac{{\sum\limits_{i = 1}^{{n_{\rm{b}}}} {{\rm{R}}{{\rm{R}}_i}} }}{{{n_{\rm{b}}}}}. $ | (7) |
在不同的测序参数下, 对算法IHapCompass、HapCompass[5]、DGS[4]和Fast Hare[12]的性能进行测试分析。每组测试参数生成100个数据集, 试验结果取100次重复测试的平均值。表 1～4对CELSIM实例进行测试。表中IIHapCompass和HapCompass简略为IHap和Hap。表 1针对覆盖率c设置9组参数, 其中n=100, fmin=3, fmax=7, Ps=0.05。随着c的逐渐增大, 4个算法的重建率均有所提高。这是因为覆盖率较大时, 重建所能运用的原始信息量较多, 从而得到较好的重建效果。同时, 在各c值下, 算法IHapCompass均获得较另3个算法更高的重建率。如表 1所示, 当c从2增加到10时, IHapCompass的重建率从0.87提升到0.99, 而其他3个算法提升至0.96。4种算法均具有较高的执行效率, 且覆盖率对算法运行时间影响不大。
|
表 1 CELSIM实例在不同覆盖率下的比较 Table 1 Comparison with different coverages for CELSIM instance |
表 2针对错误率Ps设置10组参数, 其中fmin=3, fmax=7, c=10, n=100。表 2显示, 算法IHapCompass能获得较其他3个算法更高的重建率, 且错误率对FastHare算法的影响更为显著。当错误率由0.02提高到0.2时, FastHare算法的重建率下降12.5％, 而IHapCompass算法的重建率仅下降6.06％。此外, 在各参数设置下, 4种算法均具有较高的运行效率。
|
|
表 2 CELSIM实例在不同错误率下的比较 Table 2 Comparison with different error rates for CELSIM instance |
表 3针对单体型长度n设置6组测试实例, 其中c=10, fmin=3, fmax=7, Ps=0.05。随着n的增大, 4种算法的重建率都有所下降。由于单体型重建长度越长, 由随机取值带来的误差会越大, 算法IHapCompass的去边规则减少了随机取值的几率, 因此能保证更高的重建率。
|
|
表 3 CELSIM实例在不同单体型长度下的比较 Table 3 Comparison with different haplotype lengths for CELSIM instance |
表 4针对单片段长度范围[fmin, fmax]设置3组参数, 其中c=10, Ps=0.05, n=100。随着单片段长度的减小, 片段重叠概率下降, 使得算法的重建率也呈现下降趋势。然而, 当[fmin, fmax]取值为[1, 2]时, 算法IHapCompass仍获得0.96的重建率。此外, 单片段长度对DGS算法的运行时间有影响外, 对其他3个算法影响不大。
|
|
表 4 CELSIM实例在不同单片段长度范围下的比较 Table 4 Comparison with different single fragment length ranges for CELSIM instance |
表 5～7对MetaSim实例进行测试, 在不同覆盖率、不同单体型长度和不同单片段长度设置下对算法性能进行比较。表 5针对覆盖率设置9组测试实例, 其中n=100, Pm=0.25, l=5。表 6针对不同的单体型长度n设置6组测试实例, 其中c=20, Pm=0.25, l=5。表 7针对于不同单片段设置3组测试实例, 其中n=100, c=20, Pm=0.25。
|
|
表 5 MetaSim实例在不同覆盖率下的比较 Table 5 Comparison with different coverages for MetaSim instance |
|
|
表 6 MetaSim实例在不同单体型长度下的比较 Table 6 Comparison with different haplotype lengths for MetaSim instance |
|
|
表 7 MetaSim实例在不同单片段长度下的比较 Table 7 Comparison with different single fragment lengths for MetaSim instance |
基于最少带权边删除模型, 提出一种重建二倍体个体单体型的启发式算法IHapCompass。该算法通过消除图中的冲突环基以确定单体型分型方式, 其改进之处在于:(1)当1个冲突环基中同时存在多条权值绝对值最小的边时, 增加新的去边规则, 进一步增强去边的确定性; (2)以单体型中0／1取值的概率确定孤立点的取值。
采用真实单体型进行试验, 利用CELSIM和MetaSim模拟测序生成器获得片段数据。试验结果显示, IHapCompass算法在不同的参数设置下均能获得较算法HapCompass、DGS和Fast Hare更高的单体型重建率, 且具有较高的执行效率, 是求解MWER模型的一种有效方法。
| [1] |
LANCIA G, BAFNA V, ISTRAIL S, et al. SNPs problems, complexity, and algorithms[C]//Proceedings of ESA-01. Aarhusv Denmark: Springer Press, 2001:182-193.
( 0)
|
| [2] |
LIPPERT R, SCHWARTZ R, LANCIA G, et al. Algorithmic strategies for the single nucleotide polymorphism haplotype assembly problem[J].
Briefings in Bioinformatics , 2002, 3 (1) : 23-31 DOI:10.1093/bib/3.1.23 (0)
|
| [3] |
SCHWARTZ R. Theory and algorithms for the haplotype assembly problem[J].
Communications in Information & Systems , 2010, 10 (1) : 23-38 (0)
|
| [4] |
GERACI F. A comparison of several algorithms for the single individual SNP haplotyping reconstruction problem[J].
Bioinformatics , 2010, 26 (18) : 2217-2225 DOI:10.1093/bioinformatics/btq411 (0)
|
| [5] |
AGUIAR D, ISTRAIL S. HapCompass: a fast cycle basis algorithm for accurate haplotype assembly of sequence data[J].
Journal of Computational Biology , 2012, 19 (6) : 577-590 DOI:10.1089/cmb.2012.0084 (0)
|
| [6] |
XIE Minzhu, WANG Jianxin, CHEN Xin. LGH: a fast and accurate algorithm for single individual haplotyping based on a two-locus linkage graph[J].
IEEE／ACM Transactions on Computational Biology and Bioinformatics , 2015, 12 (6) : 1255-1266 DOI:10.1109/TCBB.2015.2430352 (0)
|
| [7] |
RIZZI R, BAFNA V, ISTRAIL S, et al. Practical algorithms and fixed-parameter tractability for the single individual SNP haplotyping problem[C]／／Proceedings of WABI-02. Rome, Italy: Springer Press, 2002:29-43.
(0)
|
| [8] |
BAFNA V, ISTRAIL S, LANCIA G, et al. Polynomial and APX-hard cases of the individual haplotyping problem[J].
Theoretical Computer Science , 2005, 335 (1) : 109-125 DOI:10.1016/j.tcs.2004.12.017 (0)
|
| [9] |
XIE Minzhu, CHEN Jian′er, WANG Jianxin. Research on parameterized algorithms of the individual haplotyping problem[J].
Journal of Bioinformatics and Computational Biology , 2007, 5 (3) : 795-816 DOI:10.1142/S0219720007002710 (0)
|
| [10] |
CILIBRASI R, IERSEL L V, KELK S, et al. The complexity of the single individual SNP haplotyping problem[J].
Algorithmica , 2007, 49 (1) : 13-36 DOI:10.1007/s00453-007-0029-z (0)
|
| [11] |
WANG Ruisheng, WU Lingyun, LI Zhenping, et al. Haplotype reconstruction from SNP fragments by minimum error correction[J].
Bioinformatics , 2005, 21 (10) : 2456-2462 DOI:10.1093/bioinformatics/bti352 (0)
|
| [12] |
PANCONESI A, SOZIO M. Fast hare: a fast heuristic for single individual SNP haplotype reconstruction[C]／／Proceedings of WABI-04. Bergen, Norway: Springer Press, 2004:266-277.
(0)
|
| [13] |
GENOVESE L M, GERACI F, PELLEGRINI M. SpeedHap: an accurate heuristic for the single individual SNP haplotyping problem with many gaps, high reading error rate and low coverage[J].
IEEE／ACM Transactions on Computational Biology and Bioinformatics , 2008, 5 (4) : 492-502 DOI:10.1109/TCBB.2008.67 (0)
|
| [14] |
LEVY S, SUTTON G, NG P C, et al. The diploid genome sequence of an individual human[J].
PLOS Biology , 2007, 5 (10) : 2113-2144 (0)
|
| [15] |
BANSAL V, BAFNA V. HapCUT: an efficient and accurate algorithm for the haplotype assembly problem[J].
Bioinformatics , 2008, 24 (16) : i153-i159 DOI:10.1093/bioinformatics/btn298 (0)
|
| [16] |
MAZROUEE S, WANG W. FastHap: fast and accurate single individual haplotype reconstruction using fuzzy conflict graphs[J].
Bioinformatics , 2014, 30 (17) : i371-i378 DOI:10.1093/bioinformatics/btu442 (0)
|
| [17] |
吴璟莉.遗传多态性检测中组合优化问题的研究[D].长沙:中南大学, 2008.
WU Jingli. Research on the combinatorial optimization problem in detection of genetic diversities[D]. Changsha: Central South University, 2008. http://www.shangxueba.com/lunwen/v46573.html ( 0)
|
| [18] |
WU Jingli, LIANG Binbin. A fast and accurate algorithm for diploid individual haplotype reconstruction[J].
Journal of Bioinformatics and Computational Biology , 2013, 11 (4) : 1350010 DOI:10.1142/S0219720013500108 (0)
|
| [19] |
The The International HapMap Consortium. The international HapMap project[J].
Nature , 2003, 426 (6968) : 789-796 DOI:10.1038/nature02168 (0)
|
| [20] |
MYERS G. A dataset generator for whole genome shotgun sequencing[C]／／Proceedings of ISMB-99. Heidelberg, Germany: AAAI Press, 1999:202-210.
(0)
|
| [21] |
RICHTER D C, OTT F, AUCH A F, et al. MetaSim: a sequencing simulator for genomics and metagenomics[J].
Plos One , 2008, 3 (10) : 1-12 (0)
|