CN111128305B - 对具有已知序列的生物序列进行分析的方法和系统 - Google Patents

Abstract

本发明提出了一种对具有已知序列的生物序列进行分析的方法。该方法包括：(a)基于所述生物序列，确定所述生物序列的全部Kmer序列，所述生物序列是通过对大片段氨基酸序列或核苷酸序列进行窗口划分获得的；(b)确定所述全部Kmer序列的每一个的频数，并基于所述频数，确定至少一个高频Kmer和至少一个低频Kmer；(c)基于所述至少一个高频Kmer，确定初步重复候选区域，其中，基于所述低频Kmer与相邻高频Kmer之间的距离，确定是否将所述低频Kmer整合至所述初步重复候选区域中。

Description

对具有已知序列的生物序列进行分析的方法和系统

技术领域

本发明涉及生物信息领域，具体地，本发明涉及对具有已知序列的生物序列进行分析的方法和系统。

背景技术

研究表明，在高等生物的基因组中非编码区都占到基因组序列的绝大部分，如人的基因组～3Gb，但非编码区占有高达～97％的比例。而绝大部分非编码序列以高度重复序列的形式存在，如卫星、小卫星、微卫星、长散布元件、短散布元件等，各种重复序列的类型与它在染色体上的分布密切相关。

以前，人们认为重复序列不过是一些冗余，或“无用”DNA。然而，大量的实验及研究表明：重复序列不是垃圾，而是影响着生命的进化、遗传、变异；同时它对基因表达、转录调控、染色体的构建以及生理代谢都起着不可或缺的作用。例如，一些三核苷酸重复序列拷贝数的异常增加会导致某些人类遗传病的产生，如脆性X染色体综合症。另外，随着研究的深入，发现基因重复的蛋白功能域，常被应用于结构锚定模式，与生物聚合物平稳地相互作用。例如，具有Tetratricopeptide repeats(TPRs)，ankyrin(ANK)repeats的蛋白，重复单元分别为34和33个氨基酸，它们都形成一个helix-turn-helix结构。在原核到真核的整个进化过程中，这类家族都非常保守。这些重复的功能域已经被报道与其它的蛋白和RNA相互作用，在细胞周期调控，转录调节，转化抑制，蛋白易位上扮演着重要的角色。因此，识别非编码区核酸，或编码区蛋白序列的重复序列是分析其功能的基础。

其中，串联重复序列(tandemrepeat)是指以一定的碱基数作为重复单元，首尾相连排列在一起形成聚集区的重复序列。在此基础上，提出了周期性重复序列的概念，类似于串联重复序列，但允许以下特殊情况：不同的重复单元间存在差异(在长期的进化过程上，会出现少量的错配或Gap)；相邻的重复单元间存在其它序列(即被某些序列随机隔开)；部分蛋白序列上也存在重复单元(尤其是功能域区域)。串联重复序列是这类重复序列中的一种情况，表现为在某个区域内，出现相对集中的重复单元“簇”。

串联重复序列的识别问题，根据所采用方法的不同，现有的重复序列发现方法可以分为2类，这2种方法都能识别出基因组序列中串联重复序列出现的位置：

(1)基于字符串精确匹配的方法：可理解为判断是否存在重复单元，且其存在的形式是否为串联，如TRF。

(2)基于数字信息处理的方法：采用二进制方法来表示各个碱基，并分别求出各碱基的频谱，最后将4个碱基的频谱相加得到序列的总频谱。观察频谱图可得到序列中串联重复序列拷贝出现的频率，如SRF。

然而，序列在长期的进化过程中，存在某些位置上发生例如插入、删除、替换等突变，因此，对于串联重复序列，其重复单元可能存在差异。另外，串联重复序列识别中要处理的数据量一般都比较大，往往是整个基因组，因此，计算量是识别方法中应该考虑的重要问题。而现有的方法，

(1)基于字符串精确匹配的方法，一方面不能保证发现序列中所有可能的串联重复序列，另一方面，这种方法的计算复杂度会随着序列中串联重复序列拷贝的长度呈现指数形式增长。

(2)基于数字信号处理的方法，一方面对核酸序列采用二进制表示法，需要对每条序列做4次离散傅里叶变换才能求出核酸序列的频谱图，计算量大。另外，需要针对每个串联重复拷贝频率分别求其加窗傅里叶变换，才能得到核酸序列中所有串联重复序列出现的位置，识别灵敏度低。

(3)目前只针对核酸序列，而不能实现对蛋白序列的预测，而研究发现，蛋白序列中重复的功能域也发挥着重要作用。

(4)目前只实现从头预测的方法，而不能针对某些目标序列进行预测。由于某种需要，有可能只需要判断在特定区域内是否存在预期的具有某种特征序列的重复单元，即基于已提供的特征序列进行预测。

因此，新的重复序列发现方法仍需要进一步开发和改进。

发明内容

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

