CN105303068A - 一种基于参考基因组和从头组装相结合的二代测序数据组装方法 - Google Patents

Abstract

本发明涉及一种基于参考基因组和从头组装相结合的二代测序数据组装方法。将基于参考基因组组装和基因组从头组装这两种策略结合起来，来克服这两者的劣势，同时充分利用这两者的优势。首先，利用基于参考基因组的策略获得一个连续性和准确率较好的基因组序列。然后，利用从头组装的策略获得一个从头组装的基因组，该基因组对物种特异性序列的组装方面表现较好。最后，将这两个基因组整合到一起，产生一个在准确率、连续性和完整性都较好的基因组。

Description

一种基于参考基因组和从头组装相结合的二代测序数据组装方法

技术领域

本发明涉及存在参考基因组的情况下对第二代测序读段进行全基因组组装的方法。

背景技术

目前，根据有无参考基因组来区分，二代测序数据的组装策略主要两种，一种是基因组从头组装，另一种是基于参考基因组的组装。但这两者各有优缺点。

基因组从头组装是在不依赖参考基因组或其他基因组的情况下，完全基于测序读段进行的组装。目前，对于序列组装主要有三种算法。第一种是贪婪算法。这种算法只要序列间有重复(一致)序列，就寻找最大重复区，合并更多序列。这种算法简单，可以达到局部最优的结果，但难以实现全局最优。在实际中，由于准确率比较低，使用的这种算法的组装软件也比较少。常用的软件有SSAKE，TIGRAssembler等。第二种是OLC(Overlap-Layout-Consensus)算法。这种算法将小序列间的关系以图的形式呈现，每个节点代表一个序列片段，如果两个序列片段间有重叠就用重叠区来连接两个节点。当测序的读段较长时，该算法能得到很好的结果。在使用Sanger读段来组装时，该算法被普遍使用。常用的软件有CAP3，Phrap，Newbler等。第三种是基于deBruijin图的算法。这种算法不直接使用读段来组装基因组。首先，将读段分割成长度相同的kmer片段，相邻kmer只差一个碱基。然后构建deBruijin图，以kmer为图的节点，以相邻kmer的连接为边。最后，寻找可以包含deBruijin图中所有读段所对应的路径。在实际组装过程中，一般要有去掉图中短的分支、低覆盖度的边、合并泡等操作，以获得更加准确的简化图。这种算法不需要读段之间的两两比对，并且使用的内存与测序覆盖深度关系不大，是一个线性时间问题，具有更低的时间复杂度。现在，使用二代测序数据进行的基因组组装大部分都是使用这种算法。常用的软件有：AllPath-LG，SOAPdenovo，Velvet等。

基因组从头组装的策略不依赖参考基因组的优点是可以更好地组装出一些物种特有的基因组片段，并且已开发出许多工具和软件，这提高了计算的自动化，降低了对科研人员对计算机技术的要求，在科学研究中应用广泛，如大熊猫，甜橙等物种的基因组组装。但是，它也有明显的缺点，在基因组的复杂度高的区域，利用从头组装的策略很难组装出来。因而，对于复杂度高、重复序列高的大基因组而言，这种策略在基因组组装的连续性上表现较差，即N50或N90很小；在准确性上的表现也较差，即在重叠群(contig)和支架(scaffolds)存在很多的错误连接。

基于参考基因组的组装，这是一种比较组装的方法，通过比较和参考已知或类似的基因组来实现短序列的组装。基于参考基因组的方法可以用于组装重叠群阶段或组装支架阶段。当参考基因组用于组装重叠群阶段时，基于参考基因组组装的方法大致上有两种。一种是基于比对构建重叠群的方法，如AMOScmp，samtools，LOCAS。另一种是基于参考基因组来划分更小的区域，在局部进行组装从而降低组装的复杂度，提高组装质量。这种方法没有现成的软件，需要设计更适合自己研究的流程，但大致思想都是“化整为零”的局部组装。这种方法已应用到拟南芥的研究中。当参考基因组用在组装支架阶段时，一般是基于重叠群在参考基因组上的比对的位置信息来构建支架，有时需要额外信息的支持，如双末端读段，或物理图谱等。使用这种方法的软件有scaffol_builder，RACA，nucmer的比对排序算法等。

