CN108699601A

CN108699601A - 第三代测序比对算法

Info

Publication number: CN108699601A
Application number: CN201780010771.0A
Authority: CN
Inventors: W·H·王; P·T·阿夫沙尔
Original assignee: Leland Stanford Junior University
Current assignee: Leland Stanford Junior University
Priority date: 2016-02-11
Filing date: 2017-02-10
Publication date: 2018-10-23
Also published as: EP3414348A4; WO2017139671A1; US20190042696A1; EP3414348A1

Abstract

公开用于将读取序列与参考序列比对的方法、软件和系统。在某些实施例中，所述方法、软件和系统涉及确定在所述读取序列的区域和所述参考序列的区域之间的k‑聚体的分布的相似性以便确定所述读取序列的所述区域是否映射到所述参考序列的所述区域。

Description

第三代测序比对算法

交叉引用

本申请要求2016年2月11日提交的美国临时专利申请第62/294,205号的权益，所述申请以全文引用的方式并入本文中。

关于联邦赞助的研究或开发的声明

本发明是在政府支持下在由美国国家卫生研究院(National Institutes ofHealth)授予的合同R01HG007834下进行。政府对本发明拥有一定的权利。

背景技术

全基因组测序已彻底改变生物和医学驱动的全面表征DNA序列变化、多种物种的重新测序、微生物群落的测序、检测基因组的甲基化区域、定量转录丰度、表征存在于给定样品中的基因的不同同工型、识别mRNA转录物有效地平移的程度等。实际上，药物基因组学领域由于患者基因组序列信息的增加的可获得性而以指数方式扩增。

第一和第二代测序技术以相对低成本提供巨大吞吐量。第三代测序(TGS)技术为基于单分子测序(SMS)的测序方的下一种重要技术。与第一和第二代测序技术相比，TGS工具产生较长读段，但是测序其受主要呈插入和缺失(插入缺失)形式的较高错误率困扰。

测序DNA的过程包含三个基本阶段，包括样品制备、物理测序和任选地比对，和/或重新组装。样品制备涉及使测序的基因组片段化和扩增片段。在测序期间，依次识别在每个片段中各个碱基，创建各个读段。然后利用包含算法的生物信息学软件以比对重叠的读段，这允许原始基因组组装成连续序列。

目前，用于将各个长读段与参考序列或数据集比对的常用算法基于种子和延伸概念的修改型式。这类方法通常通过寻找在查询和参考序列之间的精确匹配起始，然后大量寻找理想种子链并且使用动态编程借助任选的急下降启发法将其延伸以避免在差区域上延伸。

在本公开中提供的方法、软件和系统提供鲁棒的定位读段的测序位置的方法，从而实现比对和组装可包含畸变(如插入和/或缺失)的序列读段。

发明内容

本公开提供用于将读取序列与参考序列比对的方法、系统、可实行软件产品和存储装置。在某些实施例中，公开用于将读取序列与参考序列片段比对的方法。方法可包含创建用于读取序列的窗口和用于参考序列片段的，其中窗口具有相同长度；计算在每个窗口内独特k-聚体出现的数量，基于在每个窗口内独特k-聚体出现的数量，计算k-聚体计数相似性值；对于跨读取序列的多个窗口和跨参考序列片段的多个窗口迭代地执行步骤(a)到(c)，由此计算多个k-聚体计数相似性值，其中在读取序列和参考序列片段中的每一个中的每个随后窗口的开始从在相应序列中的先前窗口的开始偏移距离d；通过求多个k-聚体计数相似性值的平均值，计算相似性评分；和当相似性评分高于阈值时，将读取序列与参考序列片段比对，其中将在步骤(a)的第一次执行中创建的窗口放置在每个序列的起始处。

在某些实施例中，方法可包含对于读取序列和参考序列的不同片段重复步骤(a)到(f)。

在某些实施例中，参考序列片段可为从基因组数据库获得的参考序列的区域。在某些实施例中，参考序列可为读取序列。在某些实施例中，参考序列可为从测序与获得读取序列的序列相同样品获得的读取序列。

在某些实施例中，窗口中的每一个的长度可为至少50个碱基。在某些实施例中，窗口中的每一个的长度可为在1-10,000个碱基范围内的任何整数值，其中长度保持恒定。

在某些实施例中，距离d可为至少10个碱基长。在某些实施例中，距离d的长度可在1-500个碱基范围内，其中d保持恒定。

在某些实施例中，k-聚体的长度可为2-10碱基。在某些实施例中，k-聚体的长度可为3个碱基。在某些实施例中，k-聚体的长度可为4个碱基。

此外本文公开存储在计算机可读媒体上的可实行软件产品。在某些实施例中，存储在计算机可读媒体上的可实行软件产品可含有用于进行上文所公开方法的程序指令。

还提供被配置成实行指令以进行上文所公开方法的系统。系统可包含具有进行上文所公开方法的存储的指令的存储器和耦合到存储器并且被配置成实行在存储器中的指令的处理器。

在某些实施例中，公开存储可实行用于执行上文所公开方法的指令的存储装置。

附图说明

图1描绘参考序列的参考序列片段和用于参考序列片段和用于读取序列的例示性窗口。

图2描绘用于在参考序列片段的窗口内和在读取序列的对应的窗口内计数k-聚体的实施例。

图3描绘在参考序列片段和读取序列中的多个窗口。

图4描绘用于将读取序列与参考序列的多个片段比较的示意图。

图5描绘跨参考序列比对读取序列的计算的相似性评分。

图6说明用于进行所公开方法的计算机的一个实施例。

图7描绘在大肠杆菌(E.coli)基因组中随机位置之间在其突变(仅替代)型式的情况下的余弦距离的分布，使用k＝3。

图8为图7的续图。

图9描绘在大肠杆菌基因组中随机位置之间和在其突变(仅替代)型式的情况下的余弦距离的分布，使用k＝4。

图10为图9的续图。

图11描绘来自大肠杆菌基因组的长度5000个碱基的1000个随机序列之间和相对于其突变型式(替代和插入缺失)的余弦距离的分布，使用k＝3。

图12描绘来自大肠杆菌基因组的长度5000bps的1000个随机序列之间和相对于其突变型式(替代和插入缺失)的余弦距离的分布，使用k＝4。

图13描绘跨全基因组比较的来自大肠杆菌基因组的长度5000个碱基的读取序列的余弦相似性评分，其中模拟错误率为15％和35％，其中k＝3并且d＝10。

图14集中在来自图13的预期比对位置(竖直点线)附近并且缩放，描绘与取样位置附近的大肠杆菌基因组相比，来自大肠杆菌基因组的长度5000个碱基的读取序列的余弦相似性评分，其中模拟错误率为15％和35％，其中k＝3并且d＝10。

定义

本文中(不论上文或下文)所列举的所有出版物、专利和专利申请都以全文引用的方式并入在此。

在描述本发明时，将采用以下术语并旨在如下文所指示定义。