对于周期性重复序列的预测，发明人认为，关键在于是否可以在某个区域找到相似的重复单元出现相对集中的情况，是一“簇”的概念(对于疏远重复单元的情况，不属于本申请识别的范围)；并且周期性重复序列既然是重复序列的一种，同样体现为某种序列多次重现，即重复单元，而当重复单元过短时，由于现有的比对软件基本都涉及到种子定位和延伸的问题，因此，使用自身比对的方式无法实现这种重复单元过短的周期性重复序列的识别。

针对这些问题，本发明提出了一种基于短Kmer识别重复单元相对集中的情况，确定存在周期性重复的单元序列和完整区域。

在本发明的第一方面，本发明提出了一种对具有已知序列的生物序列进行分析的方法。根据本发明的实施例，所述方法包括：(a)基于所述生物序列，确定所述生物序列的全部Kmer序列，所述生物序列是通过对大片段氨基酸序列或核苷酸序列进行窗口划分获得的；(b)确定所述全部Kmer序列的每一个的频数，并基于所述频数，确定至少一个高频Kmer和至少一个低频Kmer；(c)基于所述至少一个高频Kmer，确定初步重复候选区域，其中，基于所述低频Kmer与相邻高频Kmer之间的距离，确定是否将所述低频Kmer整合至所述初步重复候选区域中。根据本发明实施例的分析方法将序列信息转化为Kmer以及相应的频数，无论是核酸，还是蛋白序列，都不影响Kmer的生成和计算，因此，根据本发明实施例的方法既可实现核酸序列的预测，也可以实现蛋白序列的预测。根据本发明实施例的方法实现了高效、精准地预测全基因组，或基因蛋白序列中的周期性重复序列的单元序列和完整区域，进而可根据预测的周期性重复序列的位置特征，确定周期性重复序列区域内部或附近的调控元件如启动子、增加子、终止子，rRNA基因和组蛋白基因、编码基因等，分析周期性重复序列在参与顺式调控元件、基因表达、表观遗传修饰等重要过程中所发挥的作用。因此，根据本发明实施例的方法为大规模动植物进化和遗传学研究提供了强有力的技术支持。

根据本发明的实施例，上述方法还可以进一步包括如下附加技术特征至少之一：

根据本发明的实施例，所述生物序列为氨基酸序列或者核酸序列。如前所述，根据本发明实施例的分析方法将序列信息转化为Kmer以及相应的频数，无论是核酸，还是蛋白序列，都不影响Kmer的生成和计算，因此，根据本发明实施例的方法既可实现核酸序列的预测，也可以实现蛋白序列的预测。

根据本发明的实施例，所述生物序列为核酸序列，所述Kmer序列的长度为10个核苷酸。

根据本发明的实施例，所述生物序列为氨基酸序列，所述Kmer序列的长度为3个氨基酸。

发明人发现，Kmer序列的长度过长，会导致短串联重复无法被识别，而当Kmer序列的长度过短时，会导致在基因组任意位置匹配到的概率非常高，增加判断串联重复单元序列的难度。Kmer序列的长度为10个核苷酸或3个氨基酸，可进一步提高本发明方法预测的精准性。

根据本发明的实施例，对于氨基酸序列，所述生物序列的长度为200～500个氨基酸，优选300个氨基酸，对于核酸序列，所述生物序列的长度为800～1500bp，优选1000bp。

根据本发明的实施例，在步骤(b)中，将所述全部Kmer序列每一个的频数与预定的频数阈值进行比较，以便确定所述高频Kmer和所述低频Kmer，其中，所述预定阈值是所述至少一个Kmer序列中最高频数的至少0.3倍。

根据本发明的实施例，在步骤(c)中，当所述低频Kmer与相邻高频Kmer之间的距离小于预定阈值时，将所述低频Kmer的序列整合到所述初步重复候选区域中。

根据本发明的实施例，所述方法进一步包括通过下列步骤，基于所述初步重复候选区域确定重复单元：(c-1)通过对所述初步重复候选区域进行延伸处理，得到经过延伸的初步候选区域；(c-2)沿着所述经过延伸的初步候选区域的预定顺序，确定首个高频Kmer，并确定所述首个高频Kmer在所述经过延伸的初步候选区域中的位置；(c-3)基于下游Kmer的频数，沿着所述首个高频Kmer的位置向下游进行延伸，直到遇到频数实质变化Kmer，停止所述延伸，以便获得至少一个重复单元。进而可有效防止漏掉部分相隔较远的散置的周期性重复序列单元，进一步提高获得完整的周期性重复序列区域的概率。

需要说明的是，频数实质变化Kmer是指沿着延伸的方向上，特定Kmer的出现使得该特定Kmer与上游临近Kmer的频数之间存在显著差异，例如，针对两个临近的Kmer，在包括但不限于下列情形下，下游Kmer可以被认定为频数实质变化Kmer：

(1)上游Kmer为高频序列，下游Kmer为低频序列，并且该低频序列距离其下游的高频序列的距离超过了预定阈值，则该下游Kmer可以被认定为频数实质变化Kmer；

(2)上游Kmer和下游Kmer均为高频序列，但其二者之间的频数差异超过了20％，例如超过了30％，40％，50％，则认定下游Kmer为频数实质变化Kmer；