基于参考基因组组装的策略的优点是可以降低组装的复杂度，对计算资源的消耗较低。同时可以提高基因组组装的正确率和连续性。但其缺陷也非常明显。该方法受参考基因组质量的影响较大。当参考基因组质量较高时，利用参考基因组组装的质量也会较高；当参考基因组质量较低时，利用参考基因组组装的质量也偏低。同时这种策略很难将该物种一些特异性的片段组装出来，这会导致组装的完整性大打折扣。

在本研究中，我们将基于参考基因组组装和基因组从头组装这两种策略结合起来，来克服这两者的劣势，同时充分利用这两者的优势。具体而言，我们首先利用基于参考基因组的策略获得一个连续性和准确率较好的基因组序列。然后，利用从头组装的策略获得一个从头组装的基因组，该基因组对物种特异性序列的组装方面表现较好。最后，我们将这两个基因组整合到一起，产生一个在准确率、连续性和完整性都较好的基因组。

发明内容

本发明充分利用现在已有的软件和方法，巧妙设计合理的基因组组装流程，整合各软件和方法的优点，弥补各软件和方法的缺点，从而得到质量更好的基因组。利用参考基因组能获得连续性更好的组装结果，因而，本研究首先基于参考基因组来划分更小的区块，在局部进行组装从而降低组装的复杂度，从而提高组装的连续性和准确度。其次，由于在一些重复序列较多的区域，利用已有的工具无法组装出来。因此，我们整合组装获得的序列与参考基因组序列信息进行局部重组装，来获得连续性更好的序列。最后，由于基于基因组组装的基因组序列会受限于参考基因组的序列，从而导致该物种一些特异性的基因组片段无法组装出来。因此，最后我们将基于基因组组装获得的序列与利用二代测序数据从头组装的结果整合到一起，从而进一步提高基因组组装的连续性、准确性和完整性。

本发明涉及一种利用参考基因组获得基因组编码信息的方法，其特征在于，取待测样品，利用二代测序方法测序，获得未处理的编码数据；将测序得到的未处理的编码数据经过以下步骤处理：

(1)对测序读段数据进行预处理；

(2)使用短序列比对工具将处理过的读段比对到参考基因组，然后根据比对结果来统计基因组每个位点的覆盖度，基因组平均覆盖度和覆盖率；

(3)根据参考基因组位点的覆盖度，定义连续覆盖的区域为区块；将小区块合并成大区块，并且相邻大区块设置重叠区；

(4)对于大区块内的读段，进行重头组装，将组装结果的重叠群序列输出；

(5)使用序列组装工具进一步处理，去除重叠冗余的区段；

(6)再经过构建支架序列、填补空缺、校正组装错误步骤，获得基因组编码信息。

根据优选的技术方案，步骤(5)以后，进一步采用可变长度末端序列唯一匹配的方法来定位序列在参考基因组序列上的起始位置，所在序列称为查询序列，进行局部重组装，将局部重组装得到的序列与相邻序列连接。

根据优选的技术方案，1)使用NUCmer软件将所有查询序列比对到参考基因组上；2)筛选末端的比对序列，包括开始末端和结束末端。

根据优选的技术方案，采用可变长度末端序列唯一匹配的方法定位了重叠群在参考基因组序列上的起始位置后，提取一条重叠群的末端位置到相邻的另一条重叠群始端位置内加上下游各一定长度，优选100-500bp，更优选200bp，范围内的读段；使用CAP3组装读段，得到组装的序列；然后，将局部重组装得到的序列与相邻的两条重叠群比对，设定比对长度和相似性的阈值，用该序列来连接相邻的重叠群。

根据优选的技术方案，将可以组装出基因组特异性序列的从头组装获得的基因组序列与采用基于参考基因组组装的基因组序列整合起来，从而获得连续性好、完整性好的基因组序列。

根据优选的技术方案，步骤(6)所述的构建支架序列，填补空缺，校正组装错误步骤，指的是利用二代测序读段获得的重叠群序列，使用SSPACE和/或OPERA工具，来构建支架序列；使用gapCloser来填补gaps；使用iCORN软件校正基因组组装错误碱基，来获得准确率更高的基因组。

