CN111128303A - 基于已知序列确定目标物种中对应序列的方法和系统 - Google Patents

Abstract

本发明提出了基于已知序列确定目标物种中对应序列的方法和系统。方法包括：基于所述已知序列，确定所述已知序列的全部Kmer序列，以便获得种子序列Kmer序列集合；获取来自于所述目标物种的测序读段，并基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列，以便获得测序读段Kmer序列集合；确定同时存在于所述测序读段Kmer序列集合和所述种子序列Kmer序列集合的Kmer序列作为延伸起始序列；基于重叠原则，利用所述测序读段Kmer序列集合，对所述延伸起始序列进行延伸处理，以便获得所述已知序列在所述目标物种中的对应序列。

Description

基于已知序列确定目标物种中对应序列的方法和系统

技术领域

本发明涉及生物信息领域，具体地，本发明涉及基于已知序列确定目标物种中对应序列的方法和系统。

背景技术

自“大数据”被提出以来，生物大数据如火如荼地发展。通过高通量测序，人们已经积累了庞大的数据(～30Pb/年)；实现构建多种数据库，如10KP、OneKP等(https:// db.cngb.org/)。通过移动互联网，人们也从各大数据库中获得了海量的数据，如NCBI、EMBL等。这些数据的应用，将为疾病诊断、分型、医药开发、精准育种等提供新方向以及新工具。

对于生物大数据，包括各种情况：所涉及到的样品多；每个样品的数据量多；每个样品的测序数据有差异：读长、质量、数据量等。面对这种现状，不少研究者表示，这些海量数据可能会淹没现有的分析渠道，并提出前所未有的“高”要求：假如在全基因组研究中，关注的仅仅只是整个基因组中的外显子部分，它能够将需要分析的数据量减少到原来的1％，但即使在这种情况下，每年产出的数据量仍可达4000万Gb。而现在的研究，往往聚集于识别个人基因组中可扰乱基因功能的“小错误”，即所谓单核苷酸突变(single-nucleotidevariants,SNPs)，平均下来，依然有近13000个之多，而仅仅2％被预知可影响相应蛋白的变化。

那么，又假如已经确定某类疾病的具体致病基因，或者说在关心已经被预知可影响相应蛋白的位点的同时，又如何从海量的个人全基因组数据中高效精准地确定每个个体在该基因是否有变异，以实现类似于“精准医学”的概念？

为此，设计出一套具有广泛应用性，能够实现从全基因组大数据中高效、精准捕获长片段目标区域的方法就显得尤为必要。

发明内容

本申请是基于发明人对以下事实和问题的发现和认识作出的：

基于全基因组大数据，获得目标区域相应的完整序列，目前，其实现的方案主要有：

1)对全基因组的所有数据，进行全局组装，并对最后的结果进行筛选，获得相应的区域。

2)选取其近缘物种相应的区域作为参考序列，将全基因组数据比对到参考序列，并筛选相应的测序数据，进行局部组装得到相应的区域。

3)在其近缘物种相应的区域以及上下游区域设计引物进行PCR扩增，通过测序，并进行全局组装得到相应的区域。

基于上述现有技术，发明人发现，针对全基因组数据进行全局组装，因为涉及到较大的起始数据量(可能还需涉及到多种文库，才能得到较好的结果)，导致在实现的过程中所需的内存大，时间长，中间文件存储大；受其它区域的干扰，影响目标区域的组装，如高度相似的区域；结果中，除了目标区域的序列之外，其它的结果对目标区域的分析不起作用；对有用结果不能做到有效合理的筛选，也无法实现自动化操作。另外，在大数据背景下，所涉及到物种数目庞大，无法一一进行组装，如10KP，其存储和计算将面临巨大的挑战。

针对全基因组数据进行局部组装，如果没有近缘物种或者物种目标区域变化较大，通过比对的方式，只能获取高度相似区域相应的测序数据，而差异较大的区域却无法实现，在这种情况下，目标区域有可能不是完整的。

同样，依赖参考序列设计引物时，有可能因为差异较大而无法完成PCR扩增。另外，挑选相应的测序数据时，主要通过比对的方式，对于大数据，在这一步将消耗大量的时间。

本发明旨在至少在一定程度上解决相关技术中的技术问题之一。

为此，本发明提出了一种能够从全基因组大数据中高效、精准捕获长片段目标区域(目前可实现～200Kb)的方法，并从结果中可获得大量的信息供给进行下一步的判断，包括：确定物种是否存在该目标区域；确定物种目标区域是否存在变异(SNPs、INDEL等)；确定目标区域或局部区域的拷贝数(依赖测序深度)；确定基因的可变剪切类型及其表达量(针对全转录组数据，并在本发明的结果中捕获到多条序列)。

在本发明的第一方面，本发明提出了一种基于已知序列确定目标物种中对应序列的方法。根据本发明的实施例，所述方法包括：(1)基于所述已知序列，确定所述已知序列的全部Kmer序列，以便获得种子序列Kmer序列集合；(2)获取来自于所述目标物种的测序读段，并基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列，以便获得测序读段Kmer序列集合；(3)确定同时存在于所述测序读段Kmer序列集合和所述种子序列Kmer序列集合的Kmer序列作为延伸起始序列；(4)基于重叠原则，利用所述测序读段Kmer序列集合，对所述延伸起始序列进行延伸处理，以便获得所述已知序列在所述目标物种中的对应序列。根据本发明实施例的方法能够从全基因组大数据中，基于已知序列，高效、精准捕获目标物种中对应序列的长片段目标区域(目前可实现～200Kb)，并从结果中可获得大量的信息供给进行下一步的判断，为大规模动植物进化和遗传学研究提供了强有力的技术支持。