(3)上游Kmer为高频序列，下游Kmer为低频序列，并且该低频序列距离其下游的高频序列的距离没有超过预定阈值，则如果该低频序列的下游高频序列与该上游Kmer之间的频数差异超过了20％，例如超过了30％，40％，50％，则认定下游Kmer为频数实质变化Kmer。

根据本发明的实施例，所述方法进一步包括对所述重复单元进行评估，确定所述重复单元的PR分数，所述PR分数的计算方法可参考文献Mori H,et al.Nucleic AcidsResearch,doi.org/10.1093，进而根据所获得的PR分数，记录满足条件的所有重复单元，用于下一步的分析。

根据本发明的实施例，所述方法进一步包括：将所述至少一个重复单元进行比对，以便确定所述至少一个重复单元的一致性序列。对于周期性重复序列，是通过相应的单元序列的重复出现反映其周期性的，但是由于存在变异，并不是每个单元序列都是相同的。因此，为了反映该周期性重复序列的主要特征，需要通过比对，确定所有单元序列最为保守的部分，并作为周期性重复序列的特征序列。根据本发明的具体实施例，基于(c-3)的结果，提取满足条件的单元序列，借助软件(-PRCTool)，如MAFF，MUSCLE等，进行多序列比对，并将对比对结果进行如下处理：对于每个位置进行统计，包括每个碱基或氨基酸及gap的比例，当最主要的成分达不到参数-对应的每个位置相同碱基或氨基酸的最小比例(PRCratio)，则进行移除，最终获得周期性重复序列单元的一致性序列。

在本发明的第二方面，本发明提出了一种对具有已知序列的生物序列进行分析的系统。根据本发明的实施例，所述系统包括：Kmer序列确定装置，所述Kmer序列确定装置用于基于所述生物序列，确定所述生物序列的全部Kmer序列，所述生物序列是通过对大片段氨基酸序列或核苷酸序列进行窗口划分获得的；高频Kmer和低频Kmer确定装置，所述高频Kmer和低频Kmer确定装置与所述Kmer序列确定装置相连，用于确定所述全部Kmer序列的每一个的频数，并基于所述频数，确定至少一个高频Kmer和至少一个低频Kmer；初步重复候选区域确定装置，所述初步重复候选区域确定装置与所述高频Kmer和低频Kmer确定装置相连，用于基于所述至少一个高频Kmer，确定初步重复候选区域，其中，基于所述低频Kmer与相邻高频Kmer之间的距离，确定是否将所述低频Kmer整合至所述初步重复候选区域中。根据本发明实施例的系统适于执行上述对具有已知序列的生物序列进行分析的方法，既可实现核酸序列的预测，也可以实现蛋白序列的预测。

根据本发明实施例的上述对具有已知序列的生物序列进行分析的系统所具有的附加技术特征以及技术效果与上述根据本发明实施例的方法类似，在此不再赘述。

附图说明

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

图1是根据本发明实施例的对具有已知序列的生物序列进行分析的方法的流程图；

图2是根据本发明实施例的对具有已知序列的生物序列进行分析的系统的结构示意图；

图3是根据本发明再一实施例的对具有已知序列的生物序列进行分析的系统的结构示意图；

图4是根据本发明再一实施例的对具有已知序列的生物序列进行分析的系统的结构示意图；

图5是根据本发明再一实施例的对具有已知序列的生物序列进行分析的系统的结构示意图；以及

图6是根据本发明实施例的X物种的目标区域预测结果图，其中，横坐标为位置，Chr1:500000～600000，中间黑色箭头为基因，箭头左右方向表示基因位于染色体的正负链，柱状图宽度表示周期性重复序列单元一致性序列的长度，高度表示单元相应的拷贝数。

具体实施方式

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

一般定义

Kmer，长度为K的短序列。如长度为L的序列，从起始端开始，每次移动一个碱基，则可得到(L-K+1)个Kmer。

Kmer Depth，Kmer频数。如长度为L的序列，那么，该序列将产生(L-K+1)个Kmer，进而可获得每种Kmer出现的频数。

Repeat，重复序列。某特定的核酸序列在基因组中重复出现，而该发明中，除了核酸序列之外，还包括蛋白序列。

Periodic Repeat Element，周期性重复序列单元序列。在某个区域中，按位置变动，呈现重复性行为的序列。这些单元序列间存在差异。

ElementConsensus Sequence，周期性重复序列单元的一致性序列。在该发明中，它由多条周期性重复序列单元通过多序列比对，去除非保守区域后的结果。

Periodic Repeat Region，周期性重复序列区域。对于某个区域，大部分由周期性重复序列单元组成，可为串联或散置形式。

Annotation，注释。在该发明中，主要判断是否有基因包含于某个周期性重复序列区域，或者某个周期性重复序列区域包含于某个编码区中。

InFile，本发明的输入文件。本发明可适用于处理核酸或者蛋白序列，同时，输入文件可为FASTA，JSON，或GenBank格式。

Region，目标区域。本发明可针对某个目标区域进行分析，只需要提供染色体，起始和终点。

Window，滑动窗口。相邻的窗口没有重叠部分，本发明将统计该窗口内所有的Kmer的频数。默认核酸为1000，氨基酸为300。

