CN103388025B - 基于克隆dna混合池的全基因组测序方法 - Google Patents

Abstract

本发明公开了一种基于克隆DNA混合池的全基因组测序方法。该方法首先构建BAC文库并构建BAC克隆DNA混合池，对混合池DNA进行双末端的高通量测序；解析克隆的特征序列集合与k-mer集合，利用克隆的特征序列构建克隆重叠群，利用克隆重叠群将克隆k-mer集合分割成小的k-mer集合，将混合池的NGS序列进行组装，并将组装后的序列定位到克隆重叠群上；将组装后的序列与NGS序列进行双末端比对，根据比对结果连接定位后的序列，并确定它们的方向；最后利用分子标记确定克隆重叠群的相对位置与方向，将克隆重叠群的序列连接成染色体序列。本发明利用NGS高通量测序技术与克隆DNA混合池，不仅构建出全基因组的物理图谱，同时将拼装后的序列定位到物理图谱上，实现两者的整合。

Description

基于克隆DNA混合池的全基因组测序方法

技术领域

本发明涉及全基因组测序技术领域，具体地指一种基于克隆DNA混合池的全基因组测序方法。

背景技术

全基因组序列作为一种研究物种全部遗传信息的重要资源受到越来越多的关注。从2001年人类基因组草图的完成，到2003年基因组全序列的正式完成以来，越来越多的物种完成了全基因组测序。对于植物基因组，拟南芥(Arabidopsis thaliana)作为一种重要的模式植物，是第一个完成全基因组测序的高等植物。水稻作为重要的粮食作物，其全基因组序列也已完成。其中，水稻籼稻品种(Oryza sativa ssp. indica)和粳稻品种(Oryza sativa ssp. japonica)的全基因组序列均完成于2005年。另外，重要的饲料作物玉米的全基因组序列于2009年完成。

目前，全基因组测序的策略主要有两种：clone-by-clone(CBC)和全基因组鸟枪法。

CBC的策略是首先构建遗传图谱(Genetic map)和物理图谱(Physical map)，再根据物理图谱选择最小重叠路径(Minimal tiling path, MTP)上的克隆(BAC或YAC克隆)进行测序。对克隆序列组装后，再将克隆序列一步一步地回归定位到物理图谱、遗传图谱和染色体上。在整个全基因组序列构建过程中，构建物理图谱是非常耗时费力的一步。基于文库的构建物理图的常用方法大致可以分为两类，基于标记的方法和限制性酶切图谱的方法。基于标记的方法一般是将克隆混合成多维的混合池(Pool)，然后用特定的探针与混合池进行杂交或用设计好的引物对混合池进行PCR扫描，从而确定每个混合池含有的标记，再通过计算筛选得到每个克隆可能含有的标记。利用限制性酶切图谱的方法，一般是用一种或多种限制性酶将每个克隆进行完全消化，然后对消化后的DNA片段进行电泳分离，从而得到每个克隆消化后DNA片段的长度信息，再根据这些片段长度信息利用计算机软件把相互重叠的克隆连接起来，形成重叠群(Contig)。与基于标记的方法相比，限制性酶切图谱的方法具有更高的分辨率和准确性，是现在用来构建物理图谱应用较多方法。特别是利用荧光标记分辨DNA片段，并利用基因分析仪进行高通量自动化分析方法的开发，使得限制性酶切图谱的方法得到更广泛的应用。

物理图谱对于定位基因，确定染色体，构建框架和全基因组的序列组装是一种非常重要的工具，也是一种非常有用的基因组资源。同时，对于研究基因的功能和物种进化有着很大的价值。利用CBC策略完成全基因组测序的代表物种主要有人类，拟南芥，水稻的粳稻品种，玉米等。这些基因组的全序列的质量非常高，常作为很多研究中的参考序列。基于CBC策略进行全基因组测序，因为其高分辨率和高准确性，使其已经成为一种全基因组测序的“黄金法则”。但是，它也有不可避免的缺点：构建物理图谱需要耗费大量的人力物力和时间，是整个测序过程中最为复杂的一步。因此，这种全基因组测序策略的应用也受到一定限制。

另一种全基因组测序的策略是全基因组鸟枪法测序(Whole genome shotgun sequencing, WGS)。这种方法不需要构建物理图谱，直接将要测序物种的全基因组随机打断成一定长度范围的小片段，将这些小片断插入到载体中，再对这些小片段文库直接进行测序，如人类的全基因组测序时，除了CBC方法外，也用了全基因组鸟枪法；水稻的籼稻品种93-11的全基因组序列便是用这种方法完成的。然后，将得到的序列根据序列的重叠信息进行拼装，得到序列重叠群。再将这些拼装后的序列定位到染色体上。在对小片段文库进行测序时，又有两种策略，一种是进行单向测序即只测得插入片段一端的序列；另一种是进行双末端测序(Paired-end sequencing)，其两端的序列方向相反。相比于单向测序，双末端测序不仅能够利用序列之间的重叠信息来组装序列重叠群，也可以利用插入片段的长度信息来将序列重叠群连接成更长的Scaffold。尽管双末端测序能够将序列重叠群相连接，但是由于重复序列或大的缺口的存在，还有很多的序列无法连接起来。针对这个问题，常常会构建多个不同插入片段大小的测序文库如人类基因组测序时就构建了2-kb，10-kb和50-kb的测序文库，这样能够有效地提高序列拼装的连续性。由于全基因组鸟枪法测序有着简单快速并不需要构建物理图谱的优点，特别随着下一代测序技术(Next-generation sequencing technologies, NGS)的不断发展，这种方法的运用越来越广泛。