在本发明的第二方面，本发明提出了一种基于已知序列确定目标物种中对应序列的系统。根据本发明的实施例，所述系统包括：种子序列Kmer序列集合获得装置，所述种子序列Kmer序列集合获得装置用于基于所述已知序列，确定所述已知序列的全部Kmer序列，以便获得种子序列Kmer序列集合；测序读段Kmer序列集合获得装置，所述测序读段Kmer序列集合获得装置用于获取来自于所述目标物种的测序读段，并基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列，以便获得测序读段Kmer序列集合；确定延伸起始序列装置，所述确定延伸起始序列装置与所述种子序列Kmer序列集合获得装置和所述测序读段Kmer序列集合获得装置相连，用于确定同时存在于所述测序读段Kmer序列集合和所述种子序列Kmer序列集合的Kmer序列作为延伸起始序列；目标物种中的对应序列获得装置，所述目标物种中的对应序列获得装置与所述确定延伸起始序列装置相连，用于基于重叠原则，利用所述测序读段Kmer序列集合，对所述延伸起始序列进行延伸处理，以便获得所述已知序列在所述目标物种中的对应序列。根据本发明实施例的系统适于执行根据本发明实施例的基于已知序列确定目标物种中对应序列的方法，根据本发明实施例的系统可用于从全基因组大数据中，基于已知序列，高效、精准捕获目标物种中对应序列的长片段目标区域(目前可实现～200Kb)，并从结果中可获得大量的信息供给进行下一步的判断，为大规模动植物进化和遗传学研究提供了强有力的技术支持。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是根据本发明实施例的基于已知序列确定目标物种中对应序列的方法的实施流程图；

图2是根据本发明实施例的基于已知序列确定目标物种中对应序列的系统的结构示意图；

图3是根据本发明再一实施例的基于已知序列确定目标物种中对应序列的系统的结构示意图；

图4是根据本发明实施例的两个水稻OLE1基因结构比较图；

图5是根据本发明实施例的野生水稻目标区域覆盖度的结果图；以及

图6是根据本发明实施例的野生水稻ELMO基因可变剪切形式的结果图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

基于已知序列确定目标物种中对应序列的方法

在本发明的第一方面，本发明提出了一种基于已知序列确定目标物种中对应序列的方法。根据本发明的实施例，所述方法包括：(1)基于所述已知序列，确定所述已知序列的全部Kmer序列，以便获得种子序列Kmer序列集合；(2)获取来自于所述目标物种的测序读段，并基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列，以便获得测序读段Kmer序列集合；(3)确定同时存在于所述测序读段Kmer序列集合和所述种子序列Kmer序列集合的Kmer序列作为延伸起始序列；(4)基于重叠原则，利用所述测序读段Kmer序列集合，对所述延伸起始序列进行延伸处理，以便获得所述已知序列在所述目标物种中的对应序列。根据本发明实施例的方法能够从全基因组大数据中，基于已知序列，高效、精准捕获目标物种中对应序列的长片段目标区域(目前可实现～200Kb)，并从结果中可获得大量的信息供给进行下一步的判断，为大规模动植物进化和遗传学研究提供了强有力的技术支持。

根据本发明的实施例，所述已知序列来自于所述目标物种的亲缘物种。进而利用根据本发明实施例的方法捕获对应序列的成功率、精准率进一步提高。

根据本发明的实施例，步骤(1)和步骤(2)中，分别采用相同的序列长度K值和间距D值，确定所述已知序列的全部Kmer序列和所述测序读段的Kmer序列，所述K值为27-39的整数，优选地，K值为31，所述D值为1。

根据本发明的实施例，所述测序读段是双向测序系统，优选地，来自BGI SEQ 500的双向测序数据。进而如果采用低深度延伸，如测序结果中只需约30X，即可达到约99.99％的准确性。

根据本发明的实施例，在步骤(2)中，仅针对双向测序读段中的正向测序读段构建所述测序读段Kmer序列集合。

根据本发明的实施例，在步骤(2)中，在基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列之前，预先对所述测序读段的至少一部分进行优化处理。

根据本发明的实施例，所述优化处理包括删除所述测序读段起端和末端至少之一的1～10个碱基。进而可尽量减少高错误率的测序数据，获得高质量的测序读段数据，比如，要求测序读段的Q30达到85％以上。

根据本发明的实施例，在步骤(4)中，当所述延伸处理的延伸产物长度达到预定长度的1.5～2.5倍时，停止所述延伸处理，优选地，当所述延伸处理的延伸产物长度达到预定长度的2倍时，停止所述延伸处理。进而可进一步确保所获得的目标物种中对应序列的完整性。

根据本发明的实施例，在步骤(3)中，确定多个延伸起始序列，在步骤(4)中分别针对所述多个延伸起始序列分别进行所述延伸处理，以便获得多个延伸产物，并通过下列步骤从所述多个延伸产物中选择最终序列作为所述已知序列在所述目标物种中的对应序列：分别确定所述多个延伸产物的每一个中所包含的种子序列Kmer序列数目；选择包含所述种子序列Kmer序列数目最高的所述延伸产物作为所述对应序列。

根据本发明的实施例，所述方法进一步包括：(5)确定步骤(4)中进行所述延伸处理所采用的测序读段Kmer序列所对应的候选测序读段；(6)基于所述候选测序读段对所述对应序列进行修正处理。

根据本发明的实施例，所述修正处理包括：基于所述候选测序读段与所述对应序列的比对，针对所述对应序列的至少一个位点，确定所述至少一个位点的优势碱基，并利用所述优势碱基对所述对应序列进行修正。

根据本发明的实施例，当所述对应序列的至少一个位点上存在优势碱基时，则说明对应序列的至少一个位点上存在测序错误，利用优势碱基，对所述测序错误进行修正。

根据本发明的实施例，所述修正处理包括：确定双向测序读段中成对测序读段在所述对应序列上的间距，如果所述间距与预定的插入片段长度差异超过10％，则判定所述对应序列为错误序列。

