CN102831331A

CN102831331A - 基于酶切建库双末端测序的长度多态性标记的引物设计开发方法

Info

Publication number: CN102831331A
Application number: CN2012102302331A
Authority: CN
Inventors: 郑泽群; 任一; 陶晔; 胡秋萍; 黄华生
Original assignee: SHANGHAI MAJORBIO PHARM TECHNOLOGY Co Ltd
Current assignee: SHANGHAI MAJORBIO PHARM TECHNOLOGY Co Ltd
Priority date: 2012-07-04
Filing date: 2012-07-04
Publication date: 2012-12-19
Anticipated expiration: 2032-07-04
Also published as: CN102831331B

Abstract

本发明公开了基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，它采用了生物信息学分析方法，处理RAD双末端测序的测序数据，从而寻找RAD测序片段上的Indel位点信息，以突破非模式生物缺少参考序列的瓶颈，简化了基因组的复杂度，同时也减少了测序成本。

Description

基于酶切建库双末端测序的长度多态性标记的引物设计开发方法

技术领域

本发明涉及一种基因组长度多态性标记的引物设计方法。具体为基于酶切建库双末端（Pair-end）测序的长度多态性标记的引物设计开发方法；在缺少参考序列的情况下，寻找到个体间的Indel标记，并能够在两端设计引物。属于生物信息学技术领域。这对于缺少参考序列的非模式生物的研究具有重要的意义。

背景技术

InDel(insertion-deletion)***缺失标记，指的是两种亲本中在基因组上的差异，相对另一个亲本而言，其中一个亲本的基因组中有一定数量的核苷酸***或缺失。Indel位点信息的获得可以有许多重要的应用，如构建遗传图谱，基因分型，分子标记辅助育种，疾病检测等。

如今，第二代DNA测序技术是一种高通量低成本的测序技术，基本原理是边合成边测序。以solexa测序方法为例，先用物理方法将DNA链随机打断，然后在片段两端加上特定接头，接头上有扩增引物序列。测序时，DNA聚合酶合成待测片段的互补链，通过检测新合成碱基所携带的荧光信号读取碱基序列，从而获得待测片段的序列。

第二代测序技术已经广泛应用于生物科学的许多领域，特别是研究一个物种不同个体之间的多态性。传统Call Indel标记的方法是将测序个体得到的短reads通过比对软件比对回参考序列，从而得到测序个体的Indel信息。常见的流程有：使用BWA软件将reads比对回参考序列，使用SAMtools软件处理比对结果寻找Indel位点^1,2。大体过程如图1所示。

目前，有参考序列的物种都可以很方便的进行Indel标记的查找，并在两端设计引物进行实验验证。但是对于那些非模式生物而言，基本上是不存在参考序列的。而在没有参考序列的情况下，传统寻找Indel标记的方法存在着技术上的瓶颈。

1.Li H.and Durbin R.(2009)Fast and accurate short read alignment with Burrows-Wheeler Transform.Bioinformatics,25:1754-60.[PMID:19451168]

2.Li H.*,Handsaker B.*,Wysoker A.,Fennell T.,Ruan J.,Homer N.,Marth G.,Abecasis G.,Durbin R.and 1000 Genome Project Data Processing Subgroup(2009)The Sequence alignment/map(SAM)format and SAMtools.Bioinformatics,25,2078-9.[PMID:19505943]

RAD-PE测序技术采用了新的建库方式(酶切建库)，其测序具体过程如图2所示，用限制性内切酶切断DNA特定的位点，再用物理方法将酶切之后的DNA分子随机打断，通过琼脂糖胶DNA分离技术挑选特定长度的DNA分子，然后在挑选出来的DNA末端添加特定的扩增接头与测序接头，从而构建上机文库进行高通量测序。

其中RAD测序方法为本领域公知的方法，例如可参考以下文献：

(1)Michael R Miller,Tressa S Atwood,B Frank Eames,et al,RAD marker microarrays enable rapid mapping of zebrafishmutations,Genome Biology,2007,8(6):R105.1-R105.10;