鸟枪法测序与下一代测序在本质上是一样的，只是相对于鸟枪法测序，下一代测序技术具有更高的通量。下一代测序技术的测序主要有三种：Roche/454 FLX(http://www.454.com/)，Illumina/Solexa Genome Analyzer (http://www.illumina.com)和Applied Biosystems SOLiD^TM System (https://www.appliedbiosystems.com)。随着下一代测序技术不断发展，它不仅应用于基因组从头测序，越来越多的其它应用技术也逐渐被开发出来，如对已存在参考序列的物种进行基因组重测序(Genome re-sequencing)以研究不同个体基因组结构间的差异；转录组测序(Transcriptome sequencing)用来分析特定时期的表达，发现可变剪切位点及等位基因特异表达；以及染色质免疫共沉淀(Chromatin immunoprecipitation, ChIP)测序用来分析DNA和蛋白质的相互作用；还有目标序列捕获测序(Targeted re-sequencing)用于检测易被遗漏的低水平变异。由于下一代测序技术测序通量高，并且成本相对较低，现在已经有很多物种的全基因组序列是利用下一代测序技术进行的。尽管下一代测序技术有很多优点，但是也存在不可忽视的缺点。这些缺点主要表现在：序列读长(Read length)太短，序列组装过程计算量很大，完成序列组装困难，组装后的序列包含很多缺口并难以填补，在没有参考序列的情况下很难将得到的scaffold定位到染色体上并确定它们之间的相对位置，特别是当基因组含有大量的重复序列或含有很大的基因家族以及大片段的重复时，这些缺点就更加突出。

两种全基因组测序的策略各有其优缺点，为了充分利用两种策略的优点，避免其缺点，在两种策略的基础上，开发了很多其它方法。如运用鸟枪法测序的基于序列的物理图谱的构建方法，其中Milosavljevic等于2005年发表在《Genome Research》上的基于短序列标签的基因组混合池索引(Short-tag pooled genomic indexing, ST-PGI)的方法，是利用鸟枪法对BAC克隆的混合池进行测序，再结合参考序列，将序列分配到对应的克隆上，并将序列和克隆定位到染色体上；Volik等于2003年发表在《Proceedings of the National Academy of Sciences of the United States of America》上的末端序列分析(End-sequence profiling, PS)的方法，是利用BAC末端序列将BAC克隆定位到参考序列上。比较以上两种方法，ST-PGI利用了混合池，对所有的克隆都进行了测序，而EPS只对克隆进行末端测序，ST-PGI获得的克隆序列信息更多，而EPS相对较少导致很多克隆无法定位。如果所分析物种的基因组序列较复杂，含有大量的重复序列，或参考序列的质量不高，就会导致很多序列无法定位或定位错误，特别是只利用末端序列信息EPS方法更容易出现以上问题。

ST-PGI和ESP两种方法能够运用的一个共同的前提，就是必须有一个较高质量的相同或相近物种的全基因组序列作为参考序列，它们只能用来定位序列和克隆。但是现实是，并不是每一个物种都有一个高质量的全基因组序列可以参考。如果要对没有参考序列的物种进行全基因组测序并得到高质量的序列作为其它物种的参考序列，CBC的策略还是首选。而CBC策略需要解决的一个问题是如何快速地建立一个物理图谱，针对这个问题, Oeveren等在2011年开发了全基因组剖析(whole genome profiling ,WGP)的方法，该方法发表在《Genome Research》上。WGP是将BAC克隆混合成二维(Tow dimension,2D)的混合池，然后对从每一个混合池得到的DNA用一种限制性内切酶进行酶切，再利用下一代测序技术对酶切的片段进行测序，得到各个酶切片段的前11-26个碱基作为标签，将各个混合池包含的标签集合进行解析得到各个克隆所有包含的标签集合，最后将克隆的标签集合转化为物理图谱拼装软件FPC(FingerPrinted Contig,V9.4)兼容的数据进行拼装，得到克隆重叠群。WGP充分利用了下一代测序技术的优势，实现了物理作图的更高通量。同时利用酶切位点后的序列作为标签，相对于传统的酶切谱作图的方法，序列标签的数量更多，准确性更高。WGP利用下一代测序技术用来构建物理图谱，无论在思想上还是在技术上都是一个很大的进步。但是这种方法还不完美，它的主要目标还只局限在构建物理图谱，标签的数量受酶切位点的多少和分布以及酶切效果影响。如果能够在构建物理图谱的同时将序列也拼装起来并进行定位，将会在保证基因组序列的准确性的情况下，大大减少进行全基因组测序的费用，并提高全基因组测序的速度，这正是该方法开发的主要目的。

发明内容

本发明正是根据以上全基因组测序策略，针对它们的缺点提出的一种基于克隆DNA混合池的全基因组测序方法，快速高效精确完成全基因组测序的同时，有效地降低成本。包括以下步骤：

1）提取全基因组DNA，构建BAC文库；

2）构建BAC克隆混合池；

3）提取BAC克隆混合池的DNA；

4）对步骤3）中的BAC克隆混合池的DNA利用NGS进行双末端测序；

5）扫描各个混合池的序列，获得各个混合池的特征序列集合与k-mer集合；

6）根据混合池的特征序列集合与k-mer集合，解析出各个克隆的特征序列集合与k-mer集合；

7）利用克隆的特征序列集合构建克隆重叠群；

8）利用克隆重叠群将克隆的k-mer集合分割成小的k-mer集合并定位到克隆重叠群上；

9）对步骤3）中的混合池的NGS序列进行拼装得到序列重叠群；

10）将序列重叠群定位到克隆重叠群上，并利用测序的双末端信息连接序列重叠群,确定它们的方向，得到克隆重叠群的序列；

11）利用分子标记确定克隆重叠群的相对位置与方向，将克隆重叠群的序列连接成整条染色体序列，得到全基因组序列。

其中，关于克隆特征序列集合与k-mer集合的解析的总算法如下：

某一物种BAC文库的克隆总数为a，构建的混合池总数为x，则：

第κ维，索引为λ混合池的k-mer集合表示为：P_{( κ , λ )}

包含某一给定克隆的混合池的k-mer集合为：{P_{( δ )}|δ<m,P_{( δ )}∈P}

包含同一克隆的混合池的k-mer集合的交集，即克隆的IKS为：

某一克隆在混合池中的排除并集EUKS为：

某一克隆k-mer集合的所有排除并集的交集为：

某一克隆的最终k-mer集合FKS为：CF=C_τ-CI_τ。

更详细地，关于k-mer集合解析更具体的算法如下：由于特征序列集合为k-mer集合的子集，以下算法完全适用于特征序列集合。

Ⅰ.集合基本运算

对于两个集合A和B，在该方法中涉及到的集合的基本运算有：

1.A∩B：其结果集为既在集合A中的元素，同时也在集合B中的元素组成的集合，即A与B的交集；

2.A∪B：其结果集为在集合A中的元素，或在集合B中的元素组成的集合，即A与B的并集；

3.A – B：其结果集为在集合A中的元素，但是不在集合B中的元素组成的集合，即A与B的差集；

4. :表示A_m∪A_m+1∪A_m+2…A_n-2∪A_n-1∪A_n

5. :表示A_m∩A_m+1∩A_m+2…A_n-2∩A_n-1∩A_n

以上定义的集合中，不含有重复元素，即集合中的每一个元素只能出现一次。但是在计算过程中还需要考察元素出现的次数，对于这种同时记录元素出现次数的集合，定义为频数集合(Frequency set)或键值集合(Key-value set, KV-set)。对于频数集合A和B，在该方法中有以下基本运算：

1.A+B：表示频数集合A与频数集合B中的元素求并集，如果某一元素即在A中出现，也在B中出现，其频数为在A中的频数与在B中频数的和。

2.A–B：表示在频数集合A中的元素但不在B中的元素，如果某一元素即在A中，也在B中，其频数为在A中的频数减去在B中频数后的差；如果这个差小于等于0，则将这个元素在结果频数集合中去除。

3.Key(A)：表示频数集合A中所有元素的集合。

Ⅱ.克隆特征序列集合与k-mer集合的解析

将混合池进行测序后，扫描混合池序列，得到各个混合池的特征序列（定义参见实施例1中“扫描测序结果并得到混合池的特征序列集合”）集合和k-mer集合(定义参见实施例1中“扫描测序结果并得到混合池的k-mer集合”)，再计算出各个克隆的特征序列和k-mer集合。其中计算克隆的特征序列和k-mer集合，是整个方法的核心部分。

以计算克隆的k-mer集合为例，假设将所有的克隆放入一个立方体后，构建克隆混合池，混合池的维数为x(x>0)；每一维的混合池的数目为n(n>0)，克隆总数为a（a>0）。对于第κ维，索引为λ混合池的k-mer集合表示为 P_{( κ , λ )} (κ<x, λ<n)，所有的混合池的k-mer集合组合后的集合为这个基因组序列的k-mer集合的全集，表示为，即所有的混合池的k-mer集合的并集。

对这x个混合池，它们包含有相同数目的克隆数，假设一个给定的混合池的k-mer集合为{P_{( δ )}|δ<m,P_{( δ )}∈P}，则由多个混合池共同确定的一个克隆的k-mer集合为，称之为克隆的初始k-mer集合(Initial k-mer set，IKS)，C_τ (0<τ<a)表示索引为τ的克隆的k-mer的IKS。

对于在指定混合池中的一个给定克隆，计算出这个混合池中除了给定克隆的其它所有克隆的k-mer的IKS的并集，CU(κ,τ)表示第κ维且包含索引为τ的克隆的混合池中，索引为τ的克隆的排除并集。由此得到的并集称为特定克隆排除并集(Excluded union k-mer set, EUKS)。再对某一个克隆的所有排除并集求交集，最后得出克隆的最终的k-mer集合(Final k-mer set,FKS)CF=C_τ-CI_τ。

根据特征序列与k-mer的定义可以知道，当它们的长度相同时，特征序列的集合是k-mer集合的子集。所以，克隆最终的特征序列的集合(Final feature sequence set, FFS)的计算方法与求最终k-mer集合的方法一样。如图8所示的同一个立方体中的8个克隆，这是克隆数最简单的情况。

在这里考虑最简单的情况，假设有8个克隆A1-A8(克隆与克隆真实的k-mer集合用同一符号表示，如A1即表示克隆A1，也表示A1真实的k-mer集合)，将它们放到一个立方体中(图8)。然后构建3维的克隆混合池(即混合池构建策略中PX，PY，PZ)，可以得到6个混合池P1-P6。根据克隆在立方体中的位置，可以知道每个混合池中包含的克隆，如下所列：

P1={A1,A2,A5,A6}P2={A1,A2,A3,A4}P3={A3,A4,A7,A8}

P4={A5,A6,A7,A8}P5={A1,A3,A5,A7}P6={A2,A4,A6,A8}

混合池P1包含有克隆A1，A2，A5和A6，其它的混合池与之类似。

假设现要求出克隆A1的最终的k-mer集合，则按照以下步骤进行：

1、对所有的混合池进行高通量测序，扫描测序结果，得到混合池P1-P6对应的k-mer集合B1-B6.

2、根据克隆所在的混合池的信息，对混合池进行求交集，如克隆A1在P1，P2和P5中，则A1的k-mer的交的集合为A1=B1∩B2∩B5，其它克隆的k-mer的交的集合与之类似，得到各个克隆的k-mer的交的集合：

A1=B1∩B2∩B5 A2=B1∩B2∩B6 A3=B2∩B3∩B5 A4=B2∩B3∩B6

A5=B1∩B4∩B5 A6=B1∩B4∩B6 A7=B3∩B4∩B5 A8=B3∩B4∩B6

3、对于A1所在各个混合池中除A1外的所有克隆的k-mer交的集合求并集，得到各个混合池中排除克隆A1的并集：

B1'=A2∪A5∪A6 B2'=A2∪A3∪A4 B5'=A3∪A5∪A7

4、再将上一步得到的包含克隆A1的排除并集求交集：

A1-=B1'∩B2'∩B5'

5、用第2步得到的A1的k-mer交的集合减去上一步的交集，得到克隆A1的最终的k-mer集合：

A1~=A1-A1^-

对于上面的过程，可以用图形来表示。假设所有的克隆A1-A8(图9)，它们所拥有的k-mer用一个圆圈来表示，它们相互重叠的部分代表了共同拥有的k-mer。因为所有的克隆来自于同一个BAC文库，所以它们都包含有构建文库的载体的k-mer，即中心圆圈的部分。圆与圆之间的重叠部分越多，说明这两个克隆所拥有的相同的k-mer越多。

下面利用图形的相互叠加或排除计算来解释如何计算出克隆A1的最终k-mer集合。如图10A，是由8个克隆构建的6个混合池，对应到各个混合池所含有的克隆。然后对每三个混合池的k-mer集合求交集，如A1是由P1，P2和P5三个混合池交集结果，如图10B列出了所有的8个克隆计算结果的图形。其中颜色最深的部分，即为三个混合池的交集部分。对于克隆A1，存在于混合池P1，P2和P5，分别对P1，P2和P5中的不包括A1的其它所有克隆的k-mer集合求并集，得排除A1的并集P1’、P2’和P5’，如图10C所示。再对P1’、P2’和P5’求交集A1’，如图10D。最后用原A1的k-mer集合减去A1’，得到A1的最终k-mer集合A1^～，如图10E，图中深色部分即为A1的最终k-mer集合。

通过图形计算，可以直观地看出A1的最终k-mer集合并不是完整的A1真实的k-mer集合(Real k-mer set, RKS)，其中一部分的k-mer被从真实的k-mer集合中排除了。但是保留下来的k-mer都是真实k-mer集合中的元素。

经过以上计算，A1^~即为克隆A1的最终k-mer集合，其它克隆A2-A8的最终的k-mer集合的算法与之相同。

在上面的例子中，当对混合池的k-mer集合进行求交集以获得克隆A1的k-mer集合的A1，如果克隆A6，A4和A7的真实的k-mer集合包含共有的k-mer，那么混合池P1，P2和P5也会包含有这些共有的k-mer，在求A1时，这些k-mer也会被分配到A1中。但是实际上，这三个克隆共有的k-mer并不属于克隆A1，这就在A1中引入了错误的k-mer。这些错误的k-mer需要去除，所以才有了第3-5步。如果在A1的真实的k-mer集合与A6，A4和A7的真实的k-mer集合有共同的k-mer，则这些共同的k-mer也会出现在B1’、B2’和B5’中，即这三个集合相交得到的A1^-中也会含有这些k-mer。经过第5步的计算，这些k-mer会从A1中删除，也就是在A1中删除了这些真实的k-mer。所以第3-5步，在剔除假的k-mer的同时，也删除了一些真实的k-mer。尽管通过以上算法得到的克隆的最终k-mer集合与克隆真实的k-mer集合相比，少了一些k-mer，但是在得到的最终的k-mer集合中的元素都是属于真实的k-mer集合，也就是说其中不含有错误的元素。

以上推导的前提是在最理想的情况下，即每个混合池中所有克隆的序列都能被检测到，每个克隆的k-mer也都能被检测到，并且没有测序错误。但是实际上，由于实验条件，测序错误和测序的随机性等因素的影响，使得这种理想情况是不可能实现的。所以，A1中不可能包含所有的属于C1的真实的k-mer。同样，其它的克隆的k-mer集合中也不可能包含所有的属于其真实的k-mer。从而导致B1’中就会丢失一些属于混合池P1的k-mer，B2’和B5’中同理会丢失一些分别属于P2和P5的k-mer。在将B1’,B2’和B5’求交集得到A1^-时，A1^-中就会丢失一些k-mer。当这些丢失的k-mer属于A1，但又不属于P1的真实的k-mer时，这些错误的k-mer就会保留在P1的最终k-mer集合A1^~中，也就是说在A1^~中引入了假阳性。对于因测序错误而引入的错误k-mer集合，在解析克隆的k-mer集合之前，需要利用k-mer在混合池k-mer集合中出现的频数进行过滤。即如果k-mer在集合中出现的次数少于某一个数值，则认为该k-mer是因为测序错误引入的，从而去除混合池k-mer集合中大部分的错误的k-mer。过滤频数与测序深度和测序错误二者均有关，测序深度越大，测序错误越多，过滤频数则越大。如测序错误<3%,测序深度为20倍时，过滤频数选择1，即可以过滤绝大多数的错误k-mer。

考虑到克隆的最终k-mer集合的假阳性的引入，就需要增加测序深度以检测到每个混合池包含的尽可能多的真实的k-mer。要尽可能多地排除最终的k-mer集合中的假阳性，计算特定克隆的排除并集时，对并集进行平衡处理，这一步称为平衡混合池的k-mer集合。

Ⅲ.平衡混合池的特征序列集合和k-mer集合

假设在一个混合池中有个n克隆。将这个混合池中所有克隆k-mer的交的集合中的元素放入到一个KV-set中，所有的元素的频数设为1。当计算克隆的k-mer的交的集合的并集时(以上例子中的第3步)，对混合池中的所有克隆的KV-set求和，得到这个混合池的KV-set。然后比较混合池的KV-set中的键与混合池原来的k-mer集合中的元素，将原来k-mer集合中包含但是KV-set中不包含的元素添加到KV-set中，并将新加入键的频数设置为2，对于这样得到的KV-set称为这个混合池的全和KV-set(Sum KV-set)。

如果一个克隆在这个混合池中，则根据这个克隆的k-mer的交的集合构建一个KV-set C，C中元素为克隆的k-mer的交的集合中元素，所有的元素的频数都设为1。假设这个混合池的全和KV-set为P，则将P减去C，即P’=P-C，再对P’取键的集合Key(P’)，作为这个混合池排除给定克隆的并集。对于以上例子中，即为B1,B2和B5的k-mer集合分别排除A1后的B1',B2'和B5'。

关于序列定位与物理图谱的整合，更详细的算法如下：

1)克隆重叠群的分割

