CN1169967C - 一种基于重复序列识别的全基因组测序数据的拼接方法 - Google Patents

Abstract

本发明涉及一种基于重复序列识别的全基因组测序数据的拼接方法，首先计算散弹法测序数据中非重复片段和重复片段出现的概率分布，并根据这一概率分布确定重复序列的识别标准，然后用该标准屏蔽重复序列，再根据目标基因组的大小进行分组拼接，将得到的大片段中的N恢复成原有碱基，并利用同一克隆正反向测序的信息找出相关的大片段以及可能出现在它们之间的读数，并将其连接起来，所有能连接的片段都连接以后，再使用正反向测序信息把大片段排好顺序，即得到目标基因组的工作框架图。本发明的方法，具有提高效率、能处理复杂基因组、明显减少错误拼接出现的概率、减少大量前期生物学实验准备等优点。

Description

一种基于重复序列识别的全基因组测序数据的拼接方法

技术领域

本发明涉及一种基于重复序列识别的全基因组测序数据的拼接方法，属基因工程技术领域。

背景技术

基因组学是对某一生物的全套遗传物质进行全面的分析，从整体的角度去了解遗传信息的功能和作用。其最重要的一步是将该生物的全套遗传信息测定出来，即知道该生物所有的核酸碱基排列顺序，即所谓的全基因组测序分析。目前全基因组测序主要采用两种策略：1、“分级克隆法”，即先将较大基因组打断成中等大小片段(150kb~300kb)并克隆，再将中等片段打断成小片段(1kb~3kb)进行测序，最后通过计算机进行数据拼接。如人类基因组计划(lander ES，2001)即采用此方法。该方法的优点是有较高的准确性，因为现有的计算机数据拼接软件对于较小基因组的拼接，如细菌、病毒等微生物，有较高的准确性。但是该方法首先必须对待测基因组有较多的了解，如已知相当数量的分子标记才能正确区分和定位中等大小片段克隆，进而正确进行小片段拼接。因此该方法需要大量前期的实验工作基础。对于以前了解不多的基因组，该方法具有明显的缺点。2、“霰弹法”，又称鸟枪法(shotgun sequencing)，是直接将大基因组随机地打断成小片段并测序，再用计算机对数据进行自动拼接组装。如美国塞莱拉(celera)公司进行人类基因组测序时即采用这种方法。这样做的好处是：对于较大基因组来说，便于进行规模化的操作，从而大大节省时间和减少人工操作，进而节约大量的经费支出，提高效率。但目前已有的用于大规模基因组数据拼接的计算机软件都是基于微生物基因组数据的拼接方法，存在明显的局限性，因为高等动植物的基因组与微生物基因组之间具有较大的差异，如高等生物基因组中具有大量的重复序列(即相同碱基顺序的片段在基因组中多次重复出现的现象)，而微生物中则没有。因此，用已有的数据拼接软件来进行高等生物较大基因组的数据拼接时，将会产生大量的错拼现象。事实上，塞莱拉公司最后只有在参考了人类基因组计划的实验数据后才完成他们的“工作框架图”。

发明内容

本发明的目的是提出一种基于重复序列识别的全基因组测序数据的拼接方法，在分析了高等动植物基因组中的重复序列在采用“霰弹法”进行测序时所具有的规律，解决由重复序列造成的高等动植物大基因组的数据拼接错误的问题，从而为采用“霰弹法”进行高效、快速的高等动植物全基因组测序分析提供可靠的手段。

本发明提出的基于重复序列识别的全基因组测序数据的拼接方法，包括以下步骤：

(1)设定一个最小的DNA片段长度为15bp-20bp，

计算散弹法测序数据中非重复片段出现的概率分布：

下列公式中各参数的含意为：G：基因组总长，L：测序平均有效读长，N：成功测序反应数，F：识别最小片段长度，

定义一个随机变量Y_ik描述用散弹法对全基因组测序中上述指定长度DNA片段出现K次的事件：

若某点开始的片段出现次数为k，则有k个测序片段的起点应在基因组上[I-L+F，i]区间内，而其它N-k个测序片段的起点不在此区间内，这一区间的长度为L-F，若所有测序片段起点在基因组上随机分布，则根据古典概型，上述随机变量等于1的概率为：