(2)Michael R.Miller,Joseph P.Dunham,Angel Amores,et al,Rapid and cost-effective polymorphism identificationand genotyping using restriction site associated DNA(RAD)markers,Genome Research,2007，17，240-248；

(3)Nathan A.Baird1,Paul D.Etter,Tressa S.Atwood,et al,Rapid SNP Discovery and Genetic Mapping Using Sequenced RAD Markers,PLoS ONE,2008,3(10),e3376,doi:10.1371/journal.pone.0003376.

散列表（Hash table，或哈希表），是根据关键码值(Key value)而直接进行访问的数据结构。也就是说，它通过把关键码值映射到表中一个位置来访问记录，以加快查找的速度。这个映射函数叫做散列函数，存放记录的数组叫做散列表。使用哈希表对数据进行索引基本是随着数据量的上升线性增长，而且由“ATCGN”构成的字符串，键值出现冲突的可能性非常低。这样在处理海量测序数据的时候有着很好的性能。

发明内容

本发明的目的是提供一种基于酶切建库双末端（Pair-end）测序的长度多态性标记的引物设计开发方法；它是一种通过处理基于酶切建库pair-end测序（RAD-PE测序技术）得到的测序数据，在两个个体之间寻找长度多态性位点，并能够在两端侧翼序列设计引物进行实验验证的技术方案。

本发明的目的通过以下技术方案来实现：

基于酶切建库双末端（Pair-end）测序的长度多态性标记的引物设计开发方法，其步骤如下：

1）在获得RAD高通量测序技术的测序结果后，对RAD双末端测序序列进行过滤以去除不合格的测序序列。

其中，RAD高通量测序技术可以为Illumina GA测序技术，也可以为现有的其他高通量测序技术。

所述的不合格的测序序列为测序质量低于预定的低质量阈值的碱基个数超过整条序列碱基个数的50%的序列。

2）根据测序个体基因组酶切一端的测序序列，利用序列的全同性生成每个个体堆的信息。例如，将每个个体过滤后的酶切一端的测序序列信息作为哈希的键，哈希的值指向一个链表，用于存放另一端的序列信息，并计算测序深度信息。可用任何一种编程语言实现该过程。

3）过滤掉酶切一端序列测序深度为1的结果（成对过滤）。

4）两个个体内分别将酶切一端的测序序列数据进行不容许空隙的两两比对，对堆进行聚类以确定个体内在酶切一端序列上的杂合SNP信息。

所述的不容许空隙的两两比对是指比对的时候不容许开空位。

其中，只有一个堆的聚类结果表明在酶切一端测序片段上不存在杂合位点，只有两个堆的聚类结果表明在酶切一端测序片段上存在杂合位点，一般情况下这个杂合位点不会处于重复区域。对于那些堆的个数超过两个的聚类结果，通常由于酶切一端测序序列处于基因组的重复区域所造成的，因此这些聚类结果将会被过滤掉。

所使用的比对软件可以是任何一款序列比对软件，如blast、blat等。计算堆的个数为一和二的聚类结果的总深度，并进一步过滤掉低深度和高深度的聚类结果。低深度的阈值通常为平均测序深度的四分之一，高深度的阈值通常为平均测序深度的两倍。

5）在两个个体内部，对每个堆的另一端数据进行局部组装，采用的组装可以是任何一款组装软件，如基于重叠关系的组装软件phrap，如基于De Brui jn graph算法的组装软件SOAPDenovo。

利用重叠关系进行组装的时候通常不容许空位的存在，其他参数可采用软件默认设置。利用De Brui jn graph算法进行组装的时候，通常kmer的大小要在21到31之间，其他参数可采用软件默认设置。

6）利用两个个体酶切一端的测序序列信息将两个个体堆的信息相互进行两两对齐，即在个体A和个体B中，个体A的某个堆能够和个体B的某个堆对齐，当且仅当两个个体堆中的酶切一端的测序序列完全相同。对能够对齐的堆，两个个体之间的另一端的组装结果序列相互进行比对，来寻找Indel位点信息。比对的时候，两个个体之间容许的空位数为一、并且公共区域内不存在错配。使用的比对软件可以是任何一款序列比对软件，如blast、blat等。