根据克隆重叠群中克隆的相对位置，将克隆的k-mer集合分割成更小的区块。假设有两个克隆属于同一个重叠群并相互重叠，集合A为其中一个克隆的k-mer集合，集合B为另一个克隆的k-mer集合。那么这两个克隆的k-mer集合可以分解成3个子集，即共有的k-mer集合(A∩B), A特有的k-mer集合(A-B)和B特有的k-mer集合(B-A)。两者共有的k-mer集合位于它们特有的k-mer集合之间，这样三个k-mer的相对位置就确定了。根据这一算法，即可对每一个克隆重叠群中所有克隆的k-mer集合，按照它们相对位置，从一端开始逐一分割成小的子集，并且这些子集的相对位置是确定的(图7)。

2)序列定位与整合

对各个混合池的NGS序列利用短序列拼装软件(如velvet)分别进行拼装，再将拼装结果合并，利用长序列拼装软件 (如PCAP)进行进一步拼装，得到序列重叠群。然后将每一个序列重叠群分解成k-mer集合，与上一步得到的所有小区块的k-mer集合求交集，找出其中交集元素总数最多的区块，这样就将序列重叠群定位到了克隆重叠区的小区块上(图7所示9条序列seq_1到seq_9分配到对应的小区块)。如果定位到同一个克隆重叠群的序列包含有较多的相互重叠的部分，则需要分别再次拼装合并相互重叠的序列，然后再次定位合并后的序列。

