一种物种特异性内源性条形码的搜索方法及其在多样本混合测序中的应用

技术领域

本发明属于基因测序领域，尤其是一种物种特异性内源性条形码的搜索方法及其在多样本混合测序中的应用。

背景技术

近年来高通量测序技术迅猛发展，应用范围不断拓展，经常需要进行大量样本中特定DNA区域的测序工作。为了提高多样本的并行测序能力，高通量测序平台基本都提供了物理分割的测序通道，但是并行处理能力受限于通道数目，仍然无法满足大量样本的并行测序需求，因此多样本混合测序的实验方案应运而生。目前混合测序的主要方法包括DNA条形码标记和重叠混合测序。

DNA条形码(又称DNA Tag)，即一段长约数个碱基的人工DNA条形码序列，通过测序前的样本制备过程，利用PCR或连接反应，插入到对应样本序列中，并且每个DNA条形码序列唯一对应一个测序样本。之后将标记过的样本混合测序，通过获取的DNA条形码信息确定所测DNA片段的样本归属，从而实现在一次测序过程中同时对多个样本进行并行分析。

目前DNA条形码应用过程中，有两项操作是比较费时费力的：一是DNA条形码需要体外合成，并且每个样本对应的条形码要具有唯一性，因此有多少样本混合，就要体外合成多少条DNA条形码；二是在条形码合成后，需要连接对应的样本，此过程也是需要逐个样本添加。当样本数目比较多时，上述两项操作无疑需要耗费大量的人力物力。

生物内源性条形码是指生物体内能够用于辨识该物种、物种内保守并且物种间保持足够变异DNA短片段。类似于“通用产品编码”，即生活中常见的“条形码”的设计理念，生物条形码尝试用一小段DNA序列提供物种的鉴定信息，这种物种特异性的生物条形码已成为生物分类学家的一种重要研究工具，大量应用于物种鉴定、发现隐藏物种或生物多样性研究。

生物条形码来源于生物体内部，不用体外合成，并且具有很好的物种特异性，因此我们考虑在对来自不同物种的样本混合测序时，是否有可能将其作为DNA条形码标记特定样本，再进行混合测序呢？生物学分类用的条形码长度一般为几百至上千个碱基，而用于高通量混合测序的样本标记条形码受限于有效读长和测 序成本，在能够有效区分样本的前提下应当尽可能短。目前已知的各种生物条形码是否存在短的并仍可以有效区分来自不同物种的样本？即使存在这样一段区域，其两端序列是否足够保守到易于截取共同的内源性条形码DNA？这就要求一种通用有效的搜索方法，只取生物条形码的部分区域用于标记混合样本。

找到样本特异性的短的条形码序列区域时，如何简单的将其与测序目标序列进行连接呢？重叠延伸PCR(Overlap Extension PCR)技术，或称融合PCR(Fusion PCR)，利用具有碱基互补末端的引物，使PCR产物形成重叠链，从而在随后的扩增反应中通过重叠链的延伸，将不同来源的扩增片段拼接起来。利用此技术可实现多个样本一步反应内同时扩增并连接各自条形码和待测目标序列，简化了先提取测序目标序列、再逐个连接体外合成条形码的实验过程。

发明内容

发明目的：针对目前测序过程中应用外源性条形码，需要逐条合成且逐个与测序目标序列连接的实际情况，本发明提供一种物种特异性内源性条形码的搜索方法及其在多样本混合测序中的应用，通过搜索生物体内部物种特异性的短的序列片段，并在提取测序目标序列时同时进行连接，实现多个样本的同时同步处理，优化文库制备过程，提高测序效率。

技术方案：为实现上述技术目的，本发明提出了一种物种特异性内源性条形码的搜索方法，包括如下步骤：

(1)确定内源性条形码所在的候选基因组序列：根据待测序样本的特征，选择分类学上广泛应用的相应物种的生物条形码作为搜索范围，从而确定内源性条形码序列在全基因组上的一个大范围区域，通过在此范围内搜索，得到尽可能短的基因组序列作为内源性条形码所对应的区域。通常，动物界可选用线粒体基因编码的3种细胞色素氧化酶亚基之一的COI基因，植物界可选用编码叶绿体的matK、rbcL等基因，细菌可选择编码16s rRNA或编码线粒体功能性蛋白cpn60的基因；