根据优选的技术方案，步骤(3)中，所述的将小区块合并成大区块，并且相邻大区块设置重叠区，指的是如果相邻区块之间的距离小于50-300bp之间的某一设定数值，则将其合并；如果有读段两端分别位于相邻的两个区块内部，则将这两个区块合并；对于长度小于20-200KB之间的某一设定值的区块，将它与其相邻的区块合并成大区块，直到大区块长度大于所设定的数值,并且相邻大区块设置长度为1-5KB之间的某一设定值的重叠区。

根据优选的技术方案，步骤(4)中所述的从头组装指的是利用SOAPdenovo或者ZORRO等软件进行局部从头组装，组装过程中循环尝试不同的参数，选取重叠群N50最大的参数，并将其组装结果重叠群序列输出。

根据优选的技术方案，步骤(5)所述的使用序列组装工具进一步处理，去除重叠冗余的区段，指的是使用AMOS工具对上一步得到的重叠群序列进一步组装，去除不同重叠群之间可能存在的重叠冗余的区段。

根据优选的技术方案，所述的步骤(3)中，划分区块和大区块的时候，可以使用下列方法中的一种：

(a)基于读段在基因组比对的“断点”来划分区域；

(b)基于读段聚类和读段在基因组比较信息来划分区域。

附图说明

图1基于参考基因组组装的步骤(步骤1-步骤6)。

图2利用可变长度末端序列唯一匹配定位重叠群序列示意图。

图3利用局部重组装连接重叠群序列示意图。

图4全基因组从头组装序列连接重叠群序列示意图。

具体实施方式

为了更清楚的阐述本发明的方法内容，现将本发明所涉及的方法更进一步的总结如下，所涉及的实验数据或者软件等，属于本领域的常规技术，其并不对本专利的保护范围造成限制。

一种利用参考基因组获得基因组编码信息的方法，主要包括以下步骤：

(1)对测序读段进行预处理，正确率高的读段有利于提高后续分析的准确性，尤其是基因组组装的准确度和连续性。测序读段进行预处理的步骤主要包括使用软件jellyfish9和quake10对测序错误碱基进行校正，使用cutadapt去除接头序列，使用软件Trimmomatic11去除读段末端质量较低的碱基以及长度过短的读段(图1中步骤1)。

(2)使用短序列比对工具BWA12，将处理过的读段比对到参考基因组。然后根据比对结果来统计基因组每个位点的覆盖度，基因组覆盖度和覆盖率(图1中步骤2)。

(3)根据参考基因组位点的覆盖度，定义连续覆盖的区域为区块；如果相邻区块之间的距离小于100bp，则将其合并；如果有读段两端分别位于相邻的两个区块内部，则将这两个区块合并(图1中步骤3)。

(4)对于长度小于40KB的区块，将它与其相邻的区块合并成大区块，直到大区块长度大于40KB,并且相邻大区块设置长度为2KB的重叠区(图1中步骤4)。

(5)提取比对到每个大区块内的所有读段。对每个大区块内读段，使用SOAPdenovo进行局部从头组装，组装过程中循环尝试不同的参数，选取重叠群N50最大的参数，并将其组装结果重叠群序列输出(图1中步骤5)。

(6)使用AMOS工具对上一步得到的重叠群序列进一步组装，去除不同重叠群之间可能存在的重叠冗余的区段。组装软件在组装上的表现并不完美。在某些区域，尤其是一些短简单重复序列区域，组装软件不能完整地将该区域组装出来(图1中步骤6)。

(7)为了解决(6)中的问题，我们采用了基于参考基因组进行局部重组装的策略。采用可变长度末端序列唯一匹配的方法(图2)来定位其在参考基因组序列上的起始位置。具体过程及其技术要点为：1)使用NUCmer将所有查询序列比对到参考基因组上；2)筛选末端(包括开始末端和结束末端)比对序列。其中，末端序列的长度是变化的，这与查询序列最末端序列比对的长度相关。末端序列的筛选主要采用一下条件：末端序列离查询序列末端的距离；末端序列的比对长度；末端序列比对的相似度；末端序列是否为唯一比对。