利用双末端比对软件(如Bowtie2)将所有混合池的NGS序列与定位后的序列进行双末端比对(图7所示，4对双末端PE1到PE4，高匹配到相对的序列)。如果两个序列与某一条序列的双末端高匹配(如Identity值≥97%)，并且这两个序列定位在同一个克隆重叠群，它们所定位的重叠群小区块在一定范围内，则认为这两条序列是可以连接在一起的。根据测序文库片段长度和双末端比对结果，计算它们之间的缺口长度，并用“N”填补其中的缺口。同时根据两条序列所在区块的相对位置，确定连接后序列的方向(图7中seq_4与seq_5可由PE2的两个末端连接，用seq_4-seq_5表示，seq_4与seq_5分别定位在小区块bin_4与bin_5中，所以连接后的序列seq_4-seq_5的方向可以确定)。如果一条序列与其它序列没有双末端序列连接，则将它分解的k-mer集合与相邻区块的k-mer集合比较，如果它与相邻区块的k-mer集合有一定的交集，则它可以定位到两个区块的交界处，从而确定它的方向(图7中,seq_3定位在小区块bin_3中，同时也与bin_4有较大的交集，所以可以根据它与两个小区块的k-mer集合的交集，确定seq_3的方向)。最后根据这些序列所在区块的相对位置，将所有定位的双末端连接后的序列连接成整个克隆重叠群的序列。如果双末端连接后的多个序列定位在相同的区块内，则它们之间的相对位置与方向随机选择(图7中，由PE3，PE4连接的3条序列seq_6-seq_7-seq_9，与seq_8都定位在bin_6中，所以它们之间的顺序与方向都无法确定，连接时随机选择它们方向和位置)。

本发明的有益效果在于：

（1）本发明无需对每个克隆进行酶切电泳，直接根据特征序列构建物理图谱。可用的特征序列数目可根据需要进行选择，并且数目庞大，不存在电泳条带长度读取的误差

（2）本发明利用NGS高通量测序技术与克隆混合池，不仅能够构建出全基因组的物理图谱，同时能够将拼装后的序列定位到物理图谱上，实现两者的整合，并且所有的数据都来自于同一个测序数据集。

（3）构建物理图谱与序列拼装是可以同时进行的，对于大的基因组可以分割成很多小工程同时进行，以提高效率。

（4）得到的基因组全序列，由于物理图谱的存在，在大的框架上的错误很少。对于小区域里的错误，在序列组装过程中可以给出明确的位置。有利于以后对基因组全序列的缺口填补与错误纠正。