(2)搜集候选基因组序列：搜集并下载与待测序样本进化关系较近物种的已测序候选基因组序列；

(3)比对搜集的候选基因组序列：使得所有序列长度一致并对齐，以便滑动窗口在此范围内逐个碱基扫描，寻找尽可能短的样本特异性序列对应的区域；

(4)设定条形码长度参数：内源性条形码的长度与待测序的样本数量、样本间的进化关系有关，即样本数量越多，则能够保证每个样本特异性的内源性条形码的长度就越长；样本间进化关系越近，则各自内源性条形码区域之间的序列变异度就越小，因此内源性条形码序列的长度越长。所以在较大数量或者进化关系较近的样本混合测序使用本方法时，有可能找到的区分不同样本的内源性条形码太长，而太长的条形码占用测序空间，增加测序成本。因此考虑实际测序长度和目标DNA长度，本方法需设定条形码长度参数，根据内源性条形码对实际测序长度所占的比例算得，默认值为20％，即条形码长度上限。若目标DNA长度比较短，内源性条形码长度可以适当增加；反之则应适当降低；

(5)计算滑动窗口内序列的变异度和窗口两侧序列的保守度：初始化窗宽，并以条形码长度参数为最大窗口宽度，进行以下循环计算，即以固定宽度的滑动窗口在对齐后的序列区域逐位滑动，计算每个窗口内序列的变异度和窗口两侧指定长度序列的保守度，然后逐渐扩大窗口宽度，直至找到满足要求的高变异和高保守区域，或达到滑动窗口宽度的上限。其中，变异度表征当前滑动窗口内序列的变异程度，变异度要尽可能大，以便区分不同物种，成为真正的条形码。变异度是通过当前滑动窗口中和其余序列均不同的序列占序列总数的比例来表示，具体定义如下：

设序列集A为所有样本的基因组序列在对应滑动窗口下的等长序列集合，特异性序列集合B为集合A中某序列和其它序列至少有一个碱基以上不同的序列集合，则

其中card(X)为集合X中元素的个数；

保守度表征当前滑动窗口两侧指定长度序列的保守程度，保守度要尽可能高，以便于设计引物扩增内源性条形码序列。滑动窗口的两侧序列保守度要分别计算，均要保守到同一对引物能够与所有样本的两侧序列结合并扩增。通过计算序列之间的海明距离(hamming distance)来计算窗口两侧的保守度，保守度定义如下：

设序列集C为滑动窗口某一侧指定长度的基因组序列集合，设C中数目最多的相同序列为a，C中序列y组成集合D＝{y|hamming(y，a)≤3)}，

其中hamming(y，a)表示序列y与序列a之间的海明距离，两侧的保守度均要大于给定值才认为当前滑动窗口区域内序列能够被同一引物扩增；

(6)根据滑动窗口扫描计算结果，确定内源性条形码：选择的滑动窗口内序列的变异度为100％、窗口两侧序列的保守度也为100％的序列作为物种特异性的内源性条形码。为了能够准确辨别不同样本，选择的滑动窗口内序列的变异度为100％，即窗口内每一条序列都至少有一个碱基与其他序列不同；同时为了能够扩增所有样本的条形码，要求窗口两侧序列的保守度也为100％，即所有的条形码序列都能够使用共同引物进行扩增。若滑动窗宽达到设定上限仍无法找到满足要求的内源性条形码，则停止搜索，认为在设定的当前参数条件下无法找到合适的内源性条形码，需要适当提高条形码长度参数，或者将样本分组分别混合测序、通过减少每组样本数量找到满足要求的内源性条形码。

本发明进一步提出了上述的搜索方法在多样本混合测序方法中的应用。

具体地，上述应用包括如下过程：

过程一：利用上述搜索方法得到物种特异性的内源性条形码所在的基因组区域，从而确定物种特异性的内源性条形码及其扩增引物；

过程二：扩增并连接内源性条形码与目标DNA序列，上机测序：对于每个样本，基于重叠延伸PCR技术，针对内源性条形码与目标DNA序列两端的保守序列分别设计对应的引物，同时扩增出内源性条形码和目标DNA两段序列，再通过设计引物时的碱基互补序列连接内源性条形码和目标DNA，使得每个样本都形成目标DNA序列和对应的内源性条形码的连接片段，然后混合各样本的连接片段并添加测序接头，构成测序文库，送到DNA测序仪实际测序；