根据本发明的实施例，对于不存在所述优势碱基且测序深度差异低于2倍，则将所述位点标记为SNP突变。

根据本发明的实施例，对于不存在所述优势碱基且测序深度差异不低于2倍，则将所述位点标记为CNV突变。

为了便于理解，发明人将根据本发明实施例的基于已知序列确定目标物种中对应序列的方法的发明构思以及具体实施步骤进行了进一步详细阐述。

本发明的目的是为了实现：对于大量的全基因组或全转录组数据(一个物种的多个样品，多个物种的同一组织，多物种等)，如何高效、精准得获得目标物种中对应序列的长片段目标区域的完整序列，进而通过对其结果进行分析，可对相应的目标物种或个体进行判断：确定是否存在目标区域；确定目标区域是否存在变异(SNPs、INDEL等)；确定目标区域或局部区域是否存在多个拷贝数(通过测序深度实现)；确定多种可变剪切及其表达量(针对全转录组数据，是否生成多种可能序列)。因此，本发明的实现为大规模动植物进化和遗传学研究提供了技术支持。

本发明实施例所涉及到概念或参数，具体如下：

Seed(参数，string)，种子序列。即目标区域的参考序列，可为近缘或远缘物种的任意一条序列(长度应大于K，包括编码区，内含子，或基因间区。测试结果中，编码区效果最好)。本发明利用该序列获取目标区域可能的Kmer，得到Seed-Kmer哈希列表(记录了种子seed序列中所有可能的Kmer序列)。

Kmer，长度为K的短序列。如长度为L的Seed，从起始端开始，每次移动一个碱基，则可得到(L-K+1)个Kmer。取值范围27-39，优选为31。

K和D(参数，int)，Kmer的长度和相邻Kmer间距。D主要应用Seed-Kmer哈希列表的生成，使产生的Kmer差异足够大。优选D为1。

Hash，哈希。编程语言中的一种数据类型，为Key/Value对的集合。每个Key都是唯一存在的，Value可为多种形式。

Kmer哈希列表存储了参考种子所有的Kmer序列，即长度为K的序列。判断测序数据read1是否存在相应的参考种子Kmer序列，这里要求Kmer序列完全相同，存在时，依赖K-1个重叠进行延伸，并记录相应的Paired end信息，以供后续验证。

Hash Table，哈希列表。在本发明中用于储存Seed中Kmer种类(Seed-Kmer哈希列表)、测序数据Read1和Read2(Read和Read-Kmer哈希列表)。

Read(参数，string)，测序下机数据，包括Read1和Read2。目前，本发明适用于类似Illumina数据，且为Paired end数据，如BGISEQ-500。本发明适用于Illumina的Paired end数据，改进后，也可使用高准确性的“长”读长的Single end数据，例如sanger。这是因为，首先，本发明使用Kmer概念，要求测序数据足够低的测序错误，像高错误率的PacBio数据(～15％)不适合。其次，本发明通过Paired end的对应关系确定其真正的延伸序列(当有多种延伸方式时起作用)。如果测序的读长足够长的话，其实也不需要借助Paired end关系，如sanger(准确性高且长读长)。最后，项目数据基本来源于华大BGISEQ500的Paired end数据，也是针对实际数据来设计的。

Assembly，组装序列。本发明是基于Assembly，迭代延伸，直到满足条件，迭代停止，获得目标区域序列。初始的Assembly为通过Seed-Kmer定位到的Read1。

Consensus，一致性序列。获得Assembly之后，将Read重新定位到相应的位置，在高深度Read覆盖的情况下，进行纠正。

LEmax和REmax(参数，int)，左右总延伸序列长度最大限制。本发明认为：不同物种目标区域序列有可能在长度上存在差别，优选设置为2倍期望长度。例如，假设目标区域是2K，即Seed序列长度，如果不设置的话，优选设置为Seed序列长度的2倍，即2*2K。如果延伸成功的话，获得的序列应该是4K。正常的话，目标区域存在于4K序列中。发明人之所以将左右总延伸长度设置为Seed序列长度的2倍，是因为发明人考虑到：本发明要实现的是根据参考种子序列，获得目标物种对应序列的相应区域。如果是近缘物种，相应的区域应该和参考种子序列差别不大，包括长度、相似度。但是，假如物种相隔甚远，或者物种本身所处环境特殊而导致相应的区域发生巨大的变化。如果是因为片段丢失导致目标区域缩短或消失，本发明的方法依然可以识别；但如果是因为片段插入导致目标区域扩张，那么，延伸长度不足的话，会导致获得的目标区域残缺。前期分析过程中，设置2倍，足以确保目标区域的完整性，后续只需要再比对下，就可以确定其准确的坐标位置。

Emax(参数，int)，总延伸序列长度最大限制。该参数有区别于LEmax和REmax，主要针对有可能存在目标区域上下游存在复杂区域，导致其中一端无法延伸。

Rmax(参数，int),存储于内存的最大Read数目。本发明将Read1中起始位置上的Kmer(真实位置具体由RcutP确定)设为Key,将相应已储存的Read条数Num，即第Num条，设为Value，得到Read哈希列表。同时，以Num为Key，以Read1和Read2为值，得到Read-Kmer哈希列表。

RcutP和RcutL(参数，int),Read1和Read2截取起点和长度。由于测序数据在前后端存在较高的错误率，因此，在存储于内存的Read序列可进行适当截取，以提高准确性。当测序读段质量较高，可不截取。不截取能确保足够长的延伸长度，即全读长，能解决部分短的重复序列识别问题，但是，如果测序质量非常差，表现为起端和末端，建议进行适当截取，如截取测序读段起端和末端至少之一的1～10个碱基。

RQmin(参数，int)，Read1和Read2最小平均质量值。由于测序数据在前后端存在较高的错误率，因此，在存储于内存的Read序列可进行适当过滤，以提高准确性。