采用可变长度末端序列唯一匹配的方法定位了重叠群在参考基因组序列上的起始位置后，提取一条重叠群的末端位置到相邻的另一条重叠群始端位置内加上下游各200bp范围内的读段。使用CAP3组装读段，得到组装的序列。然后，将局部重组装得到的序列与相邻的两条重叠群比对，设定比对长度和相似性的阈值，用该序列来连接相邻的重叠群(图3)。

至此，可以得到基于参考基因组组装的基因组序列。采用这种策略获得的基因组序列比从头组装的基因组在组装的连续性和完整性上要更好些，但在一些物种特异性的序列区域，采用这种方法一般无法将其组装出来。

(8)为了解决基于参考基因组无法组装出基因组特异性序列的缺陷，我们将可以组装出基因组特异性序列的从头组装获得的基因组序列与上面采用基于参考基因组组装的基因组序列整合起来，从而获得连续性好、完整性好的基因组序列。使用SOAPdenovo做全基因组从头组装。使用blat工具将重叠群序列比对从头组装的重叠群基因组序列。解析比对结果，确定相邻重叠群之间的gap序列(图4)。解析比对结果时，主要采用以下的判断条件：比对的长度，比对的相似性，末端未比对序列的长度，从头组装序列上用于连接的序列的覆盖度等。

(9)利用二代测序读段获得到的重叠群序列，使用SSPACE和OPERA工具，来构建支架序列。

(10)利用二代测序读段和得到的支架序列，使用gapCloser来填补gaps。

(11)使用iCORN软件校正基因组组装错误碱基(替换、缺失、插入错误)，来获得准确率更高的基因组。

为了进一步的验证本方法的实用性，研究人员采用了水稻叶片基因组进行测序，并进行了结果表征。

实施例1