K，Kmer的长度。在本发明中，当取较大K时，会导致短串联重复无法被识别，如K大于短串联重复序列总长度，故本发明建议使用较小的K值。默认核酸为10，氨基酸为3。

HKmer，高频Kmer的最小值。本发明通过该参数，将滑动窗口中所有Kmer分为高频和低频。默认为0.3。这个参数不是设置一个整数，而是一个小数，反映为大窗口中最高Kmer频率的HKmer倍，默认为0.3倍。主要考虑到每个区域的周期性重复序列出现的频率不相同，但始终体现出比邻近区域来得高。为了防止漏掉，设置一个较小的数，时间上可能会延长。另外，可能会检测到较多的区域，但本发明后续会基于PRElen过滤掉部分异常单元序列，且最后还会将单元一致性序列比对回参考序列进行单元序列和边界修正。

HkR，高频Kmer区域的最小长度。本发明统计滑动窗口每个位置的频数，连续的高频区域表现为潜在的重复序列。默认核酸为20，氨基酸为6。

HkRGap，高频Kmer区域区域允许的Gap最大值。两个相邻的周期性重复序列单元间的非重复区域长度，应用于检测散置的重复序列。默认核酸为200，氨基酸为50。

HkRExt,高频Kmer区域的延伸长度最大值。确定高频Kmer区域之后，即重复序列区域，对上下游进行延伸，以检测完整的周期性重复序列区域。默认核酸为1000，氨基酸为300。

PRElen，周期性重复序列单元的长度范围。由于变异，周期性重复序列单元的长度是一个波动的范围。在本发明中，自动识别周期性重复序列单元长度的平均值，在平均值某范围内波动的单元才被用于判断是否为周期性重复序列。默认为0.8-1.2。

PREgap，周期性重复序列单元最大的Gap长度。由于存在变异，允许周期性重复序列单元有一定的错配或Gap。默认核酸为3，氨基酸为1。

PRscore，周期性重复序列的最小分值。在重复序列区域的基础上，对周期性重复序列进行定义，体现为单元序列出现的规律是否具有周期性。默认核酸为0.5，氨基酸为0.3。

PRCratio，定义周期性重复序列单元一致性序列时，其对应的每个位置相同碱基或氨基酸的最小比例。默认核酸为0.8，氨基酸为0.8。

PRCmerge，对相邻窗口鉴定的周期性重复序列区域的最小重叠区域进行合并。只有满足相邻窗口的周期性重复序列单元一致性序列相同时才进行操作。

PRCsite，周期性重复序列单元一致性序列或文件输入。本发明允许只检测所提供的某些特定类型的周期性重复序列，包括核酸和蛋白序列。

PRCNum，针对已提供单元一致序列判断目标区域存在周期性重复序列的最小拷贝数。

PRCTool，获得周期性重复序列单元一致性序列时所使用的工具。本发明需借助其它的软件进行多序列比对，以观察每个位置出现的碱基或氨基酸。

PRRTool，获得完整的周期性重复序列区域时所使用的工具。本发明需借助其它的软件将一致性序列重新比对到重复区域，以确定完整的周期性重复序列区域。

对具有已知序列的生物序列进行分析的方法

一方面，本发明提出了一种对具有已知序列的生物序列进行分析的方法。根据本发明的实施例，所述方法包括：(a)基于所述生物序列，确定所述生物序列的全部Kmer序列,所述生物序列是通过对大片段氨基酸序列或核苷酸序列进行窗口划分获得的；(b)确定所述全部Kmer序列的每一个的频数，并基于所述频数，确定至少一个高频Kmer和至少一个低频Kmer；(c)基于所述至少一个高频Kmer，确定初步重复候选区域，其中，基于所述低频Kmer与相邻高频Kmer之间的距离，确定是否将所述低频Kmer整合至所述初步重复候选区域中。如前所述，根据本发明实施例的方法既实现了核酸序列，也实现了蛋白序列中周期性重复序列的高效、精准地预测，为进一步研究周期性重复序列的结构和功能奠定了基础，为进行大规模动植物进化和遗传学研究提供了强有力地技术支持。

为了便于理解，申请人详细介绍了根据本发明实施例的对具有已知序列的生物序列进行分析的方法的具体执行过程，具体参见图1所示流程：

1)自动识别输入文件格式及部分参数

根据输入的序列，自动识别输入文件格式，确定其序列ID及相应的序列。根据参数-Region[Chr1,Start1,End1；Chr2,Start2,End2；…]提供的坐标位置信息或文件，自动截取相应的目标区域序列；否则，所有序列均为分析对象。序列内容为核酸，或蛋白序列。

2)检测重复序列区域

输入的序列有可能非常大，而根据本发明实施例的方法所检测的周期性重复序列只占较少的一部分。因此，在检测周期性重复序列之前，发明人先探索了潜在重复区域。