在大数据背景下，捕获目标区域的完整序列并不需要用到所有的数据，因此，本发明首要解决的技术问题是如何从大数据中获取足够的有效数据进行后续的分析，以降低使用大数据下占用大内存，大存储的问题。同时，由于涉及到过多的样品或物种数量，不能使用常规的方法实现。针对这些问题，本发明主要通过Kmer定位目标区域序列，并进行双向延伸，迭代延伸实现，具体的实施步骤(可见图1)包含以下内容：

(1)根据输入的Seed序列，需要生成的Kmer长度K以及相邻Kmer间距D，生成Seed-Kmer哈希列表；该列表的作用包括：用于识别存在于物种目标区域的Read，即存在相同的Kmer，作为初始化Assembly；无法获得初始Assembly时，本发明认为该物种或个体已完全丢失目标区域序列，程序终止；在最后的结果中，Assembly中存在数目较多的序列时，利用该表判断每条序列所含Seed Kmer的比例，最多的一条反映为真实的目标区域序列。

换句话说，即首先，通过设置K和D这两个参数，可以基于预先输入的Seed序列生成一系列Seed的Kmer，另外，通过对read进行相同的设置参数，也可以生成一系列read的Kmer；其次，将Seed的Kmer和read的Kmer进行比较，如果能找到相同的Kmer，那这个Kmer序列就作为组装序列assembly的初始序列。

这时有可能出现三种情况，第一，一个相同的Kmer也没发现，那说明seed序列对应的区域在研究的样本中已经消失了。第二，只找到了一个相同的Kmer，那直接将这个Kmer作为初始的assembly。第三，找到了多个相同的Kmer。需要进行的处理是，都让它进行延伸，当最终的结果中，其延伸的多条序列有差异，那么，通过比较所含的参考种子Kmer比例来确定哪条序列是真正的目标区域。这一步是在最后才进行处理的，在前期延伸部分不会执行，而前期也不会对存在这种情况的read进行过滤。

(2)根据输入的测序数据Read1和Read2，Reads的截取参数(RcutP和RcutL)，以及Reads的质量控制参数(RQmin)进行分析，获得较高质量的Read1和Read2(如高质量的read的Q30达到85％以上)。根据控制Reads条数的参数Rmax，开始将高质量的Reads存储到内存，得到Read哈希列表。该表同时保存了Read1和Read2的信息，因此，根据该表和延伸序列可获得每对Read间距的大小，以判断是否符合实际文库大小，对于异常的情况，其相应的延伸提前终止。

例如，一个具体的测序读段可依据下列运行方式进行判断测序质量：

第一行，包括序列信息，lane信息和坐标信息；

第二行，测序碱基；

第三行，特定符号；

第四行，包含序列的质量信息，反映序列测序错误的概率。根据公式Q＝-10log(P/1-P)，公式中的Q就是第四行每个碱基对应的字符的ASCII-33或ASCII-64得到，然后就可以推算P。根据第四行，可以计算整条read的平均错误率，如果错误率很高，那就直接过滤掉。

根据本发明的实施例，每对Read间距不是一个精确值，比如构建一个300bp的文库，这个值是一个均值，对于每对reads，它们的间距是在这个范围波动的，正常情况下+/-10％正如上面提到的，只有明显不符合才被确认删除。在本发明中，主要根据Paired end关系确定准确的延伸序列，在其它步骤中并不需要。

根据本发明的实施例，本发明采用的先延伸序列，再根据Paired end间距确定正确的延伸序列。整个过程可如下所述：

a、根据参考种子序列获取Kmer列表；

b、根据read1，识别存在与参考种子相同的Kmer，并确定初始化Assembly；

c、根据初始化Assembly序列，根据read1判断是否存在K-1个重叠进行延伸，该步骤不需要Paired end关系，出现分叉的地方，都会延伸；

d、对于延伸序列，应该c存储的Paired end关系，重新在延伸序列中定位，并确定其间距。如果其间距超过正常文库大小的+/-10％，则认为这对Paired end不是目标区域的，则该相应的延伸序列并不是真正的目标区域。

(3)在(2)确定高质量Read的基础上，以及需要生成的Kmer长度K，获得Read1起始端的第一个Kmer，生成Read-Kmer哈希列表。该列表并不包含Read2的信息，主要应用于对Assembly的双向延伸，延伸必须符合存在K-1重叠，即每次延伸一个碱基。当它无法延伸时，表明Read-Kmer哈希列表在Assembly边界区域不存在相同的Kmer可供其延伸，程序终止。另外，本发明采用低深度延伸即可，在前期测试中，结果表明只需～30X，已获得较高的准确性。测序深度受Rmax影响，该值越大，内存使用越大，运行时间越长。当达到一定范围，会出现极限值，但并没改变结果，所以并不是越大越好的。

虽然，read1和read2都可以实现延伸。在此，相当于有一个验证的过程，根据本发明实施例的方法采用的都是read1进行延伸，然后通过read1和read2的Paired end关系验证延伸部分是否符合；这种延伸的方法并不需要高深度，即使只使用read1进行延伸，但已经可以实现了；read2测序错误比read1严重，在能实现的情况下，无需引入更多的错误。

根据本发明实施例的方法，针对的是全基因组测序，正常的基因组项目中，如果只存在一个插入片段只有300bp的文库，且测序深度只有30X。那么：

a、如果使用所有数据进行整体组装，其结果基本上是不可能好的。整体组装复杂性高，时间和内存是指数分布的。

b、如果使用近缘物种作为参考，进行比对，筛选数据，局部组装。当目标物种的目标区域与参考相差比较大时，根据这种方法，其结果会有错误。

c、根据本发明的实施例，本发明尝试测试的是30X的全基组数据，但对于目标区域来说，过滤严重测序错误、重复序列等情况，最终被所使用的read不超过10X。