（5）本发明具有较为普遍的适用性，无论是全基因组鸟枪法测序(WGS)还是第二代测序技术(NGS)，或是将来的第三代测序技术，都可以适用该方法。同时，可以整合这些测序方法，充分利用它们的优点，以提高序列的精确性和完整性。

（6）具有很高的灵活性，由于特征序列是绝对数值，如果以后想进一步完善基因组序列，可以很容易添加并整合。

附图说明

图1为本发明的方法流程图；

图2为混合池构建的直角坐标系统；立方体中的每一个克隆都可以由三个坐标值(x, y, z)确定；

图3 为定序混合池的构建策略图；图中每个立方体中，相同的颜色或材质标记的克隆属于同一个混合池；

图4为基础混合池的构建图；同一行或同一列中连续的8个克隆在同一个基础混合池中；

图5为特征序列的定义图；对于一条DNA序列，在指定标记前缀序列(如“AATT”)和特征序列长度(如31)时，前缀序列上游的一定长度的碱基和前缀序列反向互补序列上游的同样长度碱基为这条DNA序列的一条(对)特征序列。前缀序列可以指定多个，特征序列的长度也可以根据需要指定。

图6为重叠群157。该克隆重叠群共有32个克隆，图中突出显示的克隆为克隆1225。

图7为分割克隆重叠群成多个小区块的示意图。每一个小区块(Bin)包含相应的k-mer集合。bin表示克隆重叠群分割后的小区块，共15个；seq表示定位到小区块的序列，共9条；PE表示NGS的双末端，共4个，PE1_1表示PE1的一个末端，PE1_2表示PE1另一个末端，每一个PE两端的符号“+”或“-”表示对应末端序列与定位序列的比对方向，如果两个末端的符号相同，则在连接两条序列时，把其中一条序列进行反向互补。

图8为在同一个立方体中的8个克隆的排列图；

图9为克隆的相互重叠关系图。每一个圆圈的面积表示对应的克隆所拥有的k-mer总数，各个圆圈相互重叠部分的面积，反应了共同拥有的k-mer数目，其中心的圆圈部分表示所有的克隆都有的k-mer，其面积反应构建BAC文库时的载体序列的k-mer数目，而中心圆圈外其它所有圆圈内的部分，表示所有的克隆都有的相同的k-mer，这一部分的k-mer很可能来自于重复序列。

图10为图形计算过程，B中颜色最深的部分为计算的结果。

具体实施方式

为了更好地解释本发明，以下结合具体实施例进一步阐明本发明的主要内容，但本发明的内容不仅仅局限于以下实施例。

实施例1

对拟南芥(Arabidopsis thaliana ecotype Columbia)进行全基因组测序，该物种共5条染色体，全基因组序列约122Mb。注意，该实施例中再现的数据仅用于本发明实施过程的说明，不用作其它途径。

构建细菌人工染色体（BAC）文库

利用限制性内切酶BamHI对全基因组DNA进行部分酶切，构建约5倍全基因组覆盖度的BAC文库，文库保存在11块384孔块中，共4224个BAC克隆；BAC克隆的平均插入片段大小为137kb，片段长度在60~300kb之间。

对BAC克隆进行编号，从0到4223。

构建克隆混合池

从BAC文库中挑选N个克隆，将这些克隆排列成一个n×n×n立方体。然后将这个立方体(图2)置于右手坐标系统(right-hand rectangular coordinate system)，X,Y,Z为坐标轴，O为原点。这样每一个克隆的位置都可以由三个坐标值(x、y和z)来确定，克隆在立方体中的位置用(x, y, z)表示。在这样一个由BAC克隆构成的立方体基础上构建克隆的混合池。同一个混合池中的克隆集合用{(x1, y1, z1)，(x2, y2, z2)，(x3, y3, z3)，…，(xn, yn, zn)}表示。

当所有克隆在一个立方体中的位置确定后，就不再变动，需要构建多维的混合池时，只是通过不同的面上的克隆重新组合成新的混合池(图3)。考虑到实际操作的简单性，混合池的立方体中，每个面上同一列或同一行的克隆永远属于同一个混合池。根据这一原则，对于一个确定的立方体，有9种混合池构建的策略，每一种策略构建的混合池称为1维。将这9种策略分成3类，分别为直交的混合池(Perpendicular crossing pools，包括PX, PY和PZ)，斜交的混合池(Oblique crossing pools, 包括OX, OY and OZ)和角交的混合池(Angle crossing pools,包括AX, AY and AZ)。如图3，用一个面与这个立方体相交，用相同的颜色或材质标记的克隆属于同一个混合池。

在PX中，垂直于X轴的面与这个立方体相交，同一个面上的克隆属于同一个混合池，即所有的克隆位置中x值一样的克隆在同一个混合池中，如x值为0时，这些克隆包括：

{(0, 0, 0)，(0, 1, 0)，(0, 2, 0)…(0, n-1, 0)，

(0, 0, 1)，(0, 1, 1)，(0, 2, 1)…(0, n-1, 1)，

(0, 0, 2)，(0, 1, 2)，(0, 2, 2)…(0, n-1, 2)，

(0, 0, 3)，(0, 1, 3)，(0, 2, 3)…(0, n-1, 3)，

…

(0, 0, n-1)，(0, 1, n-1)，(0, 2, n-1)…(0, n-1, n-1)}

共n×n个克隆。同理x值为任意值时，属于同一个混合池中的克隆为：

{(x, 0, 0)，(x, 1, 0)，(x, 2, 0)…(x, n-1, 0),

(x, 0, 1)，(x, 1, 1)，(x, 2, 1)…(x, n-1, 1),

(x, 0, 2)，(x, 1, 2)，(x, 2, 2)…(x, n-1, 2),

(x, 0, 3)，(x, 1, 3)，(x, 2, 3)…(x, n-1, 3),

…

(x, 0, n-1)，(x, 1, n-1)，(x, 2, n-1)…(x, n-1, n-1)}

这样，可以得到n个混合池。与PX相似，PY与PZ也可以分别构建出n个混合池。

在OX中，与立方体相交平面垂直于平面YOZ，且这个平面在YOZ面上投影直线的斜率为-1。当平面在YOZ面上投影为对角线(图3中OX中白色方块连线)，则图中所有的白色标记的克隆为同一混合池。当平面在YOX面上的投影不是对角线时，则将两个面与立方体相交位置的克隆合并为一个混合池，如图中黑白网格标记的克隆，对应的两个面在YOX面上的投影位置为直线y=-z+2和直线y=-z+(2n+2)。其它的同理，这样得到n个混合池，每个混合池中包含n×n个克隆。

OY，OZ与OX相似，也可以得到n个混合池，每个混合池包含n×n个克隆。

AX可以参考OX，与OX不同的是，AX中，与立方体相交的平面，在平面YOZ上的投影直线的斜率为1。AY与AZ可以分别参考OY与OZ。