通过以上的步骤，将会得到两个个体之间，高可信度的Indel位点信息，还得到了在Indel位点周围的侧翼序列信息。

7）最后可在Indel位点周围的侧翼序列上设计引物，以便于后续的大规模实验应用。如构建遗传图谱，个体基因分型等。

本发明的技术方案采用了生物信息学分析方法，处理RAD（restriction-site associated DNA，which are short fragments of DNA adjacent to each instance of a particular restriction enzyme recognition site）双末端测序的测序数据，从而寻找RAD测序片段上的Indel位点信息，以突破非模式生物缺少参考序列的瓶颈，简化了基因组的复杂度，同时也减少了测序成本。常规的过程针对的物种通常为二倍体生物，但并不仅局限于二倍体生物。

附图说明

图1为传统Call Indel标记方法的过程原理图；

图2为RAD-PE测序技术的测序具体过程原理图；图中，(A)限制性内切酶酶切基因组DNA，并加上P1接头，每一个P1接头含有不同的序列标签；(B)带有不同P1接头的样品混合在一起，打断；(C)加上接头P2；(D)扩增富集RAD tags；

图3为RAD双末端测序的例图；

图4为生成堆的聚类过程图；

图5为堆中酶切一端序列信息示意图；

图6为堆中另一端序列信息示意图；

图7为堆中测序数据过滤和整合流程示意图；

图8为局部组装和Indel位点查找流程示意图；

图9为引物设计方式示意图；

图10为RAD序列测序深度分布图。

具体实施方式

下面结合附图与具体实施例进一步阐述本发明的技术特点。

如图2所示，获得两个个体基因组的RAD双末端（pair-end）测序序列。图3为RAD双末端测序的例图。在图3中显示了用限制性内切酶Ecor1，识别DNA分子上“G^AATTC”的回文序列，并在G与A之间将DNA分子切断，将酶切后的DNA分子用物理方法打断成短的序列片段，并在含有酶切末端序列的DNA片段两端加上接头，PCR富集并进行双末端（pair-end）测序，测序读长一般为100nt，也可以为50nt。

通常，回收300bp到500bp的DNA片段进行测序，这样PE测序的左端由于酶切的原因是对齐的。另一端测序片段由于物理打断时DNA分子长度的随机性，片段之间就会存在着overlap关系，可以进行局部组装，获得较长（200~300bp）的DNA片段，并在组装结果上查找个体间的长度多态性位点，即Indel标记。

低质量阈值由具体测序技术及测序环境而定，例如设定为单碱基测序质量低于20；测序序列中测序结果不确定的碱基（如Illumina GA测序结果中的N）个数超过整条测序序列碱基个数的10%则认为是不合格序列；除样本接头序列外，与其它实验引入的外源序列比对，如各种接头序列。若序列中存在外源序列则认为是不合格序列；在酶切一端测序序列中，若起始的几个碱基不是酶切末端序列则过滤掉（如限制性内切酶Ecor1，测序序列开头若不是“AATTC”则过滤掉整个测序序列）。

2）根据测序个体基因组酶切一端的测序序列，利用序列的全同性生成每个个体堆的信息。例如，将每个个体过滤后的酶切一端的测序序列信息作为哈希的键，哈希的值指向一个链表，用于存放另一端的序列信息，并计算测序深度信息。可用任何一种编程语言实现该过程。具体过程如图4所示。堆（Stack）中酶切一端的序列信息可以以图5的方式保存，在图5中，第一列表示的是酶切一端的序列信息；第二列表示的是酶切一端序列被测序的次数，即测序深度信息；第三列是该堆的ID，用于唯一确定堆。堆（Stack）中另一端的序列信息可以以图6的方式保存，类FASTA格式。在图6中，大于号（>）之后的是堆的ID，用于唯一确定堆。大于号（>）中的内容为堆中另一端的所有测序序列信息。