必须注意，除非上下文另有清楚地规定，否则如在本说明书和所附权利要求书中所用，单数形式“一(a/an)”和“所述(the)”包含多个指示物。还要注意，权利要求书可经起草而排除任何任选的元件。因而，此陈述旨在与对所要求保护的元件的引述结合而充当使用如“仅仅(solely)”、“仅(only)”等这类排他性术语或使用“负性”限制的前提基础。

如本文所用，术语“比对”或其语法等效物是指将读取序列映射到在参考序列中的区域。

如本文所用，术语“读取序列”是指通过测序仪器从样品核酸的单个片段确定的连续核苷酸的序列。单个片段可为通过扩增待测序的基因组或基因组的一部分产生的扩增产物。来自样品核酸的单个片段的连续核苷酸的序列可表示为通过测序技术产生的数据流，数据例如借助于与测序技术相关联的碱基呼叫软件产生。例如，来自DNA测序平台的商业供应商的碱基呼叫软件。读取序列还可被称作“查询序列”或“序列读段”。

如本文所用，术语“参考序列”是指生物体的基因组或基因组的一部分的连续核苷酸的已知序列。参考序列可用作读取序列与其比对的输入序列。待使用的参考序列取决于读取序列的来源。参考序列可为来自与获得读取序列的物种相同的物种的核酸的序列。如果来自相同物种的序列生物体不可用，那么与其基因组正在被测序的生物体最紧密相关的生物体的序列可用作参考序列。参考序列可通过测序技术确定或可从序列数据库(如从美国国家生物技术信息中心(National Center for Biotechnology Information)的基因组文库获得的生物体的基因组)获得。参考序列还可为读取序列。读取序列从测序核酸样品获得的情况下，将读取序列与读取序列比对适用于在读取序列中寻找重叠区域和组装读取序列以产生较长连续读取序列。

如本文所用，术语“数据结构”是指通常在计算机或存储器装置中的信息组织。数据结构允许高效实行处理信息/数据的算法。例示性数据结构包含词典、队列、栈、链表、堆栈、杂凑表、阵列、树等。数据结构可具有对应于信息的单元或相关信息的子集的子结构。举例来说，阵列具有项的行和列；树具有节点、分支、子树和叶子；等。例示性数据结构可包含所有可能的独特k-聚体的列表和列入读取和参考序列中的独特k-聚体出现的数量的计数指示。

如本文所用，在两个序列的情况下，术语“同一性”是指分别两个多核苷酸或多肽序列的精确核苷酸与核苷酸或氨基酸与氨基酸的对应关系。同一性百分比可通过以下确定：通过比对序列直接比较在两个分子之间的序列信息，计数在两个比对的序列之间的精确匹配数量，除以较短序列的长度，并且将结果乘以100。现成的计算机程序可用于辅助分析，如ALIGN,Dayhoff,M.O.在《蛋白质序列和结构地图集(Atlas of Protein Sequenceand Structure)》M.O.Dayhoff编,5月增刊3:353-358,美国生物医学研究基金会(Nationalbiomedical Research Foundation),华盛顿特区(Washington,DC)中，其将Smith和Waterman《应用数学进展(Advances in Appl.Math.)》2:482-489,1981的局部同源算法改编用于肽分析。用于确定核苷酸序列同一性的程序在威斯康星序列分析包(WisconsinSequence Analysis Package)，第8版(可购自威斯康星州麦迪逊的遗传学计算机集团(Genetics Computer Group,Madison,WI))中可用，例如BESTFIT、FASTA和GAP程序，其也依赖于Smith和Waterman算法。用由制造商推荐和在上文提及的威斯康星序列分析包中描述的默认参数容易地利用这些程序。举例来说，特定核苷酸序列与参考序列的同一性百分比可用六个核苷酸位置的默认评分表和空位罚分使用Smith和Waterman的同源算法确定。

在本文使用术语“多核苷酸”、“核酸”和“核酸分子”以包含核糖核苷酸或脱氧核糖核苷酸的任何长度的核苷酸的聚合形式。此术语仅指分子的一级结构。因此，所述术语包含三链、双链和单链DNA以及三链、双链和单链RNA。它还包含修饰，如通过甲基化和/或通过盖帽，和未修饰形式的多核苷酸。更具体地，术语“多核苷酸”、“核酸”和“核酸分子”包含多脱氧核糖核苷酸(含有2-脱氧-D-核糖)、多核糖核苷酸(D-核糖)、为嘌呤或嘧啶碱基的N-或C-糖苷的任何其它类型的多核苷酸，和含有非核苷酸主链的其它聚合物。

如本文所用，“靶核酸”或“靶核苷酸序列”是指其中核苷酸序列待确定的所关注的任何核酸。

具体实施方式

本公开提供用于将读取序列与参考序列比对的方法、软件和系统。

应了解，为了清楚起见，在分开的实施例的情况下描述的本发明的某些特征还可以组合在单个实施例中提供。相反，为了简便起见，在单个实施例的情况下描述的本发明的各种特征也可分开或以任何合适子组合提供。具体地说，实施例的所有组合被本发明包括并且在此公开，就如同每个和每一个组合被单独且明确公开一样，到这类组合包括可操作过程和/或装置/系统的程度。此外，在描述这类变量的实施例中所列的所有子组合也具体地被本发明包括并且在此公开，就如同每个和每一个化学基团的这类子组合在此单独且明确公开一样。

方法

本公开提供用于将读取序列与参考序列的区域比对的方法。读取序列也被称作查询序列。比对方法可涉及(a)创建用于读取序列的窗口和用于参考序列的片段的窗口，所述窗口具有相同长度；(b)计算在每个窗口内独特k-聚体出现的数量，其中k-聚体具有相同长度；(c)基于在每个窗口内独特k-聚体出现的数量，计算k-聚体计数相似性值；(d)对于跨读取序列的多个窗口和跨参考序列的片段的多个窗口迭代地执行步骤(a)到(c)，其中在读取序列和参考序列的片段中的每一个中的随后窗口的开始从在相应序列中的先前窗口的开始偏移距离d，其中d在读取序列和参考序列中的对应的窗口之间相同；(e)通过求计算的多个k-聚体计数相似性值的平均值，计算相似性评分；和(f)如果相似性评分高于阈值，那么将读取序列与参考序列的片段比对，其中将在步骤(a)中创建的窗口放置在读取序列和参考序列的片段的起始处。

在某些实施例中，创建窗口的步骤(a)可涉及放置在读取序列和参考序列的片段中的每一个的第一核苷酸处起始的窗口。在某些实施例中，创建初始窗口下游的额外窗口的步骤可涉及选择放置额外窗口的在读取序列和参考序列片段中的区域或子序列。举例来说，在读取序列和参考序列的片段中的每一个中的额外窗口可从它紧接着上游的窗口偏移距离d，所述距离d可为约1或更多个碱基。对于每个窗口，偏移距离d可保持恒定。换句话说，在读取序列和参考序列的片段中的每一个中的窗口从先前窗口偏移相同距离。