每种策略都得到n个混合池，每混合池包含n×n个克隆，使每个混合池中所拥有的基因组片段总长度尽可能地接近。混合池的维数(Pool dimension)是指构建混合池时每一个克隆在所有的混合池中出现的次数总合。例如，3维的混合池(3D pool)可以包括PX,PY和PZ或其它任意的三个立方体构建的混合池；6维的混合池(6D pool)和9维的混合池(9D pool)与3 维的混合池相似。为了定位一个克隆，至少需要3个不同维的混合池，但是并不是任意的3个混合池都能够定位一个克隆。如表1中所列出的所有的任意3个混合池的84种组合结果中有69个是可以确定克隆位置的。3个混合池可以定位一个克隆，则可以用6个混合池对同一个克隆定位两次或用9个混合池对同一个克隆定位3次。如果将克隆在立方体中的位置进行重新排列，则可以得到一个新的立方体，这个新的立方体又可以构建出9种混合池，从而得到更高维的混合池以提高精度和质量。

表1.9种克隆混合池中3个混合池的组合结果。

表中“TRUE” 表示这三个混合池的交集只有一个克隆，“FALSE” 表示这三个混合池的交集多于1个克隆。在进行混合池构建时，最好选择能够定位出一个克隆的三个混合池组合。

从BAC文库中选择4096个克隆，放入一个16×16×16的立方体中，并选择PX，PY，PZ，OX，OY与OZ这6种策略构建6维的混合池(若有更高质量的要求，可以构建9维的混合池，或更高维)。每维有16个混合池，共96个混合池，每个混合池包含256个克隆。

96个混合池与其所包含的克隆编号如下表：

表2.混合池与其所包含的克隆编号对应表

混合池DNA提取

4096个克隆从文库中挑选出来，接种到96深孔培养板中，每个培养孔中含有1ml的冷冻培养基(1L冷冻培养基配方：蛋白冻10g，酵母提取物5g，NaCl 10g，K₂HPO₄.3H₂O 8.215g，KH₂PO₄ 1.795g，柠檬酸钠0.5g，(NH₄)₂SO₄ 0.848g，甘油44ml，用前加入1M/L MgSO₄ 0.4ml)，在37°C下培养18个小时，于-80°C保存。将每一块96孔培养板的每一行每一列中连续的8个克隆相混合，组成基础的混合池(Base pool)。如果同一行中连续的克隆数不到8个，则将下一块板的同一行中的克隆并入这个基础混合池，共得到1024个基础混合池(图4)。每一个培养孔中吸出150μl的菌液，这样每一个基础混合池中含有1.2ml的菌液。然后根据定序混合池的构建策略，将4096个克隆放入一个16×16×16的立方体中，构建6D混合池，每个混合池中含有16×16=256个克隆。对于每一个混合池，其含有的克隆与基础混合池相比较，则32个基础混合池可再次组合成含有256个克隆的测序混合池。将基础混合池复制到新的96孔培养板中，每一个孔中吸入100μl菌液和1.1ml的2×YT培养基。并且属于同一个测序混合池中的基础混合池在96孔培养板中是相邻的，其中第1到4列共32个基础混合池组成一个测序混合池，第9到12列组成另一个测序混合池，中间的4列空出防止交叉污染。然后将新复制的96孔培养板在37°C下培养18小时备用。最后将每一个新复制的96孔培养板上连续排列的、属于同一个测序混合池的基础混合池复制15份到5块新的培养板中，每一个新的孔中含有40μl菌液和1ml的2×YT培养基；在37°C培养18小时。将新培养的每一个测序混合池中的克隆混合在一起利用QIAGEN公司的Large-Construct Kit试剂盒提取质粒DNA。质粒DNA浓度要求≥150ng/μl，总量要求≥30ug。此过程也可以先分别提取所有克隆的DNA后，再进行混合。

混合池DNA测序

利用Illumina公司的HiSeq 2000高通量测序仪对96个混合池的质粒DNA进行双末端测序。其中PX，PY，PZ混合池的DNA构建测序文库时，片段长度为500bp；OX混合池DNA构建测序文库时，片段长度为800bp；OY混合池DNA构建测序文库时，片段长度为1200bp；OZ混合池DNA构建测序文库时，片段长度为1800bp。序列读长两端均为100bp；每个混合池得到约710Mb的序列总长，约20倍的覆盖度；序列平均错误率在3%以下，经测序质量值计算平均错误率为2.64%。

扫描测序结果并得到混合池的特征序列集合

特征序列是用来进行物理图谱的拼装，作为克隆所拥有的DNA片段的特征。它相当于传统的利用指纹图谱构建物理图谱方法中的条带(Band)长度。如图5，对于一条DNA序列，在指定标记前缀序列为“AATT”和特征序列长度为31时，前缀序列上游的31个碱基和前缀序列反向互补序列上游的31碱基为这条DNA序列的特征序列。前缀序列可以指定多个，特征序列的长度也可以根据需要指定。

在具体实施中，选择前缀序列为两个限制性酶切位点BamHI(G|GATCC)和EcoRI(G|AATTC)，特征序列长度为31bp。对每一个混合池的所有序列进行逐一扫描，将得到的所有的特征合并为一个集合，这个集合称为这个混合池的特征序列集合(Feature sequence set, FS-set)，特征序列集合中不包含有重复的元素。在扫描混合池的特征序列时，也记录每一个特征序列出现的次数。96个混合池得到96个特征序列集合，每个混合池的特征序列集合大约包含98000个特征序列。

解析克隆的特征序列集合

对得到的各个混合池的特征序列集合，删除其中只出现一次的特征序列，并对处理后的混合池的特征序列集合进行编号，其编号与混合池编号相对应：

PX0,PX1,PX2…PX15；

PY0,PY1,PY2…PY15；

PZ0,PZ1,PZ2…PZ15；

OX0,OX1,OX2…OX15；

OY0,OY1,OY2…OY15；

OY0,OY1,OY2…OY15；

分别表示PX，PY，PZ，OX，OY与OZ混合池策略中的各16个混合池的特征序列集合。

按以下步骤解析克隆的特征序列集合：

1）求克隆在所有混合池的特征序列集合的交集。如编号为1225的克隆在混合池PX0、PY15、PZ0、OX0、OY15和OZ0这6个混合池中，对这6个混合池的特征序列集合求交集PX0∩PY15∩PZ0∩OX0∩OY15∩OZ0 (运算符“∩”表求两个集合的交集)得到克隆1225的特征序列的初始交集I₁₂₂₅。其它克隆与之类似，求出所有克隆的特征序列的初始交集，用I_i表示(其中下标i表示克隆编号)，如I₁₁₉₁表示编号为1191的克隆特征序列初始交集。

2）求克隆所在混合池的排除并集。克隆所有混合池的排除并集用符号U_i(p)，其中i为克隆编号，p为克隆所在的混合池编号。以克隆1225为例，它所在的混合池PX0的排除并集为其它所有克隆(不包含克隆1225)的初始交集的并集，克隆1225在混合池PX0的排除并集