$P(Y_{\mathrm{ik}}=1)=C_N^k{\left(\frac{L-F+1}{G}\right)}^k{(1-\frac{L-F+1}{G})}^{N-k}---\left(1\right)$

一次测序中出现次数为k的片段的平均个数可表示为：

$E\left(Y_k\right)=E\left(\underset{i=1}{\overset{G}{\mathrm{\Σ}}}{Y}_{\mathrm{ik}}\right)=G\·C_N^k{\left(\frac{L-F+1}{G}\right)}^k{(1-\frac{L-F+1}{G})}^{N-k}---\left(2\right)$

使用下式作为一次测序中出现次数为k的片段出现概率的估计值；

P_k＝E(Y_k)/G                                                (3)

(2)计算重复片段出现的概率分布：

设片段为一个有m个拷贝的重复序列，出现在基因组中的m个不同位置，在散弹法测序数据集中的出现次数是所有位置出现次数的和，用G_mk表示有m个拷贝的重复序列在一次测序中的出现次数为k的概率，则上述关系可用数学公式表示为：

$G_{m0}=P_0^m$

$G_{m1}=C_m^1\·P_1\·P_0^{m-1}$

$G_{m2}=C_m^2\·P_1^2\·P_0^{m-2}+C_m^1\·P_2\·P_0^{m-1}$

$G_{m3}=C_m^3\·P_1^3\·P_0^{m-3}+C_m^2\·C_2^1\·P_1\·P_2\·P_0^{m-2}+C_m^1\·P_3\·P_0^{m-1}$

………

G_mj+＝1-G_m0-G_m1…-G_mj-1

其中G_mj+表示出现次数为j和更多的概率；

(3)重复序列的识别：

选取非重复片段出现概率最接近0.3％左右的次数为重复片段的判别标准，超过这一标准的片段就认为它属于重复序列，否则就是非重复序列；

(4)首先屏蔽重复序列，将上述散弹法测序数据中与识别出的片段相同的碱基改写为N，屏蔽后剩余长度超过一个定值的测序数据仍进入拼接过程；

(5)若目标基因组大小为1百万-3千万碱基，则屏蔽去重复序列后不分组直接进入拼接，若目标基因组明显大于上述范围，则需要按照测序读数之间的关连进行分组，例如可以将参加拼接的读数随机分为若干组，每组读数个数在5至10万个之间，每组数据进行初步拼接，对拼接得到的大片段进行比较，同源性高的聚为一组，把组成它们的读数重新拆出来放在一起，再次进行拼接；

(6)将得到的大片段中的N恢复成原有碱基，并利用同一克隆正反向测序的信息找出相关的大片段以及可能出现在它们之间的读数，并将其连接起来；

(7)所有能连接的片段都连接以后，再使用正反向测序信息把大片段排好顺序，即得到目标基因组的工作框架图。

本发明的基于重复序列识别的全基因组测序数据的拼接方法，具有提高效率、能处理复杂基因组、明显减少错误拼接出现的概率、减少大量前期生物学实验准备等优点。采用本方法进行了水稻基因组的拼接工作，结果显示本方法完全可以胜任水稻这样复杂基因组的拼接工作，在只进行了4.2倍基因组总长测序的情况下，拼接得到的大片段已经覆盖了基因组中90％以上的基因，拼接错误率在1％左右，这大致相当于已有技术中塞莱拉公司对果蝇基因组进行13倍基因组总长测序后得到的结果，果蝇基因组的重复序列实际上明显少于水稻，因此其拼接难度也大大小于水稻。另外，使用本发明的方法对1％人类基因组数据拼接结果表明，可以节省93％的计算机机时和84％的计算机内存空间。

附图说明

图1a和图1b为显示测序覆盖度为1X和4X时选择不同标准各种拷贝重复序列被选出的概率，图中每条线代表一种判断标准。

图2表示不同测序量各种拷贝数重复序列被选出的概率比较，标准的选择使拷贝1序列选出概率保持在0.3％左右。

图3显示各插入片段均取10X冗余度时对不同长度洞均可覆盖两次的概率。

具体实施方式

下面结合附图，详细介绍本发明方法的各个步骤：