窗口的长度/大小可由w表示，其可在1-10,000个碱基，例如100-10,000个碱基、10-5000个碱基、50-1000个碱基、100-1000个碱基、100-800个碱基、100-700个碱基、50-1000个碱基、50-800个碱基、50-700个碱基、50-500个碱基、100-500个碱基、300-700个碱基、400-700个碱基或400-600个碱基的范围内。在某些实施例中，对于在读取序列和参考序列片段之间的单个比对，窗口大小w为恒定的。换句话说，对于单个比对创建的所有窗口可具有相同长度。在一些情况下，读取序列和参考序列片段的长度可相似。在其它情况下读取序列和参考序列片段可具有相同长度。窗口可用于表示序列[i,i+w]的区域，其中i为表示在读取序列或参考序列片段中的位置的索引整数，并且w为窗口长度。图1说明示出其中选择片段用于与读取序列比较的参考序列的示意图的实例。在图1中描绘参考序列的片段(灰色区域)。图1还示出在位置i处起始的长度为10个碱基的窗口(通过方括号表示)。对于读取序列，相似地创建在位置i处起始的长度为10个碱基的对应的窗口。在此实例中，读取序列和参考序列片段不相同。应注意，10个碱基的窗口大小仅出于说明的目的。如在本文中所提及，放置在所描绘的窗口下游的随后窗口的长度w保持恒定。

在某些实施例中，可通过计数和保持每一种可能独特k-聚体的每个情况的轨迹计算在每个窗口内每个可能独特k-聚体出现的数量(也被称作k-聚体分布或k-聚体计数分布)。在某些实施例中，在读取序列中的窗口中的核苷酸序列可用于产生所有重叠k-聚体的列表，并且在参考序列的片段中的对应的窗口(在相同位置i处起始)中的核苷酸序列也可用于产生所有重叠k-聚体的列表。对于每个窗口可计数独特k-聚体的数量以确定在每个窗口中独特k聚体出现数量的相似性。在其它实施例中，数据结构可用于计数独特k-聚体。在一对窗口中每个独特k-聚体出现数量的相似性越高，k-聚体计数相似性值越高。

k-聚体的长度可为3、4、5、6、7、8、9、10、20、30、40、50、60、70、80、90或100个碱基；例如，k-聚体的长度可在3-100、3-80、3-50、3-80、3-10、3-5、4-100、4-80、4-50、4-80、4-10、4-5、2-4、2-10或3-4个碱基的范围内。如在本文中所提及，在窗口内，k-聚体的大小保持恒定。在某些实施例中，在为比对方法产生的所有窗口中，k-聚体大小保持恒定。在每个窗口内，连续k-聚体可重叠至少1个碱基、至少2个碱基、至少3个或最多k-1个碱基(例如，对于10个核苷酸长的k-聚体，连续k-聚体可重叠最多9个碱基)。在某些实施例中，跨在读取序列中的窗口和跨在参考序列的片段中的对应的窗口在连续k-聚体之间的重叠为恒定的。在某些实施例中，对于整个比对方法，在相邻k-聚体之间的重叠长度为恒定的。举例来说，连续3-聚体可重叠1或2个碱基，并且连续4-聚体可重叠1、2或3个碱基。在某些实施例中，连续3-聚体可重叠2个碱基，并且连续4-聚体可重叠3个碱基。在某些实施例中，连续k-聚体可不重叠。举例来说，连续k-聚体可被1-3000个核苷酸，如，50-1000个碱基、100-1000个碱基、100-800个碱基、100-700个碱基、50-1000个碱基、50-800个碱基、50-700个碱基、50-500个碱基、100-500个碱基、300-700个碱基、400-700个碱基或400-600个碱基分开。在某些实施例中，跨整个窗口，k-聚体大小可为恒定的，并且可计数跨整个窗口的k-聚体。举例来说，如图2所示，对于10个碱基的窗口长度，计数与先前k-聚体重叠3个碱基的所有4-聚体，将跨用于读取序列和用于参考序列的片段的窗口的整个长度计数七个4-聚体。

图2示出在对于读取序列和参考序列片段重叠3个碱基的4-聚体的情况下，读取序列的长度为10(用方括号示出)的第一窗口，其中相对于参考序列片段，读取序列在位置9处插入有核苷酸“A”(通过*指示)。将所有计数的4-聚体加下划线，并且将窗口内独特k-聚体的计数列入图2中示出的表中。应注意，如在此实例中，当k＝4时，存在4⁴或256个独特k-聚体。因此，对于不存在于图2中的表中的248个其它独特k-聚体，读取序列和参考序列片段将具有0计数。

在存在插入和缺失的情况下，比对的序列使其无错误片段移位偏移o，其中o＝|ins-del|。尽管由这类插入和缺失引起的移码通常对使用用于较长读取序列的直接交叉相关的比对技术的准确度是不利的，但计数k-聚体通过分析短k-聚体的含量作为局部相似性的量度而绕过此问题。以下文读取序列和参考序列片段为实例，

其中n为任何核苷酸，*为缺失，粗体字母为插入，并且竖直线限定窗口的边界(长度17)。在两个序列中简单搜索对应的相同、连续核苷酸的方法将最多发现CCCCGG。然而，本文公开的方法确定k-聚体出现的数量，即，k-聚体的分布，导致较大局部比对。

在读取序列的指定窗口中的k-聚体如下：XAAA、AAAA、AAAG、AAGC、AGCC、GCCC、CCCC、CCCG、CCGG、CGGT、GGTT、GTTT、TTTT、TTTX。在参考序列的对应或平行窗口中的k-聚体如下：XAAA、AAAA、AAAC、AACC、ACCC、CCCC、CCCG、CCGG、CGGA、GGAT、GATT、ATTT、TTTT。

在此实例中加下划线的k-聚体出现在读取序列和参考序列两者中一次。本实例说明甚至当读取序列与参考序列不相同时，如何跨整个窗口识别相同片段和如何用于将读取序列映射到参考序列的区域。

在窗口大小接近较大数量(如500)时，本文公开的方法变得相对于插入缺失更宽容，因为尽管插入和缺失增加它们之间的偏移，但相关k-聚体可仍然处于对应的窗口内。

识别k-聚体出现数量或k-聚体计数的此方法使比对方法在面对插入缺失时比其它交叉相关比对方法更鲁棒。

在一些情况下，可通过分别对于读取序列和参考序列片段创建k-聚体计数向量和计数独特k-聚体的出现数量，其中i为在其中窗口起始的序列中的位置，和

V_i＝[n₀,n₁,…,n_q]，其中0≤i＜l-w+1, [1]

其中在窗口中，n_j＝计数(v_j)

其中对于0≤j＜q，Q＝k-聚体的列表{v_j}，和

q＝||^k，其中＝{A,T,G,C}。

举例来说，当对于DNA序列计数所有独特3-聚体的出现数量时，那么k＝3，q＝64总计可能独特3-聚体。