本发明通过参数-Window设置滑动窗口，且相邻窗口没有重叠区域，那么，每个窗口为本发明处理鉴定的对象，每个窗口的处理均为独立的。对于每个窗口，设置Kmer长度为K，计算窗口内所有Kmer深度。通过设置高频Kmer的域值-HKmer，从而获得连续高频Kmer区域。由于杂合或变异位点的存在，有可能相邻高频区夹杂着较短的低频，当低频区域长度满足参数-HkRGap时，被允许当作高频区。当高频区域的长度满足-HkR时，该区域被定义为重复序列区域，用于下一步的分析。否则，程序终止。

3)评估周期性重复序列

基于2)所检测到的重复序列区域，根据参数-HkRExt，对相应的区域进行上下游延伸。延伸的主要目的是防止漏掉部分相隔较远的散置的周期性重复序列单元，从而导致无法获得完整的周期性重复序列区域。基于重复序列区域，从起始端开始，每次移动一个碱基或氨基酸，当第一次发现高频Kmer时，记录重复序列区域上所有该高频Kmer位置。对于记录的每个位置的第二位开始往后迭代，直到遇到低频区(其长度不得大于参数-PREgap，否则，迭代停止)，或者相邻的高频区，迭代停止。那么，从迭代起点到终点被定义重复序列单元。根据重复序列单元的位置信息，可以确定任意单元间长度Elen，根据参数-PRElen，根据满足条件的单元，可获得单元最高Kmer深度的Kmax。那么，该重复序列为周期性重复序列可能性为：(Elen1+Elen2+…+ElenN)/(1+2+…+Kmax/2)。当分值满足PRscore时，记录满足条件的所有单元，用于下一步的分析。否则，程序终止。

4)确定周期性重复序列一致性序列

对于周期性重复序列，是通过相应的单元序列的重复出现反映其周期性，但是由于存在变异，并不是每个单元序列都是相同的。因此，为了反映该周期性重复序列的主要特征，需要确定所有单元序列最为保守的部分，并作为周期性重复序列的特征序列。基于3)的结果，提取满足条的单元序列，借助其它软件(-PRCTool)，如MAFF，MUSCLE等，进行多序列比，对比对结果进行处理：对于每个位置进行统计，包括每个碱基或氨基酸及gap的比例，当最主要的成分达不到参数-PRCratio，则进行移除，最终获得周期性重复序列单元的一致性序列。

5)基于已提供的单元序列判断周期性重复序列区域

在某些情况下，有可能只关注具有某些特征的周期性重复序列，而无需对所有的可能的具有某种特征的周期性重复序列进行判断，即前面的从头预测的方法。因此，本发明设计中，除了从头预测之外，还可以通过设置参数-PRCsite，实现基于已知序列进行判断。所提供的序列，本发明理解为周期性重复序列的单元一致性序列，即等同于4)的结果。本发明以提供的序列长度作为Kmer的长度，从重复序列起始端开始，遍历每个位置，允许一定的错配和Gap，当其拷贝数达到参数-PRCNum时，可定义相应的区域存在周期性重复序列。

6)确定周期性重复序列完整区域

对于周期性重复序列，其单元序列可为串联，也可为散置。对于远离集中区的单元，或者单元长度小于Kmer长度，以上的处理无法获得周期性重复区域的准确信息。在此，基于4)获得的单元一致性序列，根据参数-PRCmerge对相邻窗口的周期性重复序列单元一致性序列进行整合，并借助于比对软件(-PRRTool)，如BLAT等，将其比对到2)定义的重复序列区域(延伸后)，判断比对上的区域是否与3)的结果存在差异，如漏掉，或存在坐标重叠但不完全相同(反映为可能存在更短的单元)，从而实现对单元序列进行纠正。对比对区域的始端到末端为周期性重复序列的完整区域。

7)周期性重复序列注释

基于5)的定义的周期性重复序列区域，以及序列本身的注释文件(基因，tRNA等各类元件)。判断周期性重复序列区域是否包含或包含于各类元件，相应的基因与周期性重复序列在进化上具有重要意义。如某个编码基因为周期性重复单元的一部分，那么，反映该重复序列对基因的扩张和消失具有影响；或者周期性重复序列为某编码基因的一部分，那么，该重复序列对基因功能变异产生影响，如多功能化。

根据本发明实施例的序列分析方法使用一个滑动窗口(参数-Window，如1Kb)，进行无重叠地扫描整个基因组，每个窗口是一个分析的对象。根据本发明实施例的序列分析方法具有如下优势：

(1)对于每个窗口，每次滑动一个碱基，遍历所有位置，生成所有Kmer(参数-K，如10)，并统计每个Kmer相应的频数，由于滑动窗口小，存储Kmer少，后续分析所占用的内存小；

(2)每个窗口的分析是独立的，且可以实现并行分析，因此，可以通过并行处理的方式大量减少计算时间(相对时间)，可实现在短时间内，对大数据进行预测；

(3)本发明将序列信息转化为Kmer以及相应的频数。本发明涉及到两种窗口，大窗口-Window，无重叠扫描过整个基因组；小窗口-K，遍历每个大窗口的所有位置。本发明中只需记录大窗口的坐标，当检测到大窗口存在周期性重复序列时，根据大窗口坐标可确定周期性重复序列的准确位置。无论是核酸，还是蛋白序列，都不影响Kmer的生成和计算，因此，也可以实现蛋白序列的预测；