3）过滤掉酶切一端序列测序深度为1的结果（成对过滤）。即过滤掉图5中第二列为1的堆，同时过滤图6中对应的数据。深度为1的结果通常是由测序错误导致的，过滤掉深度为1的测序序列信息，进一步减少由测序错误所带来的分析困难。

所述的不容许空隙的两两比对是指比对的时候不容许开空位。即不考虑开空位比对上的情况，例如以下两条序列的比对结果就不满足不容许空隙的两两比对条件：

序列1：AATTCATCGAC。

序列2：AA CATCGTC。

比对的时候容许的错配数随测序的长度来定，例如在测序长度小于50nt的情况下，容许的错配数为1，长度在100nt的情况下，容许的错配数为2。具体地，两条序列之间只有一个碱基不相同，则这两条序列归为一类。如果A序列与B序列之间只有一个碱基不相同，而B与C之间只有另外一个碱基不相同，则三条序列归为一类，以此类推，通过所有测序序列之间的比对，可以将所有满足比对条件的测序序列进行聚类。挑选出聚类结果中只有一个堆和两个堆的聚类结果。其中只有一个堆的聚类结果表明在酶切一端测序片段上不存在杂合位点，只有两个堆的聚类结果表明在酶切一端测序片段上存在杂合位点，一般情况下这个杂合位点不会处于重复区域。对于那些堆的个数超过两个的聚类结果，通常由于酶切一端测序序列处于基因组的重复区域所造成的，因此这些聚类结果将会被过滤掉。使用的比对软件可以是任何一款序列比对软件，如blast、blat等。计算堆的个数为一和二的聚类结果的总深度，并进一步过滤掉低深度和高深度的聚类结果。低深度的阈值通常为平均测序深度的四分之一，高深度的阈值通常为平均测序深度的两倍。按以上步骤对个体内的数据进行过滤和整合，重新生成如图5和图6的结果。酶切一端测序序列存在杂合位点的另一端信息将会被整合到一起，而酶切一端的测序序列信息将会用一致性序列来表示（如某个位点为A/T杂合，则用字母W表示），并累加深度，重命名ID以利于后续的分析。这样，处理之后的数据将更利于另一端数据的组装。过程如图8所示。

5）在两个个体内部，对每个堆的另一端数据进行局部组装，采用的组装可以是任何一款组装软件，如基于重叠关系的组装软件phrap，如基于De Bruijn graph算法的组装软件SOAPDenovo。利用重叠关系进行组装的时候通常不容许空位的存在，其他参数可采用软件默认设置。利用De Brui jn graph算法进行组装的时候，通常kmer的大小要在21到31之间，其他参数可采用软件默认设置。参数的意义和用法查看相关使用说明。

6）利用两个个体酶切一端的测序序列信息将两个个体堆的信息相互进行两两对齐，即在个体A和个体B中，个体A的某个堆能够和个体B的某个堆对齐，当且仅当两个个体堆中的酶切一端的测序序列完全相同。对能够对齐的堆，两个个体之间的另一端的组装结果序列相互进行比对，来寻找Indel位点信息。比对的时候，两个个体之间容许的空位数为一、并且公共区域内不存在错配。使用的比对软件可以是任何一款序列比对软件，如blast、blat等。过程如图8所示。通过以上的步骤，将会得到两个个体之间，高可信度的Indel位点信息，还得到了在Indel位点周围的侧翼序列信息。

7）最后可在Indel位点周围的侧翼序列上设计引物，以便于后续的大规模实验应用。如构建遗传图谱，个体基因分型等。如图9所示。

实施例数据：

瓠瓜F2群体两个亲本的RAD-PE测序数据。（说明：父本和母本杂交生成的后代叫F1，F1自交生成的后代叫F2）

材料来源：浙江省农业科学院。

实施例具体操作流程：

将两个亲本RAD-PE测序得到的测序数据，根据测序质量值，N的含量，以及是否含有酶切末端序列进行过滤，去除不合格的测序序列，得到的有效数据统计如表1所示。