根据本发明的实施例，Rmax一般设置20G即可。Rmax一般用于存储测序read1和read2信息。当Rmax越大，意味着与目标区域相关的read越多，对于某个区域，当测序深度越高，并不意味着延伸越简单，相反，会增加延伸的干扰，主要表现为测序错误的增加，导致延伸出现更多的分歧。例如，假设20G时对应的目标区域已经达到20X，延伸成功。而当达到200G时，则目标区域深度达到200X，如果存在1％的测序错误，则意味着目标区域的每个位置都会存在2个测序错误。

(4)对延伸区域进行纠正，获得一致性序列，生成新的Assembly。由于测序数据存在高错误率的情况，虽然前期已经作了大量的控制，但依然部分错误被保留。在这里，发明人将Read1和Read2重新定位到Assembly相应的位置，全基因组测序在一定程度上符合均匀测序，因此，可根据每个位置上相应的高深度信息进行纠正，无法纠正部分，反映为SNP或INDEL，同时保留各种情况的相应序列。

根据本发明的实施例，此处的错误来源于测序错误，这部分是不可避免的。根据本发明的实施例，实现纠错的方法包括：将所有用于延伸assembly的kmer对应的reads与assembly进行比对，通过各位点的测序深度优势对assembly进行校正。例如assembly的第100位是T，而通过比对后，T所对应的测序深度仅1，而A对应的测序深度是20，那么就把T修改为A。

根据本发明的实施例，无法矫正是指的无论哪个碱基都没有明显的测序深度优势，例如，如果不是低频的变异，A的测序深度是10，T的测序深度为11，那么就同时保留A和T，认为在这个位点存在两个天然的SNP。而对于低频变异，较多的一个碱基作为代表。

(5)对于新的Assembly，判断其总的延伸长度是否满足要求，满足要求输出相应的序列的同时，把能锚定到Assembly的所有原始的Read1和Read2输出，以供更深一层的后续分析。相反，重新返回(3)，重新进行双向延伸，直接满足LEmax和REmax要求，迭代停止。

根据本发明的实施例，如果针对的是全转录组数据，且出现了多种延伸序列，则很大程度上为不同的可变剪切，发明人只需要借助专门的转录组数据比对软件(hisat，tophat，bowtie等)，将多种延伸序列重新比对回原来基因组，并根据深度分布判断其可能的剪切方式，估计不同可变剪切表达趋势即可。

根据本发明的实施例，这里提到的重新返回(3)，主要依据LEmax和REmax，不由其它参数控制。在延伸时，即使存在无限延伸的可能性，也不会无限延伸，因为每次延伸一定长度时(可设置为500bp)，它会暂停并判断是否已经满足LEmax和REmax。不满足时，它会继续延伸；满足时，终止。

至此，本发明的方法的实施流程结束，获得目标区域的完整序列的同时，也获得该区域所有的原始下机数据Read1和Read2，该部分对于后续分析起到关键作用：首先，它是原始的下机数据，反映物种真实的情况。发明人可以将数据重新进行比对(如BWA)到Assembly，根据Pair end关系，可验证其准确性，并统计相应的深度信息，可判断某些局部区域是否存在特殊情况，如多拷贝。其次，如果是全转录组数据，且Rmax选取足够大(足以覆盖整个数据，全转录组数据相比全基因组数据来得小，很容易达到)，那么，它反映的是真实表达水平。如果结果中，出现多种延伸序列，很大程度上为不同的可变剪切，发明人只需要借助专门的转录组数据比对软件(hisat，tophat，bowtie等)，重新比对回原来基因组，并根据深度分布可判断其可能的剪切方式，估计不同可变剪切表达趋势。

另外，由于发明人不确定物种真实的目标区域长度、分化程度等。因此，发明人设置了一个较大的LEmax和REmax，或Emax值。那么，得到的结果，应该是包含目标区域的完整序列，及其上下游序列。所以，在最后还需要通过比对软件，确定其真实的边界，截取相应的区域之后，可借助基因组比较软件，如Mummer，检测其是否存变异。

对于基因组某个目标区域来说，它只占全基因组很小的一部分(如大小为～1G基因组，其目标区域为100Kb，即～0.01％)，即大量的测序数据中其实只存在很小一部分能被所利用。针对这种情况，根据本发明实施例的方法通过以下方式实现：

1)基于Seed序列中的Kmer定位到某些Read，即目标区域的序列，并进行双向延伸(低深度延伸，测试结果中，只需～30x，即可达到～99.99％的准确性)，成功避开目标区域之外的绝大部分数据(目标区域存在的重复序列之外)，从而大量减少长时间占用大内存，存储等问题，尤其是无法确定基因组大小的物种。另外，相比全局组装结果，本发明产生的结果简单扼要，只需要简单处理即可进行判断，包括序列比较，深度分布等。

2)对于目标区域在不同的物种中存在较大差别的情况，本发明认为：即使差别再大，但始终会保留着大小为Kmer级别的部分序列，一旦被检测到，就可以实现延伸。而如果连Kmer级别的序列都不存在的话，本发明认为相应的物种已丢失相应的目标区域。

3)由于全转录组反映的是基因转录出来的结果，虽然在基因组体现为多个不连续的区域(多个外显子的基因)，但测序的数据反映的依然是连接的区域。因此，本发明也可通过全转录组捕获目标区域序列(即基因转录序列)，进而可判断基因可变剪切及表达情况。

基于已知序列确定目标物种中对应序列的系统

在本发明的第二方面，本发明提出了一种基于已知序列确定目标物种中对应序列的系统。根据本发明的实施例，参考图2，所述系统包括：

种子序列Kmer序列集合获得装置100，所述种子序列Kmer序列集合获得装置100用于基于所述已知序列，确定所述已知序列的全部Kmer序列，以便获得种子序列Kmer序列集合；

测序读段Kmer序列集合获得装置200，所述测序读段Kmer序列集合获得装置200用于获取来自于所述目标物种的测序读段，并基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列，以便获得测序读段Kmer序列集合；