(4)通过设置Kmer的频数(参数-HKmer)，当确定窗口序列存在重复序列时(对于任意的大窗口，如果高频Kmer区域满足最小长度HkR(核酸>200b；氨基酸>50bp)，判断该大窗口存在重复序列)，根据窗口坐标进行延伸(参数-HkRExt，如1Kb)，以防止漏掉边缘的重复单元，然后根据自定义的规则(对大窗口中的重复序列进行打分：(Elen1+Elen2+…+ElenN)/(1+2+…+Kmax/2)。Elen为各单元长度，Kmax为单元最高深度，当分值大于PRscore(核酸>0.5；氨基酸>0.3)，反映单元序列在该区域出现的频率高)判断重复序列是否为周期性重复序列；

(5)当相邻窗口同时被判断存在周期性重复序列时，判断相邻窗口的周期性重复序列单元的序列特征是否一致。相同时应进行整合，可防止数目或区域的错误计算；

(6)本发明将序列转化成Kmer，并依据Kmer频数进行后续的分析，那么，假设只关注具有某种特征序列的周期性重复序列呢？本发明将其当成Kmer，然后搜索滑动窗口是否存在相对集中的Kmer，因此，也可实现基于提供的特征序列进行预测。

对具有已知序列的生物序列进行分析的系统

另一方面，本发明提出了一种对具有已知序列的生物序列进行分析的系统。根据本发明的实施例，参考图2，所述系统包括：

Kmer序列确定装置100，所述Kmer序列确定装置100用于基于所述生物序列，确定所述生物序列的全部Kmer序列，具体地，所述Kmer序列确定装置100用于确定所述生物序列的全部Kmer序列,所述生物序列是通过对大片段氨基酸序列或核苷酸序列进行窗口划分获得的。对于氨基酸序列，所述窗口划分后的生物序列的长度为200～500个氨基酸，优选300个氨基酸，对于核酸序列，所述窗口划分后的生物序列的长度为800～1500bp，优选1000bp。当所述生物序列为核酸序列时，所述Kmer序列的长度设置为10个核苷酸，当所述生物序列为氨基酸序列，所述Kmer序列的长度设置为3个氨基酸，进而可进一步提高本发明系统预测的精准性；

高频Kmer和低频Kmer确定装置200，所述高频Kmer和低频Kmer确定装置200与所述Kmer序列确定装置100相连，用于确定所述全部Kmer序列的每一个的频数，并基于所述频数，确定至少一个高频Kmer和至少一个低频Kmer，其中，确定至少一个高频Kmer和至少一个低频Kmer可通过如下具体方式实现：将所述至少一个Kmer序列每一个的频数与预定的频数阈值进行比较，以便确定所述高频Kmer和所述低频Kmer，其中，所述预定阈值是所述至少一个Kmer序列中最高频数的至少0.3倍；

初步重复候选区域确定装置300，所述初步重复候选区域确定装置300与所述高频Kmer和低频Kmer确定装置200相连，用于基于所述至少一个高频Kmer，确定初步重复候选区域，其中，基于所述低频Kmer与相邻高频Kmer之间的距离，确定是否将所述低频Kmer整合至所述初步重复候选区域中，当所述低频Kmer与相邻高频Kmer之间的距离小于预定阈值时，将所述低频Kmer的序列整合到所述初步重复候选区域中。

根据本发明的再一具体实施例，参考图3，所述系统进一步包括重复单元确定装置400，所述重复单元确定装置400与所述初步重复候选区域确定装置300相连，用于基于所述初步重复候选区域确定重复单元，所述重复单元确定装置400包括：

延伸单元410，所述延伸单元410通过对所述初步重复候选区域进行延伸处理，得到经过延伸的初步候选区域；

首个高频Kmer确定单元420，所述首个高频Kmer确定单元420与所述延伸单元410相连，用于沿着所述经过延伸的初步候选区域的预定顺序，确定首个高频Kmer，并确定所述首个高频Kmer在所述经过延伸的初步候选区域中的位置；

重复单元确定单元430，所述重复单元确定单元430与所述首个高频Kmer确定单元420相连，用于基于下游Kmer的频数，沿着所述首个高频Kmer的位置向下游进行延伸，直到遇到频数实质变化Kmer，停止所述延伸，以便获得至少一个重复单元，所述实质低频Kmer的频数低于所述预定阈值并且与下游高频Kmer的距离超过预定阈值；

根据本发明的再一具体实施例，参考图4，所述重复单元确定装置400进一步包括评估单元440，所述评估单元440与所述重复单元确定单元430相连，用于对所述重复单元进行评估，确定所述重复单元的PR分数，所述PR分数的计算方法可参考文献Mori H,etal.Nucleic Acids Research,doi.org/10.1093，进而根据所获得的PR分数，记录满足条件的所有重复单元，用于下一步的分析；

根据本发明的再一具体实施例，参考图5，所述系统进一步包括一致性确定装置500，所述一致性确定装置500与所述重复单元确定装置400相连，用于将所述至少一个重复单元进行比对，以便确定所述至少一个重复单元的一致性序列。