为进行重复序列识别，本发明首先设定一个最小的片段长度，一般设为15bp-20bp，小于这一长度的重复序列将不再考虑。为简化模型，假设所有测序读长相等，均为L。

下列公式中参数意义：

G：基因组总长，L：测序平均有效读长N：成功测序反应数，F：识别最小片段长度。

计算霰弹法测序中非重复的小片断的出现次数：

定义一个随机变量Y_ik描述用散弹法对全基因组测序中上述指定长度DNA片段出现K次的事件：

若某点开始的片段出现次数为k，则有k个测序片段的起点应在基因组上[I-L+F，i]区间内，而其它N-k个测序片段的起点不在此区间内，这一区间的长度为L-F，若所有测序片段起点在基因组上随机分布，则根据古典概型，上述随机变量等于1的概率为：

$P(Y_{\mathrm{ik}}=1)=C_N^k{\left(\frac{L-F+1}{G}\right)}^k{(1-\frac{L-F+1}{G})}^{N-k}---\left(1\right)$

一次测序中出现次数为k的片段的平均个数可表示为：

$E\left(Y_k\right)=E\left(\underset{i=1}{\overset{G}{\mathrm{\Σ}}}{Y}_{\mathrm{ik}}\right)=G\·C_N^k{\left(\frac{L-F+1}{G}\right)}^k{(1-\frac{L-F+1}{G})}^{N-k}---\left(2\right)$

使用下式作为一次测序中出现次数为k的片段出现概率的估计值；

P_k＝E(Y_k)/G                                             (3)

在实际使用重复序列识别程序时，片段长度是根据基因组大小选择的。对于水稻基因组来说，其基因组大小大约是430Mb，小片段总数大约是10⁸数量级，因此选择20bp长的小片段。此时共有约10¹²种不同小片段，可保证不会由于随机因素而在基因组中出现相同的。如果考虑的是细菌基因组，小片段总数约为10⁶数量级，其长度可缩短为15bp，仍可保证不会由于随机而出现相同。

计算重复序列中指定长度片段的出现次数：

分析推导过程不难看出上述概率都是非重复序列的概率，因为实际上假设每个片段都只出现在基因组的一个地方。如果片段是有m个拷贝的重复序列，它将出现在基因组的m个不同位置。在霰弹法测序数据集中看到的出现次数将是所有这些位置测序出现次数的和。例如测序集中出现次数为0，意味着所有m个位置出现次数都要为0；测序集中出现次数为1，则只有一个位置覆盖为1，其它都要为0；而测序集中出现次数为2，则可能只有一个位置出现次数为2，其它都为0；或有两个位置为深度1，其它均为0；等等。用G_mk表示有m个拷贝的重复序列在一次测序中的出现次数为k的概率，则上述关系可用数学公式表示为

设片段为一个有m个拷贝的重复序列，出现在基因组中的m个不同位置，在散弹法测序数据集中的出现次数是所有位置出现次数的和，用G_mk表示有m个拷贝的重复序列在一次测序中的出现次数为k的概率，则上述关系可用数学公式表示为：

$G_{m0}=P_0^m$

$G_{m1}=C_m^1\·P_1\·P_0^{m-1}$

$G_{m2}=C_m^2\·P_1^2\·P_0^{m-2}+C_m^1\·P_2\·P_0^{m-1}$

$G_{m3}=C_m^3\·P_1^3\·P_0^{m-3}+C_m^2\·C_2^1\·P_1\·P_2\·P_0^{m-2}+C_m^1\·P_3\·P_0^{m-1}$

………

G_mj+＝1-G_m0-G_m1…-G_mj-1

其中G_mj+表示出现次数为j和更多的概率；

(4)识别重复序列

根据上述概率，算出特定测序条件下(指基因组长度、总测序量、平均读长等参数)非重复序列和不同拷贝数重复序列中片段表现为一定测序出现次数的概率。然后确定一个适当的出现次数标准，超过这一标准的片段就假定它属于重复序列，否则就是非重复序列。根据各G_mk值算出不同拷贝数重复序列在这一判断标准下漏掉的概率，如果已知基因组中各种拷贝数重复序列所占比例，则还可知道挑出的重复序列中各种拷贝数所占的比例。