U_1225(PX0)=I₆₈ U I₅₀₉ U I₃₉₀₄ U … U I₃₇₆₁ U I₃₆₆₃ U I₃₀₉₁;(运算符“U”表示求两个集合的并集；初始并集的顺序与表1中PX0中包含的克隆编号出现的顺序一致，但是不包含克隆1225)

同理，克隆1225在混合池PY15、PZ0、OX0、OY15和OZ0的排除并集分别为：

U_1225(PY15)=I₁₁₉₁ U I₃₂₉₄ U I₂₀₇₃ U … U I₃₄₀₃ U I₃₂₅₃ U I₁₀₈₅ ;

U_1225(PZ0)=I₂₁₅₁ U I₄₄ U I₃₆₆₁ U … U I₄₃₆ U I₃₃₂₉ U I₁₉₈₁ ;

U_1225(OX0) = I₆₈ U I₅₀₉ U I₃₉₀₄ U … U I₈₂ U I₃₄₈₇ U I₁₄₂₈ ;

U_1225(OY15) = I₆₈ U I₅₀₉ U I₃₉₀₄ U … U I₁₅₄₇ U I₁₇₉₁ U I₉₇₃ ;

U_1225(OZ0) =I₁₁₉₁ U I₃₂₉₄ U I₂₀₇₃ U … U I₃₉₂₃ U I₃₈₀₃ U I₁₉₈₁

3）求克隆的最终特征序列集合。用符号Fi表示，其中下标i为克隆的编号。对每个克隆，用它的初始特征序列交集减去它所在混合池的排除并集的交集得到它的最终特征序列集合。如克隆1225，其最终特征序列集合为

F₁₂₂₅=I₁₂₂₅−(U_1225(PX0)∩U_1225(PY15)∩U_1225(PZ0)∩U_1225(OX0)∩U_1225(OY15)∩U_1225(OZ0))

其中运算符“−”表示求两集合的差集。

关于克隆特征序列集合解析更具体的算法参见上述克隆特征序列集合与k-mer集合解析算法。

构建克隆重叠群

将得到的所有克隆的特征序列求并集，对并集中的特征序列建立索引表，每一个特征序列对应一个索引值，索引值从1开始。把每一个克隆的最终特征序列集合中的特征序列，用对应的索引值替换，然后对索引值从小到大排序，再导出排序后的索引值到克隆拼装软件FingerPrinted Contig (FPC,V9.4)兼容的“.size”文件。导出到“.size”文件中的索引值就相当于指纹图谱构建物理图谱方法中的条带(Band)长度。所有克隆最终特征序列的索引值导出的“.size”文件内容如下：

3280133Gel

1

2

3

4

5

6

127

128

129

130

131

132

133

-1

……

1225117Gel

759

761

1839

1841

1857

…

21287

38090

38092

89023

93016

-1

……

102280Gel

273

281

301

343

17787

17788

17789

17790

17791

53407

53417

53426

66412

66414

97491

97492

-1

每一个克隆的“.size”文件内容包含三部分，第一部分包含克隆名称(编号)，索引值数目(条带数目)和“Gel”，第二部分为各索引值，第三部分为结束标记“-1”。 如克隆1225最终得到117个特征序列，对应117个索引值。

将导出的“.size”文件出入到拼装软件FPC中，进行拼装。将Tolerance设为0(这个值一定要设为0)，Cutoff设为1e-12，其它值默认，进行初次拼装产生克隆重叠群。然后，设置if >=5 Qs，Step为5，进行DQer。最后得到220个克隆重叠群，所有重叠群中包含共4021个克隆，还有75个克隆不在任意一个重叠群中。如图6为重叠群157(ctg157)中克隆的排列情况。

扫描测序结果并得到混合池的k-mer集合

k-mer的定义与基于Bruijn graph进行短序列拼装方法中的k-mer定义是一样的。对于某一条DNA序列，取一定长度的窗口，由5’向3’一个碱基一个碱基地滑动，每滑动一次，窗口内的子序列就是这条DNA序列的一个k-mer。一条序列在窗口长度为13时，列出的所有k-mer的情况，共34个k-mer。所以，对于一条长度为L的序列，在窗口长度(或k-mer长度)为k时，可以计算出产生的k-mer的最大数目为L-k+1。

实施例中，k-mer长度为31bp。对每一个混合池的所有序列进行逐一扫描，将得到的所有k-mer合并为一个集合，这个集合称为这个混合池的k-mer集合(k-mer set,K-set)，k-mer集合中不包含有重复的元素。在扫描混合池的k-mer时，也记录每一个k-mer出现的次数。96个混合池得到96个k-mer集合，每个混合池的k-mer集合大约包含95000000个特征序列。

解析克隆的k-mer集合

由于特征序列与k-mer的本质是一样的，所以解析克隆特征序列的算法同样适用于解析克隆的k-mer集合。在具体计算过程中，只是把特征序列集合用k-mer集合替换，其它过程完全一致。

分割重叠群

分割克隆重叠群的算法详细见上述“关于序列定位与物理图谱的整合”。

实施例中，如重叠群49(ctg49)被分割成27个小区块,共635264个k-mer，从一端到另一端各个小区块中k-mer的数目分别为85353、24951、3044、3992、1284、1175、172、125881、12172、8580、4305、16549、1959、26952、12668、6886、6195、6364、20721、6553、4309、4893、26776、18726、12071、83478、1680和107575。其它重叠群与之类似，都被分割成小区块。

序列组装并整合到物理图谱上

对各个混合池的NGS序列利用短序列拼装软件velvet分别进行拼装，得到各个混合池的序列重叠群(Sequence contig)。拼装时，Hash length参数值取31，其它参数均为默认值，得到96个混合池的初次拼装序列，共约3.4Gb，其中最小的序列长为62bp。然后，对初次拼装后的序列进行过滤，去除其中小于100bp的序列，将剩下的序列合并到一起，共得到1.59Gb的“fasta”格式的序列，并压缩成“gzip”格式。根据合并的序列，产生对应的“fasta”格式的质量文件，所有碱基的质量值都取20，并将产生的质量文件也压缩成“gzip”格式。

将合并后的序列利用长序列拼装软件PCAP再次进行拼装，拼装参数均为默认参数，拼装后得到约270Mb的序列。

将PCAP拼装后的每一条序列，分解k-mer集合，与所有克隆重叠群分割后的小块的k-mer集合求交集，找出其中交集元素最多一个，把这个序列定位到这个小块的位置。经过对定位到同一个克隆重叠群的序列进行分析，很容易找出其中含有较多的相互重叠的序列，所以需要对这些序列再次拼装合并，再次定位。将定位到同一个克隆重叠群上的序列用序列拼装软件Phrap分别再次进行拼装，参数为默认参数。再将拼装后的序列再分解成k-mer集合，再次定位到重叠群的小区块上。经过这一步，定位的序列共有103.1Mb。