根据本发明实施例的系统适于执行上述根据本发明实施例的对具有已知序列的生物序列进行分析的方法，其优势和效果与上述方法类似，在此不再赘述。

下面将结合实施例对本发明的方案进行解释。本领域技术人员将会理解，下面的实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件的，按照本领域内的文献所描述的技术或条件(例如参考J.萨姆布鲁克等著，黄培堂等译的《分子克隆实验指南》，第三版，科学出版社)或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品，例如可以采购自Illumina公司。

实施例1

发明人选取了四个植物为测试对象，分别预测其周期性重复序列。为比较本发明的性能，发明人将其与TRF软件进行比较，结果如表1。

表1：TRF和SearchPRE预测重复序列结果比较

注：对于以上的两个软件，均只使用1个CPU，其它均使用默认参数。

通过测试，可以看到，本发明对于四种植物均预测到更多周期性重复序列，包括数目和长度，且结果包含大部分TRF的预测结果，表明具有较高的准确性，而其它部分主要是因为本发明允许重复单元存在差异，以及相邻的重复单元间存在其它序列的特殊情况。而本发明在运行时间，内存使用方面都低于常规方法。因此，说明本发明在全基因组的大数据中能够高效、精准地预测周期性重复序列。

实施例2

在本实施例中，以X物种区域Chr1:500000～600000为研究对象，探索该区域是否存在周期性重复序列，具体步骤如下所述：

1)由于所分析的区域是局部，可通过参数-Region设置实现：-Region Chr1,500000,600000；

2)对于该区域相应的核酸序列，使用本发明进行预测，运行命令及参数，如下：

SearchPRE-InFile Ref.fasta-Region Chr1,500000,600000-Window 1000–K10–Hkmer 0.3–HkR 20–HkRGap 200-HkRExt 1000-PRElen 0.8-1.2–PREgap 3-PRscore0.5-PRCratio 0.8-Outdir All

从结果来看，如图6和表2，该区域只检测到1个周期性重复序列区域，位于编码区中。为了验证编码区的周期性重复序列的真实性，对该区域的所有基因蛋白序列进行预测。

3)根据目标区域的坐标(Chr1,500000,600000)，获取这个区域上所有基因的蛋白序列。

4)对提取出来的蛋白序列，使用本发明进行预测，运行命令及参数，如下：

SearchPRE-InFile Pep.fasta-Window 300–K 3–Hkmer 0.3–HkR 6–HkRGap 50-HkRExt 300-PRElen 0.8-1.2–PREgap 1-PRscore 0.3-PRCratio 0.8–Outdir Gene

从结果来看，如图6(上面表示基于核酸序列预测的周期性重复序列，下面表示基于基因的蛋白序列预测的周期性重复序列)和表2，对该区域的66个基因的蛋白序列，检测到5个基因存在周期性重复序列，但2)的结果并没有出现在该结果中。主要是因为在基因组上，基因存在多个外显子的情况，即跨越多个区域的序列，重复序列被内含子隔开，故无法检测到。相反，由于蛋白序列是连接所有外显子，并转化为氨基酸序列，与实际的基因组序列有差异，因此，出现蛋白序列有结果，而基因组序列没有的情况。

同时，使用TRF软件进行预测串联重复序列，并与SearchPRE结果进行比较。

5)由于TRF软件没有截取功能，因此根据目标区域的坐标(Chr1,500000,600000)，截取相应的序列并存储于Ref.Cut.fa文件。

6)对于提取的核酸序列，使用TRF软件进行预测，运行命令如下：trf Ref.Cut.fa2 7 7 80 10 50 2000 -d –h。

从结果来看，如表2，TRF预测到4个串联重复序列，其中有2个与SearchPRE结果重叠，另外2个拷贝数只为2.1和2.8。

表2：周期性重复序列预测结果比较

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

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

Claims (21)

1.一种对具有已知序列的生物序列进行分析的方法，其特征在于，包括：

（a）基于所述生物序列，确定所述生物序列的全部Kmer序列,所述生物序列是通过对大片段氨基酸序列或核苷酸序列进行窗口划分获得的；

（b）确定所述全部Kmer序列的每一个的频数，并基于所述频数，确定至少一个高频Kmer和至少一个低频Kmer；

（c）基于所述至少一个高频Kmer，确定初步重复候选区域，其中，基于所述低频Kmer与相邻高频Kmer之间的距离，确定是否将所述低频Kmer整合至所述初步重复候选区域中；当所述低频Kmer与相邻高频Kmer之间的距离小于预定阈值时，将所述低频Kmer的序列整合到所述初步重复候选区域中；

步骤（c）进一步包括：

（c-1）通过对所述初步重复候选区域进行延伸处理，得到经过延伸的初步候选区域；

（c-2）沿着所述经过延伸的初步候选区域的预定顺序，确定首个高频Kmer，并确定所述首个高频Kmer在所述经过延伸的初步候选区域中的位置；

（c-3）基于下游Kmer的频数，沿着所述首个高频Kmer的位置向下游进行延伸，直到遇到频数实质变化Kmer，停止所述延伸，以便获得至少一个重复单元。