为节省篇幅，只给出图1显示测序覆盖度为1X和4X时选择不同标准各种拷贝重复序列被选出的概率。实际工作中选择标准最重要的指标是拷贝数为1的序列(即非重复序列)被选出的概率。由于非重复序列在基因组中一般要占到总长的2/3或更多，它被选出的概率必须充分小，以便保证选出的绝大多数确实是重复序列。在不同测序覆盖度下，本发明确定这一概率保持在0.3％左右。表1.表示不同测序量下的重复序列识别标准。图2表示在这一标准下不同测序量各种拷贝数重复序列被选出的概率比较。从图2中可见测序量达到4X以上时进一步增加测序量对重复序列识别的改善已不太明显，此时拷贝数在5以上的重复序列基本上都可识别。

                表1.不同测序量下的重复序列识别标准

测序覆盖度	2X	4X	6X
				重复序列判断标	出现次数7以上	出现次数11以上	出现次数13以上

有了上述统计模型，就可以根据3-4X测序数据识别出绝大部分拷贝数在5以上的重复序列。这使我们可以建立一套测序流程，处理高等生物复杂基因组全基因组鸟枪法测序数据。该流程主要步骤如图3所示。

屏蔽重复序列的方法是把测序数据中与识别出的片段相同的碱基改写为N。屏蔽后剩余长度超过100bp的测序数据仍进入拼接过程。屏蔽掉重复序列后大大减低了拼接的复杂性，使phrap等常用软件可以处理更大的数据集。如果目标基因组大小在数百万到数千万碱基的范围，屏蔽掉重复序列后可不分组直接投入拼接。但若目标基因组有数亿或更多碱基，则需要按照测序读数之间的关连进行分组，然后再逐组进行拼接。拼接后可利用屏蔽前的测序数据恢复contig中的重复序列，并利用正反向测序信息进行拼接和恢复正确性的检验。由于分组不会完全合理，还需要用blast等软件对各contig进行比较以便去除冗余。然后进一步利用正反向测序信息恢复一些较长的重复序列，并构建各contig之间的顺序关系。这样就基本完成了工作框架图的拼接工作。

设计测序插入片段长度的分布：

由于在上述拼接过程中屏蔽掉了重复序列，这样就会在基因组序列中留下一些洞(gap)。为了正确地填补这些洞，必须首先构建正确的片段框架(scaffold)。在没有详细的物理图谱等额外信息的情况下，就要精心设计测序插入片段的长度分布，为保证所有屏蔽重复序列后留下来的洞都能被适当长度的插入片段克隆覆盖两次以上。这就需要对插入片段长度分布进行设计。

插入片段覆盖指定长度的洞的概率仍可用(3)式计算，只是为了构建框架片段需要做一些小的修改。主要是要在洞的两边各留下50bp以上的序列，以便进行匹配识别。因此(3)式变为：

$P_k=P(Y_{\mathrm{ik}}=1)=C_N^k{\left(\frac{L-F-100}{G}\right)}^k{(1-\frac{L-F-100}{G})}^{N-k}---\left(4\right)$

上式表示起点为i，长度为F的洞被长度为L的插入片段覆盖k次的概率。其中N为插入片段总数，G为基因组总长。若忽略基因组起始和结尾的区域，(4)式对任何点都成立。因此可略去下标i。

若要求洞被覆盖两次以上，概率为：

           P₂₊＝1-P₀-P₁                                (5)

由于当插入片段长度远大于洞长时对构建片段框架不利，而接近洞长时覆盖效率又太低，根据经验，本发明规定长度L的片段只用于覆盖0.2L至0.6L的洞。用人和水稻已有的序列进行检验，发现屏蔽重复序列后留下的最大洞在20-25kb左右，据此选定插入片段长度及其负责覆盖洞长的数值见表2。

               表2.插入片段长度及负责覆盖洞的长度

插入片段长度(kb)	3	8	20	50
					覆盖洞下限(kb)	0.6	1.6	4	10
覆盖洞上限(kb)	1.8	4.8	12	30

测序成功率一般在90％左右，因此即使对所有插入片段都进行两端测序，也仍会有20％左右克隆只能得到一端测序数据。所以在本发明的模型中规定20％测序片段没有克隆另一端测序片段的数据。