过程三：判断测序片段的样本来源：测序完成后对测序结果进行分析，根据各个样本中内源性条形码的特征，溯源测序结果中的序列片段来自于哪个样本。

具体地，过程二通过如下步骤实现：

(1)分别设计内源性条形码序列和待测目标序列的扩增引物：根据内源性条形码所在的基因组区域，设计引物扩增；对于待测目标序列，按照常规方法或引物设计软件在待测目标序列两侧寻找保守区域并设计引物进行扩增，然后根据重叠延伸PCR技术的原理，内源性条形码序列的3’端引物和目标DNA序列的5’端引物除了要与各自目标区域互补外，还需各自延长15～25bp的互补区域；

(2)所有样本同时扩增并连接各自内源性条形码与目标DNA序列：对于每一个样本，第一轮PCR反应同时扩增出条形码和目标DNA两段序列，第二轮PCR通过设计引物时的互补序列连接条形码和目标DNA，即每个样本都形成目标DNA序列和对应的内源性条形码的连接片段。

本发明的思路概括如下：对来自不同物种的多个样本的某一段DNA区域测序时，首先找到另外一段DNA区域作为内源性条形码，该区域内的序列具备样本特异性且易于扩增，然后在扩增待测目标DNA序列时，通过重叠延伸PCR方法将各样本的内源性条形码序列与目标DNA序列连接，再混合后添加测序接头，之后上机测序。通过分析测序结果，根据各个样本中内源性条形码的序列特征，区分目标DNA序列来自于哪个样本。

有益效果：与现有技术相比，本发明具有如下优点：

(1)根据本发明提出的利用内源性条形码标记样本的混合测序方法，能够充分利用生物体天然的种属序列特性作为样本标记，避免逐条体外合成DNA条形码，并且利用重叠延伸PCR技术，实现了内源性条形码与目标序列的同时扩增和连接，优化了文库制备过程，提供了一种全新的高通量混合测序思路。

(2)本发明提出了基于基因组序列物种特异性搜索满足要求的内源性DNA条形码的方法，为内源性DNA条形码的实际混合测序应用提供了生物信息学的分析基础。由于基因组序列物种特异性的生物条形码本身具备物种内保守的特性，将其用于多样本测序时，各样本尽可能属于不同物种或属。目前的搜索范围仅限于生物分类学上的生物条形码区域，但是在全基因组范围上完全有可能存在个体特异性的高变区域，从而形成其它内源性条形码。随着测序技术的普及，被测序的物种越来越多，内源性条形码的搜索原理能够应用于更广的基因组范围。随着高通量测序技术的进步，测序读长会不断增加，第三代单分子测序技术的发展甚至有望打破读长的限制，基因组内源性条形码在将来高通量测序中一定会有 更大的应用空间。

附图说明

图1为本方法的原理示意图；

图2为搜索内源性条形码的程序流程图；

图3为使用滑动窗口搜索生物内源性条形码示意图；

图4为重叠延伸PCR示意图；

图5为45bp的滑动窗口在39个同科不同属样本的16s rRNA基因序列中寻找高可变区的结果图，其中虚线框为16s rRNAgene的V3高可变区；

图6为39个同科不同属样本的16s rRNA基因序列的最大变异度随滑动窗口宽度变化的结果图。

具体实施方式

本发明提出了一种利用内源性条形码标记样本的混合测序方法(如图1所示)，包括以下过程：

过程一，确定物种特异性条形码。搜索能够区分来自不同物种的所有样本、并且易于扩增的内源性条形码区域(方法流程如图2所示)。

1.1确定内源性条形码所在的候选基因组序列：针对来自不同物种样本的混合测序，确定其内源性条形码序列在全基因组上的一个大范围区域，通过在此范围内搜索，得到尽可能短的基因组序列作为内源性条形码所对应的区域。根据待测序样本的特征，可选择分类学上广泛应用的相应物种的生物条形码作为搜索范围。动物界可选用线粒体基因编码的3种细胞色素氧化酶亚基之一的COI基因，植物界可选用编码叶绿体的matK、rbcL等基因，细菌可选择编码16s rRNA或编码线粒体功能性蛋白cpn60的基因。