利用双末端比对软件Bowtie2将所有混合池的NGS序列与定位后的序列进行双末端比对。如果两个序列与某一条序列的双末端高匹配(实施例中选择Identity值≥97%)，并且这两个序列定位在同一个克隆重叠群，它们所定位的重叠群小区块在一定范围内，则认为这两条序列是可以连接在一起的，根据测序文库片段长度和双末端比对结果，计算它们之间的缺口长度，并用“N”填补其中的缺口。同时根据两条序列所在区块的相对位置，确定序列方向。如果一条序列与其它序列没有双末端序列连接，则将它分解的k-mer集合与相邻区块的k-mer集合比较，如果它与相邻区块的k-mer集合有交集，即它可以定位到两个区块的交界处，则可以确定它的方向。对于定位到区块内而无法确定方向的序列，则随机选择它的方向。最后根据这些序列所在区块的相对位置，将所有定位的并经双末端连接后的序列连接成整个克隆重叠群的序列。如果双末端连接后的多个序列定位在相同的区块内，则它们之间的相对位置与方向随机选择。

其中对于测序文库片段长为500bp的混合池的NGS序列，与定位后的序列进行双末端比对时，匹配率(Identity)要求≥97%，两个序列在重叠群小块区的相对位置在10以内(相对位置范围称为区块范围，如从重叠群的左端开始，最左端的第一个区块的相对位置为1，第二个为2，依次类推)。如果不一样，刚根据两条序列与双末端的匹配位置，以及默认长度500，计算出它们之间的缺口长度。如果它们之间的缺口长度小于100bp(这个长度称为最小缺口长度)，再将两条序列邻近末端的13个碱基进行比较，如果这13个碱基是一样的，则直接将这两条序列连接成一条序列，否则用对应数目的“N”填补上缺口。其它测序文库的NGS序列处理与之类似，只是区块范围(500bp测序文库为10)和最小缺口长度不一样。800bp，1200bp，1800bp测序文库的区块范围分别为17，22，31，最小缺口长度分别为120bp，150bp，180bp。经过序列双末端的连接并填补序列之间的缺口得到220个重叠群的序列，共得到109.9Mb。

重叠群定位并连接成染色体

如果已有分子标记，将每一个分子标记的序列分解成k-mer集合，与每一个混合池的k-mer集合(已删除其中只出现一次的k-mer)求交集，如果混合池的k-mer集合包含标记序列k-mer集合，则将这个序列分配到这个混合池。如此得到混合池的标记集合，运用解析克隆特征序列集合的算法，解析出克隆的标记集合。如果通过Southern印迹杂交或PCR，则可以直接确定BAC文库中的克隆所包含的标记(Marker)。这样就确定了包含标记的克隆所在的染色体和它们的相对位置。再根据克隆所有的重叠群与克隆包含的标记所在的染色体，将克隆重叠群分配到染色体上。由重叠群中克隆的相对位置与克隆包含的标记在染色体上的相对位置，确定重叠群的方向。最后定位在同一染色体上的克隆重叠群的序列连接起来成为整条染色体。

实施中，共用了2072个已知分子标记(来自http://www.arabidopsis.org/)，过滤掉其中序列小于100bp的标记，剩下1515个标记。如果一个混合池的k-mer集合(已删除其中只出现一次的k-mer)包含某一标记序列k-mer集合的90%元素，则将该标记分配到该混合池，再解析出克隆的标记集合，利用PCR验证各个克隆是否包含解析出的分子标记，其中2953个克隆含有一个或多个分子标记。通过克隆的标记集合，将211个克隆重叠群定位到了染色体上，再将同一染色体上的克隆重叠群序列连接起来，相邻重叠群之间填补50kb长度的“N”。最终得到5条染色体序列，共118.4Mb。

上述实施例1是以拟南芥为实验模型，在实际过程中，根据不同物种全基因组构建BAC文库，建立不小于3维的混合池，通过以上的计算方法得到全基因组序列。

Claims (1)

1.基于克隆DNA混合池的全基因组测序方法，其特征在于：包括以下步骤：

1)提取全基因组DNA，构建BAC文库；

2)构建BAC克隆混合池；

3)提取BAC克隆混合池的DNA；

4)对步骤3)中的BAC克隆混合池的DNA利用NGS进行双末端测序；

5)扫描各个混合池的序列，获得各个混合池的特征序列集合与k-mer集合；

6)根据混合池的特征序列集合与k-mer集合，解析出各个克隆的特征序列集合与k-mer集合；

7)利用克隆的特征序列集合构建克隆重叠群；

8)利用克隆重叠群将克隆的k-mer集合分割成小的k-mer集合并定位到克隆重叠群上；其中，

克隆重叠群的分割：

根据克隆重叠群中克隆的相对位置，将克隆的k-mer集合分割成更小的区块，假设有两个克隆属于同一个重叠群并相互重叠，其中集合A为其中一个克隆的k-mer集合，集合B为另一个克隆的k-mer集合，那么这两个克隆的k-mer集合分解成3个子集，即共有的k-mer集合(A∩B)，A特有的k-mer集合(A-B)和B特有的k-mer集合(B-A)；两者共有的k-mer集合位于它们特有的k-mer集合之间，这样三个k-mer的相对位置就确定了，根据这一算法，对每一个克隆重叠群的所有的克隆，按照它们相对位置，从一端开始逐一分割，从而将一个克隆重叠群中所有克隆的k-mer集合分割成小的子集，并且这些子集的相对位置是确定；

9)对步骤4)中混合池的NGS序列进行拼装得到序列重叠群，其中，序列定位与整合：

对各个混合池的NGS序列利用短序列拼装软件如velvet分别进行拼装，再将拼装结果合并，利用长序列拼装软件如PCAP进行进一步拼装，得到序列重叠群，然后将每一个序列重叠群分解成k-mer集合，与上一步得到的所有小区块的k-mer集合求交集，找出其中交集元素总数最多的区块，这样就将序列重叠群定位到了克隆重叠区的小区块上，如果定位到的同一个克隆重叠群的序列包含有相互重叠的部分，则需要分别再次拼装合并相互重叠的序列，然后再次定位合并后的序列；

10)将序列重叠群定位到克隆重叠群上，并利用测序的双末端信息连接序列重叠群，确定它们的方向，得到克隆重叠群的序列；

11)利用分子标记确定克隆重叠群的相对位置与方向，将克隆重叠群的序列连接成整条染色体序列，得到全基因组序列；

其中，关于克隆特征序列集合与k-mer集合的解析的总算法如下：

某一物种BAC文库的克隆总数为a，构建的混合池总数为x，则；

第κ维，索引为λ混合池的k-mer集合表示为：P_(κ,λ)；

包含某一给定克隆的混合池的k-mer集合为：{P_(δ)|δ<m,P_(δ)∈P}；

包含同一克隆的混合池的k-mer集合的交集，即克隆的IKS为：

某一克隆在混合池中的排除并集EUKS为：

某一克隆k-mer集合的所有排除并集的交集为：

某一克隆的最终k-mer集合FKS为：CF＝C_τ-CI_τ。