本研究中，所用的测序数据为华中农业大学水稻研究团队提供的MH63和ZS97的叶片全基因组shotgun测序数据(籼稻品种明恢63和珍汕97)。取样时期为水稻三叶期，测序平台IlluminaHiseq2000，PE100测序，共有300bp，5k，10k三个不同insert-size的文库(表1)。此外，还使用了日本晴基因组参考基因组IRGSP-1.0(http://rapdb.dna.affrc.go.jp/)

表1测序数据统计

^a根据读段比对到日本晴基因组统计。

^b根据日本晴基因组大小估计。.

我们采用了基于日本晴参考基因组进行组装的策略。我们基于日本晴序列划分出多个区域，在每个区域内进行局部的从头组装。再先后利用日本晴序列，全基因组从头组装的序列和BES序列将重叠群延长。然后，利用300bp，5k，10k文库和BES数据来搭建支架。最后，基于支架在日本晴上的比对位置构建了假染色体序列。最终，得到组装的基因组序列。具体步骤如下：

(1)测序读段预处理。使用软件jellyfish和quake对测序错误碱基进行校正，使用cutadapt去除接头序列，使用软件Trimmomatic去除读段末端质量较低的碱基以及长度过短的读段。

(2)比对。使用BWA(参数：-q15)将处理好后的读段比对到日本晴基因组序列上，同时计算基因组每个位点的覆盖度，基因组覆盖度和覆盖率，统计比对结果。

(3)确定区块。根据参考基因组位点的覆盖度计算结果，定义连续覆盖的区域为区块；如果相邻区块之间的距离小于100bp，则将其合并；如果有读段两端分别位于相邻的两个区块内部，则将这两个区块合并。

(4)确定大区块。对于长度小于40KB的区块，将它与其相邻的区块合并成大区块，直到大区块长度大于40KB,并且相邻大区块设置长度为2KB的重叠区。

(5)收集大区块内的读段。收集比对到每个大区块内的所有读段，包括短片段文库和两个长片段文库的读段。

(6)局部从头组装。对每个大区块内读段，使用SOAPdenovo进行从头组装，组装过程中循环尝试不同的参数(参数：SOAPdenovo-63merall–K[29|31|..|63]-R)，选取重叠群N50最大的参数，并将其组装结果重叠群序列调出。

(7)使用AMOS工具对上一步得到的重叠群序列进一步组装，去除不同重叠群之间可能存在的重叠冗余的区段。

(8)Gap区重新组装。使用NUCmer比对，定位重叠群在reference上的坐标，确定相邻重叠群之间的Gap区域和该区域的比对覆盖范围和覆盖深度；提取gap区加上下游(gap读段Extend)各200bp范围内比对到参考基因组上的读段(读段长度和质量控制)。使用Cap3组装读段,得到重叠群序列；使用blat比对重叠群序列到gap区域两侧各300bp(gapRefExtend)范围内的相邻重叠群序列。

(9)整合abinitio组装结果。使用SOAPdenovo做abinitio全基因组组装。将第(6)、(8)步得到的重叠群合并，用blat工具比对abinitio组装的支架到第(8)步得到的重叠群。解析比对结果，确定相邻重叠群之间的gap序列。

(10)整合BES数据。比对BES序列到第(9)步得到的重叠群，将可以连接到一起的重叠群连接到一起。

(11)校正碱基错误。使用iCORN软件校正第(10)步得到的重叠群长度<＝3bp的组装错误碱基(替换、缺失、插入错误)。

(12)构建支架。利用三个片段文库的双端序列，对每条染色体使用SSPACE进行scaffolding操作；紧接着使用Opera利用5KB,10KB文库的双端序列进一步做scaffolding；最后使用SSPACE利用BES双端数据做scaffolding，最终得到每条染色体上的scaffold序列。

(13)填补gaps。利用三个片段文库的双端序列对第(12)步得到的所有的支架，使用gapCloser来填补gaps。

基于日本晴参考基因组局部组装得到重叠群。然后先后利用日本晴序列、全基因组从头组装的序列和BES序列将重叠群延伸。重叠群延伸之后，MH63组装得到了28,136条长度大于或等于500bp的重叠群序列，N50为107,523bp，序列总长度351,756,295bp(表2)；ZS97组装得到了22,041条长度大于或等于500bp的重叠群序列，N50为188,515bp，序列总长度350,033,712bp(表3)。连续性好的重叠群将有利于构建支架，同时也有利于后续更加容易、更加准确地对其进行基因组注释和比较基因组学分析。

在组装得到重叠群之后，我们进一步用pair-end读段和mate-pair读段，使用SSPACE和OPERA来搭建支架，然后再使用SSPACE软件。在得到支架后，使用gapCloser来填补gaps，最终得到支架序列。MH63基因组组装得到了11,534条支架序列(>＝500bp)，N50为514,605bp，序列总长度398,308,932bp(表4)；ZS97基因组组装得到了8,501条支架序列(>＝500bp)，N50为910,477bp，序列总长度386,061,070bp(表5)。

表2MH63基因组序列重叠群组装结果统计^a

a只对长度大于或等于500bp序列进行统计。

b长度大于或等于N50序列的总长大于或等于组装总长一半时，这些序列中最短的长度。

c长度大于或等于N50的序列的个数。

表3ZS97基因组序列重叠群组装结果统计^a

a，b，c同表2。

表4MH63基因组序列支架组装结果统计^a

a，b，c同表2。

表5ZS97基因组序列支架组装结果统计^a

a，b，c同表2。

我们将基于参考基因组得到的序列和全基因组从头组装得到的序列、日本晴序列、BES序列整合到一起，相互补充不足。将我们获得的基因组和直接SOAPdenovo组装的结果比对，可见我们采用的方法在基因组组装上有非常大的提高(表6)。这是首次明确地综合利用多种相关序列信息来组装二代数据，而不是局限于仅仅利用测序的读段来组装。这为他人的研究提供了很好的方法和思路。

表6不同基因组组装方法组装结果比较^a

^a对长度大于或等于500bp的序列进行统计。

^b长度大于或等于N50序列的总长大于或等于组装总长一半时，这些序列中最短的长度。

本发明中所涉及的软件为商业化的常规软件，本领域技术人员在了解其作用的基础上，可以使用其它软件或者方法替代。其并不会影响本发明的方法的实施效果。

Claims (10)

1.一种利用参考基因组获得基因组编码信息的方法，其特征在于，取待测样品，利用二代测序方法测序，获得未处理的编码数据；将测序得到的未处理的编码数据经过以下步骤处理：

(1)对测序读段数据进行预处理；

(2)使用短序列比对工具将处理过的读段比对到参考基因组，然后根据比对结果来统计基因组每个位点的覆盖度，基因组平均覆盖度和覆盖率；

(3)根据参考基因组位点的覆盖度，定义连续覆盖的区域为区块；将小区块合并成大区块，并且相邻大区块设置重叠区；

(4)对于大区块内的读段，进行重头组装，将组装结果的重叠群序列输出；

(5)使用序列组装工具进一步处理，去除重叠冗余的区段；

(6)再经过构建支架序列、填补空缺、校正组装错误步骤，获得基因组编码信息。

2.根据权利要求1所述的利用参考基因组获得基因组编码信息的方法，其特征在于步骤(5)以后，进一步采用可变长度末端序列唯一匹配的方法来定位序列在参考基因组序列上的起始位置，所在序列称为查询序列，进行局部重组装，将局部重组装得到的序列与相邻序列连接。

3.根据权利要求2所述的利用参考基因组获得基因组编码信息的方法，其特征在于1)使用NUCmer软件将所有查询序列比对到参考基因组上；2)筛选末端的比对序列，包括开始末端和结束末端。