在某些实施例中，包含未知的(一个或多个)碱基的k-聚体可随机映射到q种可能性中的一种，由此保持可能的独特k-聚体池限于的组合。在其中存在包含未知的(一个或多个)碱基的k-聚体的某些实施例中，对于v_j可考虑q+1种可能性，其中q含有具有未知的(一个或多个)碱基的k-聚体并且对于读取序列和参考序列两者计数独特k-聚体。这样做将类似于对于未知的(一个或多个)k-聚体的集合将另一个维度引入k-聚体空间。

在一些实施例中，可基于在读取序列和参考序列片段的对应的窗口内独特k-聚体出现的数量，计算k-聚体计数相似性值。可通过在读取序列和参考序列片段的k-聚体计数向量之间使用以下余弦相似性公式，计算k-聚体计数相似性值(还可被称作k-聚体分布相似性值或k-聚体计数分布相似性值)：

此k-聚体计数相似性值或评分(S_i)表示分别在读取序列和参考序列片段中在局部比对位置处序列片段f¹和f²的局部相似性，其中片段f¹由窗口限定并且f²由对应的窗口限定。与其它度量相比，使用余弦相似性评分提供可使用快速傅里叶变换(FFT)高效实施全域相似性评分(方程4)的优点。在一些实施例中，可使用其它相似性度量，如欧几里德距离(Euclidean distance)。

在某些实施例中，对于跨读取序列和参考序列片段的多个窗口创建窗口、计数独特k-聚体和迭代地计算k-聚体计数相似性值的步骤提供多个k-聚体计数相似性值。创建窗口、计数独特k-聚体和计算k-聚体计数相似性值的步骤可执行到读取序列的整个长度已与参考序列的片段比较。在某些实施例中，创建窗口、计数独特k-聚体和计算k-聚体计数相似性值步骤可进行到读取序列的整个长度已与另一个读取序列比较。在某些实施例中，创建窗口、计数独特k-聚体和计算k-聚体计数相似性值的步骤可进行到读取序列的至少500个核苷酸长的段已与参考序列比较。因而，在某些实施例中，可计算读取序列的至少500个碱基长的段的k-聚体计数相似性值，例如700个碱基、1000个碱基、3000个碱基、5000个碱基、7000个碱基、10,000个碱基、13,000个碱基、15,000个碱基、18,000个碱基、20,000个碱基或至多50,000个碱基。虽然迭代地执行计算k-聚体相似性值的步骤，但是对于不同对窗口(在读取和参考序列中对应的窗口)不必依次计算k-聚体相似性值。在某些实施例中，对于多对窗口，可同时执行步骤(a)到(c)。如在本文中所提及，在读取序列和参考序列片段中的每一个中的每个随后窗口从在相应序列中的先前窗口的开始偏移距离d。此窗口偏移距离d大小可变化并且大小的长度可在1-1000个碱基、1-500个碱基、1-100个碱基、5-100个碱基、10-100个碱基、50-100个碱基、1-20个碱基或1-10个碱基、10-50个碱基范围内。在某些实施例中，对于整个比对，在窗口之间的偏移距离d为恒定的。在某些情况下，相邻窗口可重叠。在某些情况下，相邻窗口可不重叠。在某些情况下，未重叠的相邻窗口可彼此紧邻，即，没有被间插核苷酸分开。在某些实施例中，窗口大小可为500个碱基，并且偏移d可为10个碱基。在其它实施例中，窗口大小和偏移的不同组合是可能的，如窗口大小为50个碱基与偏移为5个碱基，窗口大小为100个碱基与偏移为10个碱基，窗口大小为500碱基与偏移为50个碱基，窗口大小为1000个碱基与偏移为100个碱基，或窗口大小为5000个碱基与偏移为100个碱基。图3示出从在位置i(括号)处起始的先前窗口偏移距离d＝5的随后窗口(粗体方括号)，其中对于读取序列和参考序列片段两者，w＝10。在此实例中，读取序列和参考序列片段不相同。应注意，窗口大小为10个碱基并且偏移为5仅出于说明的目的。

如在本文中所提及，在读取序列和参考序列片段中的对应的窗口之间，d相同或相等。举例来说，对于读取序列和参考序列片段中的每一个创建的第二窗口可从相应序列中的每一个中的第一窗口偏移10个碱基。如所解释，使用上述方法迭代地创建窗口并且对于每个窗口计算k-聚体计数相似性值。在某些实施例中，通过求多个k-聚体计数相似性值的平均值来计算读取序列与参考序列片段之间的总相似性评分，所述k-聚体计数相似性值可使用计算相似性评分(S_i)的公式[2]计算。通过以下方程计算总相似性评分：

其中

γ＝{(i¹,i²),(i¹+d,i²+d),…,(i¹+pd,i²+pd)}

并且i¹和i²分别表示读取序列和参考序列片段的起始索引。使用余弦相似性作为度量而跨参考序列计算方程3可公式化为交叉相关并且可使用FFT来高效计算。

在一些情况下，对于读取序列和参考序列的不同片段或区域，可迭代地重复相似性评分的计算。图4示出读取序列可如何与在参考序列中的多个片段比较。在其中读取序列与在参考序列中的多个片段比较的某些方面，对于沿参考序列读取序列的每一个放置，可比较读取序列，并且随后计算相似性评分；换句话说，比较读取序列与在参考序列的每一个取样位置i处起始的参考序列，其中0≤i＜l_参考-l_读取+1。在这类实施例中，读取序列和参考序列可被认为在从参考序列的起始的比对偏移m处进行比较。在这类实施例中，可通过以下计算余弦相似性评分

其中m＝i²-i¹并且和分别为在读取和参考序列中长度为w的比对的第一窗口的起始位置。通过以下计算在偏移m的情况下比对的两个序列之间的总相似性评分

其中

并且出于全域比对的目的，可被设定成0，并且

对于可使用快速傅里叶变换(FFT)高效计算方程4以寻找读取序列对参考序列的总相似性评分(全域比对)。对于长度l的序列，对于i＜l^V,(l^V＝l-w-1，即对于k聚体计数向量的起始位置的有效长度)，V_i如在方程1中所定义并且否则V_i＝[0,…,0]。对于i＝jd，因为计算分开d的相邻窗口的V_i。

其中

★：交叉相关运算

*：卷积运算。

DFT()_N：大小为N的区傅里叶变换

IDFT()_N：大小为N的逆区傅里叶变换

可使用快速傅里叶变换(FFT)算法高效计算DFT。对于与实际DFT大小(N)比较较大的l_max，可使用重叠相加或重叠保留技术。

在存在插入和缺失的情况下，在偏移m处比对的序列根据之前插入和缺失碱基的数量，具有在不同偏移下移位的无错误片段其原因是，代替观测在固定偏移m下核苷酸碱基的正合序列匹配，比较在每个碱基位置处起始的窗口中短k-聚体的含量作为局部相似性的量度。