表1、瓠瓜RAD-PE测序有效数据统计

根据测序个体基因组酶切一端的测序序列，利用序列的全同性分别生成两个亲本堆的信息。过滤掉酶切一端序列测序深度为1的结果。图10示出RAD序列测序深度分布图，结果统计如表2所示。

瓠瓜样品	RAD序列数目	RAD序列平均测序深度
			父本	198,367	15
母本	171,266	12

表2、瓠瓜RAD测序统计

两个亲本内分别将酶切一端的测序序列数据进行不容许空隙的两两比对，对堆进行聚类以确定个体内在酶切一端序列上的杂合SNP信息。并根据比对结果和深度信息进行过滤和整合。

在两个亲本内部，对每个堆的另一端数据进行局部组装，并根据酶切一端的序列信息将酶切一端序列进行对齐，另一端进行比对查找Indel位点信息。然后在存在Indel位点周围的侧翼序列设计PCR引物。综上，通过以上步骤的处理，瓠瓜父本和母本两个个体中，总共找到了658个Indel位点，并且583个位点可以设计出PCR引物用于后续的分析。该实施例中，PCR引物将用于后代的基因分型。

Claims

1.基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：其步骤如下：

1）在获得RAD高通量测序技术的测序结果后，对RAD双末端测序序列进行过滤以去除不合格的测序序列；

2）根据测序个体基因组酶切一端的测序序列，利用序列的全同性生成每个个体堆的信息；将每个个体过滤后的酶切一端的测序序列信息作为哈希的键，哈希的值指向一个链表，用于存放另一端的序列信息，并计算测序深度信息；

3）过滤掉酶切一端序列测序深度为1的结果；

4）两个个体内分别将酶切一端的测序序列数据进行不容许空隙的两两比对，对堆进行聚类以确定个体内在酶切一端序列上的杂合SNP信息；

5）在两个个体内部，对每个堆的另一端数据进行局部组装；

6）利用两个个体酶切一端的测序序列信息将两个个体堆的信息相互进行两两对齐，即在个体A和个体B中，个体A的某个堆能够和个体B的某个堆对齐，当且仅当两个个体堆中的酶切一端的测序序列完全相同；对能够对齐的堆，两个个体之间的另一端的组装结果序列相互进行比对，来寻找Indel位点信息；进而得到两个个体之间高可信度的Indel位点信息，还得到在Indel位点周围的侧翼序列信息；

7）最后在Indel位点周围的侧翼序列上设计引物，应用于后续的大规模实验。

2.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤1）中，RAD高通量测序技术为Illumina GA测序技术，或现有的其他高通量测序技术。

3.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤1）中，所述的不合格的测序序列为测序质量低于预定的低质量阈值的碱基个数超过整条序列碱基个数的50%的序列。

4.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤3）中，过滤掉酶切一端序列测序深度为1结果的方式为成对过滤。

5.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤4）中，只有一个堆的聚类结果表明在酶切一端测序片段上不存在杂合位点；只有两个堆的聚类结果表明在酶切一端测序片段上存在杂合位点；堆的个数超过两个的聚类结果被过滤掉。

6.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤4）中，计算堆的个数为一和二的聚类结果的总深度，并进一步过滤掉低深度和高深度的聚类结果；其中，低深度的阈值为平均测序深度的四分之一，高深度的阈值为平均测序深度的两倍。

7.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤5）中，组装选自基于重叠关系的组装软件phrap或基于De Brui jn graph算法的组装软件SOAPDenovo；利用重叠关系进行组装的时候不容许空位的存在；利用De Brui jn graph算法进行组装的时候，kmer的大小要在21到31之间。

8.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤7）中，比对的时候，两个个体之间容许的空位数为一、并且公共区域内不存在错配。

9.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：所使用的比对软件选自blast或blat序列比对软件。

10.根据权利要求1所述的基于酶切建库双末端测序的长度多态性标记的引物设计开发方法，其特征在于：步骤7）中，应用于构建遗传图谱或个体基因分型。