如果没有关于洞的长度分布的信息，本发明建议各长度插入片段的覆盖度均取为10X，这样对各种长度的洞覆盖两次的概率都在99％左右，就能取得较好的效果，见图4。如果能根据已有的部分序列得到洞的长度分布，则可从上述初值出发，用公式(4)和(5)得到覆盖各种长度洞的概率，再结合洞长分布可得到各长度洞未能被覆盖的期望值，从而可计算出在这种插入片段分布下单位基因组长度上会留下多少洞。以此为目标函数，以总测序量为约束条件，可用非线性规划的方法得到一定测序量下使遗留洞长达到最小的插入片段长度分布。一般情况下插入片段长度分布不需要这样高的精度，也可采用excel等表格软件辅助进行手动调整，其结果通常已可满足使用需要。

Claims (1)

1、一种基于重复序列识别的全基因组测序数据的拼接方法，其特征在于该方法包括以下步骤：

(1)设定一个最小的DNA片段长度为15bp-20bp，

计算散弹法测序数据中非重复片段出现的概率分布：

下列公式中各参数的含意为：G：基因组总长，L：测序平均有效读长，N：成功测序反应数，F：识别最小片段长度，

定义一个随机变量Y_ik描述用散弹法对全基因组测序中上述指定长度DNA片段出现K次的事件：

若某点开始的片段出现次数为k，则有k个测序片段的起点应在基因组上[I-L+F，i]区间内，而其它N-k个测序片段的起点不在此区间内，这一区间的长度为L-F，若所有测序片段起点在基因组上随机分布，则根据古典概型，上述随机变量等于1的概率为：

$P(Y_{\mathrm{ik}}=1)=C_N^k{\left(\frac{L-F+1}{G}\right)}^k{(1-\frac{L-F+1}{G})}^{N-k}---\left(1\right)$

一次测序中出现次数为k的片段的平均个数可表示为：

$E\left(Y_k\right)=E\left(\underset{i=1}{\overset{G}{\mathrm{\Σ}}}{Y}_{\mathrm{ik}}\right)=G\·C_N^k{\left(\frac{L-F+1}{G}\right)}^k{(1-\frac{L-F+1}{G})}^{N-k}---\left(2\right)$

使用下式作为一次测序中出现次数为k的片段出现概率的估计值；

P_k＝E(Y_k)/G                                              (3)

(2)计算重复片段出现概率分布：

设片段为一个有m个拷贝的重复序列，出现在基因组中的m个不同位置，在散弹法测序数据集中的出现次数是所有位置出现次数的和，用G_mk表示有m个拷贝的重复序列在一次测序中的出现次数为k的概率，则上述关系可用数学公式表示为：

$G_{m0}=P_0^m$

$G_{m1}=C_m^1\·P_1\·P_0^{m-1}$

$G_{m2}=C_m^2\·P_1^2\·P_0^{m-2}+C_m^1\·P_2\·P_0^{m-1}$

$G_{m3}=C_m^3\·P_1^3\·P_0^{m-3}+C_m^2\·C_2^1\·P_1\·P_2\·P_0^{m-2}+C_m^1\·P_3\·P_0^{m-1}$

………

G_mj+＝1-G_m0-G_m1…-G_mj-1

其中G_mj+表示出现次数为j和更多的概率；

(3)重复序列的识别：

选取非重复片段出现概率为0.3％的次数为重复片段的判别标准，超过这一标准的片段就认为它属于重复序列，否则就是非重复序列；

(4)首先屏蔽重复序列，将上述散弹法测序数据中与识别出的片段相同的碱基改写为N，屏蔽后剩余长度超过50bp的测序数据仍进入拼接过程；

(5)若目标基因组大小为1百万-3千万碱基，则屏蔽去重复序列后不分组直接进入拼接，若目标基因组大于上述范围，则需要按照测序读数之间的关连进行分组，然后进行拼接；

(6)将得到的大片段中的N恢复成原有碱基，并利用同一克隆正反向测序的信息找出相关的大片段以及出现在它们之间的读数，并将其连接起来；

(7)所有能连接的片段都连接以后，再使用正反向测序信息把大片段排好顺序，即得到目标基因组的工作框架图。

Images