在某些实施例中，如果相似性评分高于阈值，那么读取序列可与参考序列的片段或区域比对。在某些实施例中，阈值可为高于平均值或中位值至少1.5倍标准差(SD)或中值绝对偏差(MAD)的值，如SD或MAD的2倍或3倍。在某些实施例中，在其中对于读取序列和参考序列的不同片段已经计算许多相似性评分(S)的情况下，可使用以下公式计算阈值：平均(S)+f×std(S)或中值(S)+f×mad(S)。在一些情况下，f＝1、2或3。图5描绘读取序列与参考序列的比对。在图5中，高于阈值的在读取序列和参考序列的片段之间的相似性评分可见为峰，指示读取序列映射到参考序列的片段。在某些实施例中，当相似性评分低于阈值时，读取序列可不与参考序列的片段或区域比对。在这类实施例中，比对方法可包含对于参考序列的不同片段进行步骤(a)到(f)。

在某些实施例中，迭代地执行方法直到读取序列的整个序列已与参考序列的片段比较。在某些实施例中，迭代地执行方法直到读取序列的整个序列已与整个参考序列比较(例如当参考序列为另一个读取序列并且比较读取序列以识别重叠读取序列时)。

在某些实施例中，读取序列被分成较短的序列，并且对较短的序列执行所述方法。在某些实施例中，长度为7000个碱基或更多的读取序列可被拆分成2或更多个同样大小的(在可能时)子序列或片段。在其它实施例中，长度为5000个碱基或更多、6000个碱基或更多、8000个碱基或更多或10,000个碱基或更多的读取序列可被拆分成子序列。在某些实施例中，本文所描述的方法可使用已被分成约1000-7000个碱基(如1000-2000个碱基、1000-3000个碱基、1000-4000个碱基、1000-5000个碱基或1000-6000个碱基)的较短序列的读取序列执行。在某些情况下，疑似包含插入和/或缺失的读取序列可被分成较短序列。