2.根据权利要求1所述的方法，其特征在于，所述生物序列为氨基酸序列，所述Kmer序列的长度为3个氨基酸。

3.根据权利要求1所述的方法，其特征在于，所述生物序列为核酸序列，所述Kmer序列的长度为10个核苷酸。

4.根据权利要求2所述的方法，其特征在于，对于氨基酸序列，所述生物序列的长度为200~500个氨基酸。

5.根据权利要求2所述的方法，其特征在于，对于氨基酸序列，所述生物序列的长度为300个氨基酸。

6.根据权利要求3所述的方法，其特征在于，对于核酸序列，所述生物序列的长度为800~1500bp。

7.根据权利要求3所述的方法，其特征在于，对于核酸序列，所述生物序列的长度为1000bp。

8.根据权利要求1所述的方法，其特征在于，在步骤（b）中，将所述全部Kmer序列每一个的频数与预定的频数阈值进行比较，以便确定所述高频Kmer和所述低频Kmer，其中，所述预定阈值是所述至少一个Kmer序列中最高频数的至少0.3倍。

9.根据权利要求1所述的方法，其特征在于，进一步包括对所述重复单元进行评估。

10.根据权利要求1所述的方法，其特征在于，所述方法进一步包括：将所述至少一个重复单元进行比对，以便确定所述至少一个重复单元的一致性序列。

11.一种对具有已知序列的生物序列进行分析的系统，其特征在于，包括：

Kmer序列确定装置，所述Kmer序列确定装置用于基于所述生物序列，确定所述生物序列的全部Kmer序列，所述生物序列是通过对大片段氨基酸序列或核苷酸序列进行窗口划分获得的；

高频Kmer和低频Kmer确定装置，所述高频Kmer和低频Kmer确定装置与所述Kmer序列确定装置相连，用于确定所述全部Kmer序列的每一个的频数，并基于所述频数，确定至少一个高频Kmer和至少一个低频Kmer；

初步重复候选区域确定装置，所述初步重复候选区域确定装置与所述高频Kmer和低频Kmer确定装置相连，用于基于所述至少一个高频Kmer，确定初步重复候选区域，其中，基于所述低频Kmer与相邻高频Kmer之间的距离，确定是否将所述低频Kmer整合至所述初步重复候选区域中；

所述初步重复候选区域确定装置适于执行以下操作：

当所述低频Kmer与相邻高频Kmer之间的距离小于预定阈值时，将所述低频Kmer的序列整合到所述初步重复候选区域中；

重复单元确定装置，所述重复单元确定装置与所述初步重复候选区域确定装置相连，用于基于所述初步重复候选区域确定重复单元，所述重复单元确定装置包括：

延伸单元，所述延伸单元通过对所述初步重复候选区域进行延伸处理，得到经过延伸的初步候选区域；

首个高频Kmer确定单元，所述首个高频Kmer确定单元与所述延伸单元相连，用于沿着所述经过延伸的初步候选区域的预定顺序，确定首个高频Kmer，并确定所述首个高频Kmer在所述经过延伸的初步候选区域中的位置；

重复单元确定单元，所述重复单元确定单元与所述首个高频Kmer确定单元相连，用于基于下游Kmer的频数，沿着所述首个高频Kmer的位置向下游进行延伸，直到遇到频数实质变化Kmer，停止所述延伸，以便获得至少一个重复单元。

12.根据权利要求11所述的系统，其特征在于，所述生物序列为氨基酸序列或者核酸序列。

13.根据权利要求12所述的系统，其特征在于，所述生物序列为核酸序列，所述Kmer序列的长度为10个核苷酸。

14.根据权利要求12所述的系统，其特征在于，所述生物序列为氨基酸序列，所述Kmer序列的长度为3个氨基酸。

15.根据权利要求12所述的系统，其特征在于，对于氨基酸序列，所述生物序列的长度为200~500个氨基酸。

16.根据权利要求12所述的系统，其特征在于，对于氨基酸序列，所述生物序列的长度为300个氨基酸。

17.根据权利要求12所述的系统，其特征在于，对于核酸序列，所述生物序列的长度为800~1500bp。

18.根据权利要求12所述的系统，其特征在于，对于核酸序列，所述生物序列的长度为1000bp。

19.根据权利要求11所述的系统，其特征在于，所述高频Kmer和低频Kmer确定装置适于执行以下操作：

将所述全部Kmer序列每一个的频数与预定的频数阈值进行比较，以便确定所述高频Kmer和所述低频Kmer，其中，所述预定阈值是所述至少一个Kmer序列中最高频数的至少0.3倍。

20.根据权利要求11所述的系统，其特征在于，所述重复单元确定装置进一步包括评估单元，所述评估单元与所述重复单元确定单元相连，用于对所述重复单元进行评估，确定所述重复单元的PR分数。

21.根据权利要求11所述的系统，其特征在于，所述系统进一步包括一致性确定装置，所述一致性确定装置与所述重复单元确定装置相连，用于将所述至少一个重复单元进行比对，以便确定所述至少一个重复单元的一致性序列。