确定延伸起始序列装置300，所述确定延伸起始序列装置与300所述种子序列Kmer序列集合获得装置100和所述测序读段Kmer序列集合获得装置200相连，用于确定同时存在于所述测序读段Kmer序列集合和所述种子序列Kmer序列集合的Kmer序列作为延伸起始序列；

目标物种中的对应序列获得装置400，所述目标物种中的对应序列获得装置400与所述确定延伸起始序列装置相连300，用于基于重叠原则，利用所述测序读段Kmer序列集合，对所述延伸起始序列进行延伸处理，以便获得所述已知序列在所述目标物种中的对应序列。

根据本发明实施例的系统适于执行根据本发明实施例的基于已知序列确定目标物种中对应序列的方法，根据本发明实施例的系统可用于从全基因组大数据中，基于已知序列，高效、精准捕获目标物种中对应序列的长片段目标区域(目前可实现～200Kb)，并从结果中可获得大量的信息供给进行下一步的判断，为大规模动植物进化和遗传学研究提供了强有力的技术支持。

具体地，所述测序读段是来自双向测序数据，优选地，来自BGI SEQ 500的双向测序数据。

具体地，所述测序读段Kmer序列集合获得装置仅针对双向测序读段中的正向测序读段构建所述测序读段Kmer序列集合。

具体地，分别采用相同的序列长度K值和间距D值，确定所述已知序列的全部Kmer序列和所述测序读段的Kmer序列，所述K值为27-39的整数，优选地，K值为31，所述D值为1。

具体地，所述测序读段Kmer序列集合获得装置200中包括：

优化处理单元，所述优化处理单元用于基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列之前，预先对所述测序读段的至少一部分进行优化处理。

可选地，所述优化处理单元适于执行以下操作：

删除所述测序读段起端和末端至少之一的1～10个碱基。

具体地，所述目标物种中的对应序列获得装置400适于执行以下操作：当所述延伸处理的延伸产物长度达到预定长度的1.5～2.5倍时，停止所述延伸处理，优选地，当所述延伸处理的延伸产物长度达到预定长度的2倍时，停止所述延伸处理。

具体地，所述确定延伸起始序列装置300还适于确定多个延伸起始序列，相应地，所述目标物种中的对应序列获得装置400还适于执行以下操作：

分别针对所述多个延伸起始序列分别进行所述延伸处理，以便获得多个延伸产物，并通过下列步骤从所述多个延伸产物中选择最终序列作为所述已知序列在所述目标物种中的对应序列：

分别确定所述多个延伸产物的每一个中所包含的种子序列Kmer序列数目；

选择包含所述种子序列Kmer序列数目最高的所述延伸产物作为所述对应序列。

具体地，参考图3，所述系统进一步包括：

候选测序读段确定装置500，所述候选测序读段确定装置500用于确定进行所述延伸处理所采用的测序读段Kmer序列所对应的候选测序读段；

修正装置600，所述修正装置600用于基于所述候选测序读段对所述对应序列进行修正处理。

可选地，所述修正处理包括以下操作：

基于所述候选测序读段与所述对应序列的比对，针对所述对应序列的至少一个位点，确定所述至少一个位点的优势碱基，并利用所述优势碱基对所述对应序列进行修正。

可选地，所述修正处理包括以下操作：

确定双向测序读段中成对测序读段在所述对应序列上的间距，如果所述间距与预定的插入片段长度差异超过10％，则判定所述对应序列为错误序列。

可选地，对于不存在所述优势碱基且测序深度差异低于2倍，则将所述位点标记为SNP。

可选地，对于不存在所述优势碱基且测序深度差异不低于2倍，则将所述位点标记为CNV。

下面将结合本发明的部分具体实施例，详细描述本发明的实施例，所述实施例的示例在附图中示出。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

实施例1

本实施例中，发明人以水稻为测试对象，选取一个小片段的WGS文库(包含～5％的叶绿体、～1％的线粒体数据)，插入片段～200bp，测序数据为～120Gb(1条lane)。为比较本发明的性能，分别尝试捕获叶绿体、线粒体、核基因组上的一个基因(随机选取)，并将其与华大自主研发组装软件(SOAPdenovo2)进行比较，结果如表1。

表1：SeedExt和SOAPdenovo2捕获目标区域序列结果比较

注：对于SOAPdenovo2(https://sourceforge.net/projects/soapdenovo2/)，采用的局部组装的方式，即先比对，筛选目标区域可能的Reads，再进行组装，相应的Duration包括比对时间。全局比对会消耗更多的时间，内存，存储，在这里不进行测试。在这里，本发明所使用的Seed均为拟南芥相应的基因。如果基因存在多外显子，那么，结果应该是包括内含子序列的完整基因序列，比较其准确性时，需要同时考虑内含子和外显子。

通过测试，可以看到，无论是叶绿体，线粒体，还是核基因组，捕获其相应的基因的完整序列时，本发明在运行时间，内存使用，存储等方面都低于常规方法。而对于组装的结果，本发明经以一条contig的形式出现，很大程度上得益于线性延伸算法。另外，通过与水稻自身的序列进行比对，发现本发明在覆盖度，准确性远高于常规方法。因此，说明本发明在全基因组的大数据中能够高效、精准地捕获目标区域序列。

总之，相比于常规方法，排除硬件条件，本发明具有以下优势：高覆盖度，高准确性；可实现大样品量处理。

实施例2

本实施例中，发明人以寄生植物Santalumalbum(SRR5150442)为研究对象，探索该物种叶绿体基因组中NAD(P)H dehydrogenase家族基因情况，包括ndhA,ndhB,ndhC,ndhD,ndhE,ndhF,ndhG,ndhH,ndhI,ndhJ和ndhK的11个基因。

具体探索过程如下所述：

1)确定物种目标区域是否残缺。