在其它实施例中，长度l>2000、3000、4000、5000、6000、7000和8000个碱基读取序列可被拆分成子序列。这样做是因为随读取序列长度增加，读取序列相对于参考序列的绝对值(#插入-#缺失)可变得与朝向读取结束的窗口长度w相当，并且在读取查询窗口和参考序列窗口之间的交叉相关变得无效和嘈杂。在其中读取序列被拆分成子序列的这类实施例中，原始读取序列的每个子序列单独地与参考序列比对，对于每个子序列，重复创建窗口、计数k-聚体、计算k-聚体计数相似性值、计算相似性评分和将子序列与参考序列片段比对的步骤。

在另外的实施例中，方法进一步包括将读取子序列合并回一起以确认读取序列与参考序列的一个比对。在某些实施例中，每个读取子序列在峰位置处与参考序列的区域比对。来自读取子序列的相容的峰位置合并回一起。在一些情况下，在范围([p-o,p+l+o])内的读取序列和选择的参考序列片段之间使用带状动态编程(BDP)计算顶部选择的峰的精确起始位置，其中p为检测的峰位置，l为读取长度，并且o为由于峰位置检测不准确度考虑的界限。在某些实施例中，o＝2×d。在一些情况下，对于BDP默认评分设定为：匹配：+5，错配：-4，开放空位：-10，延伸空位：-1。如对于某些实施例所提到，基于跨参考序列的S的平均和标准差检测顶部峰。在某些实施例中，在BDP平台中也考虑在最大峰值g×std(S)内的其它位置，至多N_max位置的最大数量，其中g>0并且N_max可在1-1000个峰的范围内。如果对于读取没有检测到相当大的峰，那么在BDP平台中选择顶部N_max用于合并读取子序列。在所公开的方法中的分析步骤可以任何合适的编程语言实施，如C、C++、Java、C#、Fortran、Pascal等。

基于在局部窗口内的k-聚体的计数计算评分的优点为双重的。第一，算法对插入缺失鲁棒，因为它不寻求在序列的长段之间的精确匹配。相反，记录重叠k-聚体的计数，使插入和缺失将读取序列片段从对应的参考序列片段移位，而不损失在那些片段之间的相似性信息。第二，可通过组合计算比对相似性评分的方法和带状动态编程技术以将读取子序列拼接在一起来准确比对相对长读取查询。使用交叉相关的这类鲁棒和准确比对方法提供种子和延伸算法的有用的替代方案，其试图找到在查询和参考之间精确匹配的集群，用于将序列准确和高效地映射到较大数据库。

在某些实施例中，方法为计算机实施方法。在某些实施例中，算法和/或结果(例如在读取和参考序列之间的理想比对)存储在计算机可读媒体上，和/或显示在屏幕上或在纸印刷品上。在某些实施例中，进一步分析结果，例如以识别基因变异体，以识别序列信息的一个或多个来源，以识别在个体或物种之间保留的基因组区域，或确定在两个个体之间的亲缘关系。

在计算机上实施的所公开的方法的功能方面可使用任何适当实施环境或编程语言(如C、C++、Cobol、Pascal、Java、JavaScript、HTML、XML、dHTML、程序集或机器代码编程、RTL等实施或完成。

在某些实施例中，计算机可读媒体可包括硬盘驱动器、辅助存储器、外部存储器、服务器、数据库、便携式存储器装置(CD-R、DVD、ZIP盘、闪存卡等)等的任何组合。

读取和参考序列

用于测序的任何高吞吐量技术可用于实践本文公开的方法。DNA测序技术包含使用标记的终止子或引物和在平板或毛细管中的凝胶分离的双脱氧测序反应(桑格法(Sanger method))、使用可逆终止的标记核苷酸的边合成边测序、焦磷酸测序、454测序、使用等位基因特定杂交到标记克隆文库随后接合的边合成边测序、在聚合步骤期间并入标记的核苷酸的实时监视、聚合酶克隆测序、SOLID测序等。这些测序方法因此可用于测序所关注的靶核酸和获得查询读取序列。参考序列可同样测序(例如参考序列可为对从核酸样品的测序获得的其它(一个或多个)读取序列比对的读取序列)或可通过公用数据库(如国家DNA数据库)获得，并且可采取一个或多个序列的形式，比如基因组。在一些实施例中，参考序列为用于在参考数据库(如)中的靶核酸的序列。

可获得针对任何受试者的核酸的读取序列。受试者可为生物体，如单细胞生物体(例如细菌、古细菌、原虫、单细胞海藻和单细胞真菌)或多细胞生物体(例如海绵、刺细胞动物、扁虫、节肢动物、棘皮动物、脊索动物、脊椎动物、蕨类植物、被子植物和裸子植物)。在某些情况下，读取序列可从感染性生物体、病原体(如，奈瑟菌属(Neisseria)、HIV、大肠杆菌、沙门氏菌(Salmonella)等)获得。

读取和参考序列可从相同物种、亚种、菌株或最紧密相关的生物体获得。举例来说，来自人类的读取序列可与来自另一个人类的参考序列(如人类基因组的版本)比较。在某些实施例中，(一个或多个)参考序列可来自与来自获得读取序列的生物体进化或生物紧密相关的生物体，使得可实现高比对准确度。

在某些实施例中，所公开的方法可应用于寻找读取重叠(即读取序列的成对比对)。在这类情况下，参考序列将为另一个读取序列。

如本文所论述，读取序列为通过测序仪器从样品核酸的单个片段确定的连续核苷酸的序列。在某些实施例中，读取序列不通过组装具有重叠区域的单独读取序列预组装，在重叠区域处核苷酸序列高度相似或相同。换句话说，读取序列可为从测序由生物体的基因组产生的单个核酸片段获得的连续核苷酸的序列。读取序列长度可在1-20,000个碱基、1-15,000个碱基、50-15,000个碱基、100-15,000个碱基、100-10,000个碱基、100-9000个碱基、100-8000个碱基、100-7000个碱基、100-6000个碱基、100-5000个碱基、100-2500个碱基、500-10,000个碱基、500-7500个碱基、500-5000个碱基和500-2500个碱基长度范围内变化。

计算软件和系统

如在本文中所提及，在本申请中提供的方法可在硬件和/或软件中实施。在一些实施例中，方法的不同方面可在任一客户端逻辑或服务器侧逻辑中实施。在某些情况下，用于实施所公开的方法的组件可体现在含有逻辑指令和/或数据的固定媒体编程组件，当逻辑指令和/或数据加载到适当配置的计算装置时引起装置执行方法步骤。含有逻辑指令的固定媒体可递送到查看者的计算机或含有逻辑指令的固定媒体可驻存在查看者通过通信媒体访问以便下载编程组件的远程服务器上。

在某些方面，本文所提供的方法为计算机实施方法，其中方法的至少一个或多个步骤通过计算机程序进行。用于实施本发明计算机实施的方法的计算机系统可包含如常用于所属领域的组件的任何布置。在特定实施例中，所公开的方法可整体或部分体现为记录在固定媒体上的软件。计算机系统可为包含存储器、处理器、输入和输出装置(I/O)、数据储存库、网络接口、存储装置、电源等的任何电子装置。存储器或存储装置可被配置成存储指令，所述指令使得处理器能够通过处理和实行存储在存储器或存储装置中的指令实施本发明计算机实施的方法。计算机还可包含用于有线和/或无线通信的网络接口。

处理器控制计算机的操作并且可读取来自存储器和/或数据储存库的信息并且相应地实行实施前述实施例的指令。术语“处理器”旨在包含一个处理器、多个处理器或具有多个核的一个或多个处理器。

I/O可包含任何类型的输入装置，如键盘、鼠标、麦克风等，和任何类型的输出装置，例如，如监视器和打印机。在其中计算机包括服务器的实施例中，输出装置可耦合到本地客户端计算机。

存储器可包括任何类型的非暂时性、静态或动态存储器，包含闪存、DRAM、SRAM等。存储器可存储可用于如本文所描述的序列比对方法的程序和数据。

数据储存库可存储包含存储读取序列、参考序列、k-聚体计数向量等的一个或多个数据库的若干数据库。在一个实施例中，数据储存库可驻存在计算机内。在另一个实施例中，数据储存库可经由网口或外部驱动器连接到计算机。数据储存库可包括分开的服务器或任何类型的存储器存储装置(例如光盘或磁性媒体、固态动态或静态存储器等)。数据储存库可任选地包括多个辅助存储器装置，例如用于输入序列(例如读取序列)、序列信息、计算结果和/或其它信息的分开的存储。如在所属领域中理解，其后计算机可使用到直接服务器或用户端逻辑的信息以体现所公开的方法的方面。

在操作中，操作员可经由呈现在显示屏上的用户界面与计算机交互以指定通过各种软件程序要求的读取序列和其它参数。一旦调用，通过处理器实行在存储器中的程序以实施本发明方法。

在计算机实施的方法的一个实施例中，用户可访问在计算机系统上的文件，其中文件含有(一个或多个)读取序列和(一个或多个)参考序列数据，以及进行所公开的方法的用户和计算机可实行方法。在另外的实施例中，过程的结果可任选地进一步包括质量信息、技术信息(例如峰特性、预期的错误率)、交替(例如第二或第三最佳)共有测定、可信度度量等。

图6说明包括其中存储用于进行所公开的方法的指令的存储器的计算机的一个实施例。计算机的处理器实行存储的指令以执行比对。此计算机系统包含用于实行存储在主存储器105中的指令的CPU 101、用于显示界面的显示器102、键盘103，和指向装置104，存储各种程序的主存储器105和存储装置，如可存储输入序列109和比对结果110的辅助存储器108。装置不限于个人计算机，但是可为用于与远程数据应用交互的任何信息电气设备，并且可包含如数字支持电视、蜂窝电话、个人数字助理等的这类装置。驻存在主存储器105和辅助存储器108中的信息可用于编程这类系统并且可表示盘动态或静态存储器等。在特定实施例中，所公开的方法可整体或部分体现为记录在此固定媒体上的软件。存储在主存储器上的各种程序可包含使用本文公开的方法将读取序列与参考序列比对的程序106。连接CPU 101、主存储器105和辅助存储器108的线可表示任何类型的通信连接。举例来说，辅助存储器108可驻存在装置内或可经由例如网口或外部驱动器连接到装置。辅助存储器108可驻存在任何类型的存储器存储装置(例如服务器或媒体，如CD或软驱)上，并且可任选地包括多个辅助存储器装置，例如用于输入序列、比对结果、结果解释的结果和/或其它信息的分开存储。

比对分析的输出可以任何方便的形式提供。在一些实施例中，输出提供在用户界面、印刷品上、在数据库中等，并且输出可呈表、图、光栅曲线、热图等形式。在某些实施例中，比对方法的实施的输出可包含用于每个读取序列与在参考序列中、在多个参考序列或另一个读取序列中的位置的比对列表。在某些实施例中，过程的结果可任选地进一步包括技术信息(例如峰特性、预期的错误率)、交替(例如第二或第三最佳)比对、可信度度量等。在将读取序列与参考序列比对过程期间和在其之后，此处理的进程和/或结果可保存到存储器和数据储存库和/或通过I/O输出用于在显示装置上显示和/或保存到额外存储装置(例如CD、DVD、蓝光、闪存卡等)、传输或印刷。

实例

以下为进行本发明的特定实施例的实例。实例仅出于说明性目的提供，并且不旨在以任何方式限制本发明的范围。

已作出努力以确保关于所使用的数字的准确度(例如量、温度等)，但允许一些实验性误差和偏差应为理所当然的。

实例1

使用大肠杆菌基因组展示余弦相似性度量的效果。

余弦相似性为用于通过测量在两个向量之间的角度的余弦确定在两个向量之间的相似性的度量。为了展示此度量的效果，长度为5000个碱基的1000个序列各自选自在大肠杆菌基因组中的随机位置。对于每个序列，在不同长度w＝50、100、500、1000和5000个碱基的不重叠窗口之间，和在每个窗口的序列和平均取代率为15％和35％的其10个随机突变型式之间计算距离余弦(1-余弦相似性)。图7-8和9-10分别呈现对于k＝3和k＝4的距离余弦分布。图7-8和9-10说明在随机位置处的短k-聚体计数向量之间的距离余弦的分布如何与其突变型式可区分。此外，正如期望，对于较高突变率，分布重叠变得相当大，较高突变率增加在窗口和其突变型式之间的不同。在图7-10的图中，左边分布为在不同突变率情况下的“在其突变型式情况下的距离”，并且右边分布为在不同突变率情况下的“随机位置之间的距离”。

除了计算在具有仅由取代组成的突变图案的单个窗口之间的余弦相似性之外，测试多个窗口和插入缺失的添加。此外，长度为5000个碱基的1000个读取序列选自在大肠杆菌基因组中的随机位置。计算在每个读取序列和10个其它随机不重叠读取序列之间，和在每个读取序列和使用窗口移位大小d＝10在存在15％和35％的错误率的情况下其10个随机突变型式之间的相似性距离(1-S[0])。使用10％取代、28％缺失和62％插入的错误图案。图11和12分别示出对于k＝3和k＝4的此测试的相似性分布。在图11和12的图中，左边分布为具有不同错误率的“在其突变型式情况下的距离”，并且右边分布为具有不同错误率的“随机位置之间的距离”。由于相对较高插入缺失率，在较小窗口大小之间的不配对对相似性评分具有负影响。相反地，使用大(w>1000个碱基)窗口大小在每个位置处失去k-聚体的地点信息。

为了找到理想窗口大小，从在取样位置2,720,230-2,725,230处起始的大肠杆菌K12区域提取读取序列并且在15％和35％错误率下模拟。作为一个实例，图13示出对于使用k＝3和d＝10设定的跨整个大肠杆菌基因组的读段的总相似性评分(S[m])。通过在图14的图中的峰指示的高相似性评分S[m]可检测靠近取样位置。图14说明取样位置(以2,720,230为中心x轴)附近的评分，示出窗口大小选择的平衡点。用于较短窗口大小(图14a、14b，w＝100个碱基)，峰变得嘈杂，并且对于较大窗口大小(图14e、14f，w＝1000个碱基)，峰变得更宽，两者降低检测恰当起始位置的准确度。

综上，数据示出可优化本发明所公开的算法以利用在100-1000个碱基之间的窗口大小，以便甚至在存在相对高错误率(包含插入缺失)的情况下准确和高效地识别在两个序列之间的相似性。这些数据还示出，余弦相似性评分为准确度量并且与其它距离公式相比更高效，因为它可使用快速傅里叶变换(FFT)高效实施。

实验代码以C编程语言实施，以评估在GPU设置上的算法。在使用中的GPU设置为Nvidia Titan X(12G GDDR5)。一般来说，对于长度l的序列，计算归一化的k-聚体计数向量及时采用O(l)，并且每个k-聚体的计算其快速傅里叶变换采用O(l_DFTlog(l_DFT))和总(ql_DFTlog(l_DFT))。计算包含乘以参考和查询FFT向量的总评分S和计算其逆FFT在O(l_DFT(q+log(l_DFT))中，并且对于带状动态编程(BDP)的计算时间在(O(l_读取band_len))中(其中band_len为O(l_读取abs(插入率-缺失率))。可通过将大参考序列拆分成理想转换大小N的短片段和使用重叠保留(或重叠相加)技术(Oppenheim等人，2009)高效计算FFT和IFFT步骤。FFT和BDP运算使用NVIDIA cuFFT和NVBIO文库实施。

实例2

使用大肠杆菌基因组的准确度和性能分析

使用平均长度为5kbps和10kbps和不同序列精度为85％、75％、65％和55％的来自大肠杆菌基因组的20x模拟读取数据集评估此方法的准确度和性能。读取序列使用具有选项(--数据-类型CLR--深度20--model_qc model_qc_clr--准确度-min 0.5--长度-平均[5000|10000]--长度-sd 2000--准确度-平均[.85|.75|.65|.55]--准确度-sd 0.02)的PBSIM(Ono等人，2013)模拟。

在默认设置(w＝500，L_t＝7500，f＝2，g＝1，max-num-顶-峰＝10,max-fft-块-大小＝32768下对于平均序列长度为5kbps和10kbps的数据集分别在表1和2中报告k＝3、4的性能。甚至在～45％错误率的情况下，k＝4具有几乎完美的准确度。如从表2预期，特别在定位覆盖长重复区域的读段时，较长读段导致总比对率较高。如果l_读取<w(考虑在模拟数据集中测序列长度的分布这很少出现)，那么读取标记为跳过。表3还报告执行～45,000个模拟读段(平均5Kbps长)与人类染色体1比对。

表1：在不同错误率的情况下来自大肠杆菌基因组的20X模拟数据集的比对准确度。平均读取序列长度为5kbps。

表2：在不同错误率的情况下来自大肠杆菌基因组的20X模拟数据集的比对准确度。平均读取序列长度为10kbps。

表3：在不同错误率的情况下来自人类染色体1的1X模拟数据集的比对准确度。平均读取序列长度为5kbps。

TGS读段达到几十kpbs并且它们大部分准确度>70％。然而用较短片段(多kbps长)实现高敏感性在对于如组装的应用的原始读段的成对比对中变得更重要，其中读段部分重叠并且错误率为原始读段的2×。存在开发的使用短第二代测序(SGS)序列修正TGS读段的混合方法。但是这些方法需要多个文库准备和测序运行，此外它们引起解析短读段与源自重复区域的嘈杂长读段的复杂度。如所解释，本文公开的方法可用于执行读取序列的成对比对以识别重叠的读取序列。

序列表

<110> 斯坦福大学托管董事会

<120> 第三代测序比对算法

<130> STAN-1259WO

<150> 62/294,205

<151> 2016-02-11

<160> 8

<170> PatentIn version 3.5

<210> 1

<211> 16

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<400> 1

atgcttcgaa tcgata 16

<210> 2

<211> 16

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<400> 2

atccttcgga tcgtta 16

<210> 3

<211> 15

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<400> 3

attattagtc ctagg 15

<210> 4

<211> 15

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<400> 4

attattagat cctag 15

<210> 5

<211> 21

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<400> 5

atgcttcgaa tcgatagcct a 21

<210> 6

<211> 21

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<400> 6

atccttcgga tcgttagtct a 21

<210> 7

<211> 31

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<220>

<221> Misc_Feature

<222> (1)..(8)

<223> n is A T, C, or G

<220>

<221> Misc_Feature

<222> (24)..(31)

<223> n is A T, C, or G

<400> 7

nnnnnnnnaa aagccccggt tttnnnnnnn n 31

<210> 8

<211> 31

<212> DNA

<213> Artificial Sequence

<220>

<223> Synthetic sequence

<220>

<221> Misc_Feature

<222> (1)..(8)

<223> n is A, T, C, or G

<220>

<221> Misc_Feature

<222> (24)..(31)

<223> n is A, T, C, or G

<400> 8

nnnnnnnnaa aaccccggat tttnnnnnnn n 31

Claims

1.一种用于将读取序列与参考序列片段比对的方法，所述方法包括：

a.创建用于所述读取序列的窗口和用于所述参考序列片段的窗口，其中所述窗口具有相同长度；

b.计算在每个窗口内独特k-聚体出现的数量，

c.基于在每个窗口内独特k-聚体出现的所述数量，计算k-聚体计数相似性值；

d.对于跨所述读取序列的多个窗口和跨所述参考序列片段的多个窗口迭代地执行步骤(a)到(c)，由此计算多个k-聚体计数相似性值，

其中在所述读取序列和所述参考序列片段中的每一个中的每个随后窗口的开始从在相应序列中的先前窗口的开始偏移距离d；

e.通过求所述多个k-聚体计数相似性值的平均值，计算相似性评分；和

f.当所述相似性评分高于阈值时，将所述读取序列与所述参考序列片段比对，

其中将在步骤(a)的第一次执行中创建的所述窗口放置在每个序列的起始处。

2.根据权利要求1所述的方法，其中对于所述读取序列和所述参考序列的不同片段重复步骤(a)到(f)。

3.根据权利要求1或2中任一项所述的方法，其中所述参考序列片段为从基因组数据库获得的参考序列的区域。

4.根据权利要求1或2中任一项所述的方法，其中所述参考序列为读取序列。

5.根据权利要求4所述的方法，其中为读取序列的所述参考序列是根据测序同一样品而获得，根据权利要求1所述的读取序列的所述序列是从所述同一样品获得。

6.根据权利要求1到5中任一项所述的方法，其中所述窗口中的每一个的所述长度为至少50个碱基。

7.根据权利要求1到5中任一项所述的方法，其中所述窗口中的每一个的所述长度能够为在1-10,000个碱基范围内的任何整数值，其中所述长度保持恒定。

8.根据权利要求1到7中任一项所述的方法，其中所述距离d为至少10个碱基长。

9.根据权利要求1到7中任一项所述的方法，其中所述距离d的长度能够在1-500个碱基范围内，其中d保持恒定。

10.根据权利要求1到9中任一项所述的方法，其中所述k-聚体的长度为2-10个碱基。

11.根据权利要求10所述的方法，其中所述k-聚体的长度为3个碱基。

12.根据权利要求10所述的方法，其中所述k-聚体的长度为4个碱基。

13.一种可实行软件产品，其存储在含有用于将读取序列与参考序列片段比对的方法的程序指令的计算机可读媒体上，所述方法包括：

b.计算在每个窗口内独特k-聚体出现的数量，

14.根据权利要求13所述的可实行软件产品，其中对于所述读取序列和所述参考序列的不同片段重复步骤(a)到(f)。

15.根据权利要求13或14中任一项所述的可实行软件产品，其中所述参考序列片段为从基因组数据库获得的参考序列的区域。

16.根据权利要求13或14中任一项所述的可实行软件产品，其中所述参考序列为读取序列。

17.根据权利要求16所述的可实行软件产品，其中为读取序列的所述参考序列是根据测序同一样品而获得，根据权利要求13所述的读取序列的所述序列是从所述同一样品获得。

18.根据权利要求13到17中任一项所述的可实行软件产品，其中所述窗口中的每一个的所述长度为至少50个碱基。

19.根据权利要求13到17中任一项所述的可实行软件产品，其中所述窗口中的每一个的所述长度能够为在1-10,000个碱基范围内的任何整数值，其中所述长度保持恒定。

20.根据权利要求13到19中任一项所述的可实行软件产品，其中所述距离d为至少10个碱基长。

21.根据权利要求13到19中任一项所述的可实行软件产品，其中所述距离d的长度能够在1-500个碱基范围内，其中d保持恒定。

22.根据权利要求13到21中任一项所述的可实行软件产品，其中所述k-聚体的长度为2-10个碱基。

23.根据权利要求22所述的可实行软件产品，其中所述k-聚体的长度为3个碱基。

24.根据权利要求22所述的可实行软件产品，其中所述k-聚体的长度为4个碱基。

25.一种用于将读取序列与参考序列片段比对的系统，其包括：

存储器；和

处理器，所述处理器耦合到所述存储器并且被配置成实行存储在所述存储器中的指令，所述指令包括用于以下的指令：

b.计算在每个窗口内独特k-聚体出现的数量，

26.根据权利要求25所述的系统，其中对于所述读取序列和所述参考序列的不同片段重复步骤(a)到(f)。

27.根据权利要求25或26中任一项所述的系统，其中所述参考序列片段为从基因组数据库获得的参考序列的区域。

28.根据权利要求25或26中任一项所述的系统，其中所述参考序列为读取序列。

29.根据权利要求28所述的系统，其中为读取序列的所述参考序列是根据测序同一样品而获得，根据权利要求25所述的读取序列的所述序列是从所述同一样品获得。

30.根据权利要求25到29中任一项所述的系统，其中所述窗口中的每一个的所述长度为至少50个碱基。

31.根据权利要求25到29中任一项所述的系统，其中所述窗口中的每一个的所述长度能够为在1-10,000个碱基范围内的任何整数值，其中所述长度保持恒定。

32.根据权利要求25到31中任一项所述的系统，其中所述距离d为至少10个碱基长。

33.根据权利要求25到31中任一项所述的系统，其中所述距离d的长度能够在1-500个碱基范围内，其中d保持恒定。

34.根据权利要求25到33中任一项所述的系统，其中所述k-聚体的长度为2-10个碱基。

35.根据权利要求34所述的系统，其中所述k-聚体的长度为3个碱基。

36.根据权利要求34所述的系统，其中所述k-聚体的长度为4个碱基。

37.一种存储装置，其存储能够实行以执行包括以下的操作的指令：

a.创建用于读取序列的窗口和用于参考序列片段的窗口，其中所述窗口具有相同长度；

b.计算在每个窗口内独特k-聚体出现的数量，

38.根据权利要求37所述的存储装置，其中对于所述读取序列和所述参考序列的不同片段重复步骤(a)到(f)。

39.根据权利要求37或38中任一项所述的存储装置，其中所述参考序列片段为从基因组数据库获得的参考序列的区域。

40.根据权利要求37或38中任一项所述的存储装置，其中所述参考序列为读取序列。

41.根据权利要求40所述的存储装置，其中为读取序列的所述参考序列是根据测序同一样品而获得，根据权利要求37所述的读取序列的所述序列是从所述同一样品获得。

42.根据权利要求37到41中任一项所述的存储装置，其中所述窗口中的每一个的所述长度为至少50个碱基。

43.根据权利要求37到41中任一项所述的存储装置，其中所述窗口中的每一个的所述长度能够为在1-10,000个碱基范围内的任何整数值，其中所述长度保持恒定。

44.根据权利要求37到43中任一项所述的存储装置，其中所述距离d为至少10个碱基长。

45.根据权利要求37到43中任一项所述的存储装置，其中所述距离d的长度能够在1-500个碱基范围内，其中d保持恒定。

46.根据权利要求37到45中任一项所述的存储装置，其中所述k-聚体的长度为2-10个碱基。

47.根据权利要求46所述的存储装置，其中所述k-聚体的长度为3个碱基。

48.根据权利要求46所述的存储装置，其中所述k-聚体的长度为4个碱基。