1.2搜集候选基因组序列：搜集并下载与待测序样本进化关系较近物种的已测序候选基因组序列。

1.3比对搜集的候选基因组序列并整理：通过对比和对齐候选基因组序列，以便滑动窗口在此范围内逐个碱基扫描，寻找尽可能短的样本特异性序列对应的区域(如图3所示)。

1.4设定条形码长度参数：考虑实际测序长度和目标DNA长度，本方法需设定条形码长度参数，根据内源性条形码对实际测序长度所占的比例算得，默认 值为20％，即条形码长度上限。若目标DNA长度比较短，内源性条形码长度可以适当增加；反之则应适当降低。

1.5计算滑动窗口内序列的变异度和窗口两侧序列的保守度：初始化窗宽，并以条形码长度参数为最大窗口宽度，进行以下循环计算，即以固定宽度的滑动窗口在对齐后的序列区域逐位滑动，计算每个窗口内序列的变异度和窗口两侧指定长度序列的保守度，然后逐渐扩大窗口宽度，直至找到满足要求的高变异和高保守区域，或达到滑动窗口宽度的上限。其中，变异度定义如下：

设序列集A为所有样本的基因组序列在对应滑动窗口下的等长序列集合，特异性序列集合B为集合A中某序列和其它序列至少有一个碱基以上不同的序列集合，则

其中card(X)为集合X中元素的个数。

保守度表征当前滑动窗口两侧指定长度序列的保守程度，保守度要尽可能高，以便于设计引物扩增内源性条形码序列。滑动窗口的两侧序列保守度要分别计算，均要保守到同一对引物能够与所有样本的两侧序列结合并扩增。通过计算序列之间的海明距离(hamming distance)来计算窗口两侧的保守度，具体表示如下：

设序列集C为滑动窗口某一侧指定长度的基因组序列集合，设C中数目最多的相同序列为a，C中序列y组成集合D＝{y|hamming(y，a)≤3)}，

其中hamming(y，a)表示序列y与序列a之间的海明距离。两侧的保守度均要大于给定值才认为当前滑动窗口区域内序列能够被同一引物扩增。

1.6根据滑动窗口扫描计算结果，确定内源性条形码：选择的滑动窗口内序列的变异度为100％、窗口两侧序列的保守度也为100％的序列作为物种特异性的内源性条形码。若滑动窗宽达到设定上限仍无法找到满足要求的内源性条形码，则停止搜索，认为在设定的当前参数条件下无法找到合适的内源性条形码， 需要适当提高条形码长度参数，或者将样本分组分别混合测序、通过减少每组样本数量找到满足要求的内源性条形码。

过程二，扩增并连接内源性条形码与目标DNA序列，上机测序。对于每个样本，基于重叠延伸PCR技术，针对内源性条形码与目标DNA序列两端的保守序列分别设计对应的引物，在合适的实验条件下同时扩增出条形码和目标DNA两段序列，再通过设计引物时的碱基互补序列连接条形码和目标DNA(如图4所示)，即每个样本都对应形成目标DNA序列和对应的内源性条形码的连接片段。在保证每个样本目标DNA都连接了特异性内源性条形码的基础上，混合并添加测序接头，构成测序文库，送到DNA测序仪实际测序。

过程三，判断测序片段的样本来源：测序完成后对测序结果进行分析，根据各个样本中内源性条形码的特征，溯源测序结果中的序列片段来自于哪个样本。

下面通过具体的实施例详细说明本发明，本发明的保护范围并不仅仅局限于本实施方式的描述。