下载该物种相应的全基因组数据，使用本发明的方法，以拟南芥相应的基因为Seed，每个基因独一搜索，发现延伸结果均为空，反映为该寄生植物完全丢失该家族的所有基因，也没有出现该家族基因转移到线粒体或核基因组的情况。

其中，数据路径如下：

ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR515/SRR5150442/SRR5150442.sra

运行命令如下：

SeedExt-Seed ndh_seed.fa-K 31-D 5-Emax 5000-Rmax 20G-Read SRR5150442_1.fq.gz,SRR5150442_2.fq.gz

为证明该情况是否属实，确定Santalumalbum最近缘，且NCBI已发布相应的叶绿体基因组Taxilluschinensis(KY996492)，通过比较，发现该物种叶绿体也完全丢失这个家族的基因。为进一步证明，发明人将全基因组原始数据比对到最近缘物种葡萄叶绿体基因组(NC_023790)，该物种存在这个家族的所有基因，但发现完全没有相应的数据支持，确定其完全丢失。

2)确定目标区域是否存在变异(SNPs、INDEL等)。

以一株栽培稻和一株野生稻为研究对象，由华大建库，BGISEQ-500测序完成，并选取相应的小片段文库数据(～200bp，～120G，1条lane)，检测两者在OLE1基因(Chr4,27362086-27362954，-)是否存在变异。

使用本发明的方法，以拟南芥相应基因作为Seed，每个个体独一搜索，发现野生稻在基因的起始位点前后位置丢失66个核苷酸，对其结果进行比较，并根据输出的原始下机数据进行比对，存在较好的Paired End支持，确定其真实情况(参考图4)。其中，图4中的S1为栽培稻，S2为野生稻。两物种相应的序列都由本发明延伸获得，横坐标为位置(起点为Chr4，27361586)，中间灰色区域为两者比对结果，剪头部分为基因，方向为3’->5’，反映为负链上的基因，即基因起点在右端，从图5上可以看出S1是存在整个基因序列，而S2无法比对到基因的起始端，则证明S2缺失相应的序列，即起始端部分，即S2在基因起始位发生缺失。从数据的覆盖来看，Paired End很好地支持延伸的结果，且较为均匀，反映其准确性。

3)确定目标区域是否存在多拷贝或重复片段。

以一株野生稻为研究对象，由华大建库，采用BGISEQ500RS-1Ad-WGS试剂盒建库，主要步骤分为打断，加接头，PCR，环化。BGISEQ-500测序完成，并选取相应的小片段文库数据(～200bp，～120G，1条lane)，检测某个区域片段是否存在变异的情况(Chr10，2990853-2995534，+)是否存在变异。

使用本发明的方法，以拟南芥相应区域序列作为Seed进行搜索，对其结果进行比较，发现野生稻在该区域并没有发生结构性变化，并根据输出的原始下机数据进行比对，统计其相应的测序深度，发现在2993853-2994853区域，其测序尝试为其它的～2倍，可能为一个CNV(参考图5)。其中，图5中S3为野生稻。该物种相应的序列由本发明延伸获得，纵坐标为覆盖深度，横坐标为位置(起点为Chr10，2990853)。从数据的覆盖来看，整体深度较为均匀，但在2993853-2994853位置上，存在较高的深度(～2倍)，反映为基因组上的其它区域存在与该区域高度相似的一段序列。

4)确定目标基因是否存在多种可变剪切及其表达量(针对全转录组数据)。

以一株野生稻全转录组为研究对象，由华大建库，BGISEQ-500测序完成，并选取相应的小片段文库数据(～300bp，～6G)，检测其ELMO基因是否存在多种变剪切(Chr11，1360742-1364120，+)。

使用本发明的方法，以拟南芥相应基因作为Seed进行搜索，延伸结果存在2条序列，对其结果进行比较，发现中间一段序列存在差异，并根据输出的原始下机数据进行比对(hisat)，并统计每个位置的测序深度，确定其存在两种可变剪切，但是由于两者重叠部分较多，无法确定两者各自的表达量(参考图6)。其中，图6中标注的S4为野生稻。该物种相应的序列由本发明延伸获得，纵坐标为覆盖深度，横坐标为位置(起点为Chr11，1359742)。由于该数据为全转录组，本发明得到的延伸序列反映基因编码区，将其比对到基因组，体现为基因结构，连线为内含子，方框为外显子，在基因的前后依然有数据支持，反映为5’-UTR和3’-UTR。从数据的覆盖来看，整体深度较为均匀，且支持存在两种可变剪切的存在。

在本发明中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接或彼此可通讯；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (10)

1.一种基于已知序列确定目标物种中对应序列的方法，其特征在于，包括：

(1)基于所述已知序列，确定所述已知序列的全部Kmer序列，以便获得种子序列Kmer序列集合；

(2)获取来自于所述目标物种的测序读段，并基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列，以便获得测序读段Kmer序列集合；

(3)确定同时存在于所述测序读段Kmer序列集合和所述种子序列Kmer序列集合的Kmer序列作为延伸起始序列；

(4)基于重叠原则，利用所述测序读段Kmer序列集合，对所述延伸起始序列进行延伸处理，以便获得所述已知序列在所述目标物种中的对应序列。

2.根据权利要求1所述的方法，其特征在于，所述已知序列来自于所述目标物种的亲缘物种。

3.根据权利要求1所述的方法，其特征在于，步骤(1)和步骤(2)中，分别采用相同的序列长度K值和间距D值，确定所述已知序列的全部Kmer序列和所述测序读段的Kmer序列，所述K值为27-39的整数，优选地，K值为31，所述D值为1。

4.根据权利要求1所述的方法，其特征在于，所述测序读段是来自于双向测序系统；

优选地，所述测序读段是来自BGI SEQ 500的双向测序数据；

任选地，在步骤(2)中，仅针对双向测序读段中的正向测序读段构建所述测序读段Kmer序列集合。