4.根据权利要求2到3所述的利用参考基因组获得基因组编码信息的方法，其特征在于，采用可变长度末端序列唯一匹配的方法定位了重叠群在参考基因组序列上的起始位置后，提取一条重叠群的末端位置到相邻的另一条重叠群始端位置内加上下游各一定长度，优选100-500bp，更优选200bp范围内的读段；使用CAP3组装读段，得到组装的序列；然后，将局部重组装得到的序列与相邻的两条重叠群比对，设定比对长度和相似性的阈值，用该序列来连接相邻的重叠群。

5.根据权利要求2到4任一项所述的利用参考基因组获得基因组编码信息的方法，其特征在于将可以组装出基因组特异性序列的从头组装获得的基因组序列与采用基于参考基因组组装的基因组序列整合起来，从而获得连续性好、完整性好的基因组序列。

6.根据权利要求1-5任一项所述的利用参考基因组获得基因组编码信息的方法，其特征在于步骤(6)所述的构建支架序列，填补空缺，校正组装错误步骤，指的是利用二代测序读段获得的重叠群序列，使用SSPACE和/或OPERA工具，来构建支架序列；使用gapCloser来填补空缺；使用iCORN软件校正基因组组装错误碱基，来获得准确率更高的基因组。

7.根据权利要求1-6任一项所述的利用参考基因组获得基因组编码信息的方法，其特征在于步骤(3)中，所述的将小区块合并成大区块，并且相邻大区块设置重叠区，指的是如果相邻区块之间的距离小于50-300bp之间的某一设定数值，则将其合并；如果有读段两端分别位于相邻的两个区块内部，则将这两个区块合并；对于长度小于20-200KB之间的某一设定值的区块，将它与其相邻的区块合并成大区块，直到大区块长度大于所设定的数值,并且相邻大区块设置长度为1-5KB之间的某一设定值的重叠区。

8.根据权利要求1-7任一项所述的利用参考基因组获得基因组编码信息的方法，其特征在于步骤(4)中所述的从头组装指的是利用SOAPdenovo或者ZORRO软件进行局部从头组装，组装过程中循环尝试不同的参数，选取重叠群N50最大的参数，并将其组装结果重叠群序列输出。

9.根据权利要求1-8任一项所述的利用参考基因组获得基因组编码信息的方法，其特征在于步骤(5)所述的使用序列组装工具进一步处理，去除重叠冗余的区段，指的是使用AMOS工具对上一步得到的重叠群序列进一步组装，去除不同重叠群之间可能存在的重叠冗余的区段。

10.根据权利要求1-9任一项所述的利用参考基因组获得基因组编码信息的方法，其特征在于，所述的步骤(3)中，划分区块和大区块的时候，可以使用下列方法中的一种：

(a)基于读段在基因组比对的“断点”来划分区域；

(b)基于读段聚类和读段在基因组比较信息来划分区域。