本实施例以来自肠杆菌科不同属的39个样本为例，描述用于混合测序的内源性条形码的搜索过程。设测序目的片段为ccmD基因，与细胞色素C生物合成蛋白质和电子传递相关，长约210bp，在NCBI(http：//www.ncbi.nlm.nih.gov)的gene数据库中能够搜到671条相关基因的序列信息，其中，大肠杆菌K-12菌株的ccmD基因号GeneID为12931490，对应的核酸序列为：ATGACCCCTGCATTTGCTTCCTGGAATGAATTTTTCGCAATGGGCGGTTACGCCTTTTTTGTCTGGCTGGCGGTGGTGATGACCGTTATTCCGCTGGTGGTTTTGGTCGTGCACTCGGTGATGCAACATCGCGCAATTCTGCGTGGCGTGGCGCAACAGCGGGCGCGTGAGGCGCGTTTACGTGCTGCGCAACAGCAGGAGGCTGCATGA，利用引物设计软件得到正向引物为：5’-GAGGCCGTAAATGACCCC，反向引物为：5’-GGCAATCCACAAGCGGT。

步骤1：搜索物种特异性的内源性条形码。

1.16s rRNA基因是细菌编码16s rRNA所对应的DNA序列，对于细菌是一种较为理想的生物内源性条形码。16s rRNA基因具有高度的保守性和特异性，同样能够较为容易的进行PCR扩增。因此我们选择在16s rRNA基因区域内进行滑动窗口扫描，搜索适合作为内源性条形码的高可变区域。

2.搜集候选生物条形码序列。在NCBI网站以关键词“16s rRNA gene”和肠杆菌科不同的属名在“Nucleotide”数据库中进行搜索，收集到39条同科不同属的基因序列。

3.在Mega软件中运用ClustalW和Muscle比对算法分别对齐39条序列，综合两种算法，得到合理的比对结果。

4.设定参数。illumina公司的主流二代测序平台目前使用双端测序的测序读长为500bp，内源性条形码的长度上限默认值为20％，即100bp，因此滑动窗口宽度上限为100bp。

5.以固定宽度的滑动窗口在对齐后的序列区域逐位滑动，并根据变异度公式计算当前窗口内序列的变异度，同时根据保守度公式计算当前窗口两端序列的保守度(如图3所示)。记录当前宽度的滑动窗口所能取得最大变异度的窗口位置，然后逐渐扩大窗口宽度，直至找到满足要求的高变异区域或达到滑动窗口宽度的上限。

6.经过滑动窗口扫描，得到39条不同属16s rRNA基因在窗宽为45bp滑动时变异度结果如图5所示，其中16s rRNA基因高可变区V3内起始位点为269、长为45bp的序列达到100％的个体特异性(45bp小于长度上限100bp)，并且两侧20bp区域内的序列保守度达到100％，意味着39个样本的此45bp区域的序列能够用同一引物扩增，显示这段区域内序列具备作为内源性条形码的潜力(大肠杆菌所在的埃希氏菌属对应的内源性条形码序列为：TACTTTCAGCGGGGAGGAAGGGAGTAAAGTTAATACCTTTGCTCA)。图6显示了滑动窗口所能达到的最大变异度随窗宽变化的情况，在滑动窗口宽度达到45bp以上，对应的最大变异度才为100％，即完全能够区分39个样本的内源性条形码至少长45bp。

步骤2，通过两轮PCR完成重叠延伸PCR，连接内源性条形码和ccmD基因(见图4)。第一轮PCR使用两对引物扩增，其中正向引物(5’-CCCAACATTTCGTGAAAGTC)和反向引物(5’-GCTGGCACGGAGTTAGC)扩增出63bp的高变区(即条形码)，其中，反向引物5’端连接10对ct碱基的重复片段(共计20个碱基)；正向引物(5’-GAGGCCGTAAATGACCCC)和反向引物(5’-GGCAATCCACAAGCGGT) 扩增ccmD基因，其中正向引物5’端连接10对ag碱基的重复片段(共计20个碱基)；根据重叠延伸PCR技术原理，第二轮PCR使用扩增条形码的正向引物(5’-CCCAACATTTCGTGAAAGTC)和扩增ccmD基因的反向引物(5’-GGCAATCCACAAGCGGT)，在PCR的实验条件下，内源性条形码序列的3’端引物的20bp碱基重复片段和目标DNA序列的5’端引物20bp碱基重复片段互补，使得条形码和ccmD基因通过互补序列连接，从而各样本的ccmD基因具有样本特异性的序列条形码标识。39个样本的重叠延伸PCR过程可以在分离的反应环境中同步进行，与外源性条形码相比大大提高了文库制备效率。

步骤3，在保证每个样本都连接了特异性内源性条形码的基础上，混合并添加测序接头，构成测序文库，然后上机测序。测序完成后对测序结果进行分析，根据各个样本中内源性条形码的特征，溯源测序结果中的序列片段来自于哪个样本。

以上详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换，这些等同变换均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。