5.根据权利要求1所述的方法，其特征在于，在步骤(2)中，在基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列之前，预先对所述测序读段的至少一部分进行优化处理；

优选地，所述优化处理包括删除所述测序读段起端和末端至少之一的1～10个碱基。

6.根据权利要求1所述的方法，其特征在于，在步骤(4)中，当所述延伸处理的延伸产物长度达到预定长度的1.5～2.5倍时，停止所述延伸处理，优选地，当所述延伸处理的延伸产物长度达到预定长度的2倍时，停止所述延伸处理。

7.根据权利要求1所述的方法，其特征在于，在步骤(3)中，确定多个延伸起始序列，在步骤(4)中分别针对所述多个延伸起始序列分别进行所述延伸处理，以便获得多个延伸产物，并通过下列步骤从所述多个延伸产物中选择最终序列作为所述已知序列在所述目标物种中的对应序列：

分别确定所述多个延伸产物的每一个中所包含的种子序列Kmer序列数目；

选择包含所述种子序列Kmer序列数目最高的所述延伸产物作为所述对应序列。

8.根据权利要求4所述的方法，其特征在于，进一步包括：

(5)确定步骤(4)中进行所述延伸处理所采用的测序读段Kmer序列所对应的候选测序读段；

(6)基于所述候选测序读段对所述对应序列进行修正处理。

9.根据权利要求8所述的方法，其特征在于，所述修正处理包括：

基于所述候选测序读段与所述对应序列的比对，针对所述对应序列的至少一个位点，确定所述至少一个位点的优势碱基，并利用所述优势碱基对所述对应序列进行修正；

优选地，所述修正处理包括：

确定双向测序读段中成对测序读段在所述对应序列上的间距，如果所述间距与预定的插入片段长度差异超过10％，则判定所述对应序列为错误序列；

任选地，对于不存在所述优势碱基且测序深度差异低于2倍，则将所述位点标记为SNP；

任选地，对于不存在所述优势碱基且测序深度差异不低于2倍，则将所述位点标记为CNV。

10.一种基于已知序列确定目标物种中对应序列的系统，其特征在于，包括：

种子序列Kmer序列集合获得装置，所述种子序列Kmer序列集合获得装置用于基于所述已知序列，确定所述已知序列的全部Kmer序列，以便获得种子序列Kmer序列集合；

测序读段Kmer序列集合获得装置，所述测序读段Kmer序列集合获得装置用于获取来自于所述目标物种的测序读段，并基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列，以便获得测序读段Kmer序列集合；

确定延伸起始序列装置，所述确定延伸起始序列装置与所述种子序列Kmer序列集合获得装置和所述测序读段Kmer序列集合获得装置相连，用于确定同时存在于所述测序读段Kmer序列集合和所述种子序列Kmer序列集合的Kmer序列作为延伸起始序列；

目标物种中的对应序列获得装置，所述目标物种中的对应序列获得装置与所述确定延伸起始序列装置相连，用于基于重叠原则，利用所述测序读段Kmer序列集合，对所述延伸起始序列进行延伸处理，以便获得所述已知序列在所述目标物种中的对应序列；

任选地，所述已知序列来自于所述目标物种的亲缘物种；

任选地，分别采用相同的序列长度K值和间距D值，确定所述已知序列的全部Kmer序列和所述测序读段的Kmer序列，所述K值为27-39的整数，优选地，K值为31，所述D值为1；

优选地，所述测序读段是来自双向测序系统；

任选地，所述测序读段是来自BGI SEQ 500的双向测序数据；

任选地，所述测序读段Kmer序列集合获得装置仅针对双向测序读段中的正向测序读段构建所述测序读段Kmer序列集合；

任选地，所述测序读段Kmer序列集合获得装置中包括：

优化处理单元，所述优化处理单元用于基于所述测序读段的至少一部分，确定所述测序读段的Kmer序列之前，预先对所述测序读段的至少一部分进行优化处理；

任选地，所述优化处理单元适于执行以下操作：

删除所述测序读段起端和末端至少之一的1～10个碱基；

任选地，所述目标物种中的对应序列获得装置适于执行以下操作：当所述延伸处理的延伸产物长度达到预定长度的1.5～2.5倍时，停止所述延伸处理，优选地，当所述延伸处理的延伸产物长度达到预定长度的2倍时，停止所述延伸处理；

任选地，所述确定延伸起始序列装置还适于确定多个延伸起始序列，所述目标物种中的对应序列获得装置还适于执行以下操作：

分别针对所述多个延伸起始序列分别进行所述延伸处理，以便获得多个延伸产物，并通过下列步骤从所述多个延伸产物中选择最终序列作为所述已知序列在所述目标物种中的对应序列：

分别确定所述多个延伸产物的每一个中所包含的种子序列Kmer序列数目；

选择包含所述种子序列Kmer序列数目最高的所述延伸产物作为所述对应序列；

任选地，所述系统进一步包括：

候选测序读段确定装置，所述候选测序读段确定装置用于确定进行所述延伸处理所采用的测序读段Kmer序列所对应的候选测序读段；

修正装置，所述修正装置用于基于所述候选测序读段对所述对应序列进行修正处理；

任选地，所述修正处理包括以下操作：

基于所述候选测序读段与所述对应序列的比对，针对所述对应序列的至少一个位点，确定所述至少一个位点的优势碱基，并利用所述优势碱基对所述对应序列进行修正；

任选地，所述修正处理包括以下操作：

确定双向测序读段中成对测序读段在所述对应序列上的间距，如果所述间距与预定的插入片段长度差异超过10％，则判定所述对应序列为错误序列；

任选地，对于不存在所述优势碱基且测序深度差异低于2倍，则将所述位点标记为SNP；

任选地，对于不存在所述优势碱基且测序深度差异不低于2倍，则将所述位点标记为CNV。

Images