CN105637098B

CN105637098B - 用于比对序列的方法和系统

Info

Publication number: CN105637098B
Application number: CN201480055430.1A
Authority: CN
Inventors: 丹尼斯·库拉尔
Original assignee: Seven Bridges Genomics Inc
Current assignee: Seven Bridges Genomics Inc
Priority date: 2013-08-21
Filing date: 2014-08-21
Publication date: 2021-04-30
Anticipated expiration: 2034-08-21
Also published as: JP6902073B2; CA2921645C; AU2014308794A1; EP3036344A1; CN105637098A; KR20160047506A; US11211146B2; EP3702472A1; JP2020042813A; US20150057946A1; JP6586088B2; KR102371953B1; US20180357367A1; JP2016536698A; WO2015027050A1; EP3036344A4; US9898575B2; CA2921645A1; EP3036344B1; EP3702472B1

Abstract

本发明包含用于比对读数(例如，核酸读数、氨基酸读数)与参考序列构建物的方法、用于构建该参考序列构建物的方法、以及使用该比对方法和该构建物以产生序列的系统。该方法是可扩展的，并且可以用来比对数百万读数与数千碱基或氨基酸长度的构建物。本发明另外包含用于基于核酸读数与该构建物中的位置的比对来识别疾病或基因型的方法。

Description

用于比对序列的方法和系统

相关申请案

本申请案主张2013年9月3日提交的第14/016,833号和2013年8月21日提交的第61/868,249号美国专利申请案的优先权，所述两个申请案的全部内容以引用的方式并入本文中。

技术领域

本发明涉及用于对序列(例如，核酸序列、氨基酸序列)彼此比对以产生对应于样本(例如，遗传样本、蛋白质样本)的连续序列读数的方法和系统。本发明另外涉及用于识别样本中的变异的方法。

背景技术

遗传学已经从分析科学演变为信息科学。然而，科学家此前一直努力研究如何提取和识别核酸，此类技术现在看来并非那么重要。下一代测序(例如，全转录组鸟枪测序、焦磷酸测序、离子半导体测序、使用合成法测序)可以在仅几天内产生覆盖全基因组的数百万读数。为了实现此产出量，NGS测序在较小核酸序列上使用大规模并行计算，其一起组成大量遗传信息，例如，染色体或基因组。从遗传样本开始，核酸(例如，DNA)被分裂、扩增、并以极快速度读取。考虑到这些能力，科学家现在努力研究如何(以低成本)比对读数以识别序列中指示疾病或疾病风险的基因座。

当前技术发展水平的比对方法使用大量计算能力来比对重叠读数与参考以产生可探测用于重要遗传信息或结构信息(例如，用于疾病的生物标志物)的序列。最终，序列比对的目标是组合由定序器产生的核酸读数集以实现较长读数(即，重叠群)或甚至基于来自受试者的遗传样本的该受试者的全基因组。因为来自下一代定序器的序列数据通常包括一起表示目标序列的总数的数百万短序列，所以比对读数复杂且在计算上昂贵。另外，为了使由随机测序误差(即，不正确的测序仪输出)引起的序列失真减到最少，对探测的序列的每个部分多次(例如，2次到100次或更多)测序，以使任何随机测序误差对所产生的最终比对和输出序列的影响减到最少。最后，一旦收集了对应于所有核酸读数的所有数据，就比对该读数与单个参考序列(例如，GRCh37)，以便确定所有(或一部分)受试者序列。在许多情况下，实际上不显示个别读数，而是将比对序列组装为一个序列，并作为数据文件提供该序列。

通常，通过聚集序列信息的两个线性字符串之间的成对比对来构建序列比对。作为比对的实例，可以将两个字符串S1(序列编号12：AGCTACGTACACTACC)和S2(序列编号13：AGCTATCGTACTAGC)与彼此进行比对。S1通常对应于读数，而S2对应于参考序列的一部分。S1和S2可关于彼此通过替代、删除和插入构成。通常，相对于将字符串S1转换为字符串S2来定义这些术语：当用S1中相同长度的不同字母或序列替代S2中的字母或序列时发生替代，当在S1的对应区段中“跳过”S2中的字母或序列时发生删除，并且当在S1中的两个位置(这两个位置在S2中为相邻位置)之间出现字母或序列时发生插入。例如，可以对两个序列S1和S2比对如下。以下比对指出有十三处匹配，一处删除长度一，一处插入长度二以及一处替代：

(S1)AGCTA-CGTACACTACC(序列编号12)

(S2)AGCTATCGTAC--TAGC(序列编号13)

本领域的技术人员将了解，存在序列比对的精确算法和近似算法。精确算法将找出最高得分的比对，但是在计算上会昂贵。两个最著名的精确算法是尼德曼-翁施(Needleman-Wunsch)算法(分子生物学杂志(J Mol Biol)48(3)：443-453，1970)和史密斯-沃特曼(Smith-Waterman)算法(分子生物学杂志(J Mol Biol)147(1)：195-197，1981；数学进展(Adv.in Math.)20(3)，367-387，1976)。后藤(Gotoh)对史密斯-沃特曼算法的进一步改进(分子生物学杂志(JMol Biol)162(3)，705-708，1982)减少了从O(m²n)到O(mn)的计算时间，其中m和n比较的序列大小，该改进更能改善并行处理。在生物信息学领域，正是后藤的改良算法通常被称为史密斯-沃特曼算法。史密斯-沃特曼方法用于比对较大序列集与较大参考序列，因为可更普遍且更便宜地获得并行计算资源。参考例如，在http://aws.amazon.com可获得的Amazon.com的云计算资源。所有上述期刊论文的全部内容以引入的方式并入本文中。

史密斯-沃特曼(SW)算法通过奖励序列中的碱基之间的重叠并处罚序列之间的空位来比对线性序列。史密斯-沃特曼算法还与尼德曼-翁施算法不同，不同之处在于SW不要求短序列跨越描述长序列的字母组成的字符串。也就是说，SW不假定一个序列是另一个序列的全部内容的读数。此外，因为SW并不一定找出横跨字符串的全长的比对，所以局部比对可以在两个序列内的任何地方开始和结束。

根据以下方程式(1)，对于表示长度n和m的两个字符串的n×m矩阵H，易于表示SW算法：

H_ko＝H_ol＝0(对于0≤k≤n且0≤l≤m) (1)

H_ij＝max{H_i-1,j-1+s(a_i,b_j),H_i-1,j-W_in,H_i,j-1-W_del,0}

(对于1≤i≤n且1≤j≤m)

在以上方程式中，s(a_i,b_j)表示匹配奖励值(当a_i＝b_j时)或不匹配罚分(当a_i≠b_j时)，并且对插入和删除分别给出罚分W_in和W_del。在大多数例子中，所得矩阵具有许多为零的单元。这种表示使得更容易在矩阵中从高到低、从右到左回溯，因此识别比对。

一旦已经用得分完全填充矩阵，SW算法就执行回溯以确定比对。开始于矩阵中的最大值，算法将基于三个值中的哪个(H_i-1,j-1、H_i-1,j或H_i,j-1)曾用来计算每个细胞的最终最大值来进行回溯。当达到零时回溯停止。见例如图3(B)，其不表示现有技术，而是示出回溯的概念以及在读取回溯时的对应局部比对。因此，如通过算法确定的“最佳比对”可以含有超过最小可能数目的插入和删除，但是将含有远少于最大可能数目的替代。

当作为SW或SW-后藤应用时，该技术使用动态规划算法来执行分别具有大小m和n的两个字符串S和A的局部序列比对。此动态规划技术采用表或矩阵来保存匹配得分并避免对于连续细胞的重新计算。可以相对于序列的字母为字符串的每个单元编索引，也就是说，如果S是字符串ATCGAA，那么S[1]＝A、S[4]＝G等。替代将最优比对表示为H_i,j(上文)，可以将最优比对表示为以下方程式(2)中的B[j,k]：

B[j,k]＝max(p[j,k],i[j,k],d[j,k],0)(对于0<j≤m、0<k≤n) (2)

在以下方程式(3)到(5)中概述最大值函数B[j,k]的变量参数，其中MISMATCH_PENALTY、MATCH_BONUS、INSERTION_PENALTY、DELETION_PENALTY和OPENING_PENALTY都是常数，并且除MATCH_BONUS以外均为负数。匹配变量参数p[j,k]由以下方程式(3)得出：

若S[j]≠A[k]，则p[j,k]＝max(p[j-1,k-1],i[j-1,k-1],d[j-1,k-1])+MISMATCH_PENALTY

若S[j]＝A[k]，则p[j,k]＝＝max(p[j-1,k-1],i[j-1,k-1],d[j-1,k-1])+MATCH_BONUS (3)

插入变量参数i[j,k]由以下方程式(4)得出：

i[j,k]＝max(p[j-1,k]+OPENING_PENALTY,i[j-1,k],d[j-1,k]+OPENING_PENALTY)+INSERTION_PENALTY (4)

且删除变量参数d[j,k]由以下方程式(5)得出：

d[j,k]＝max(p[j,k-1]+OPENING_PENALTY,i[j,k-1]+OPENING_PENALTY,d[j,k-1])+DELETION_PENALTY (5)

对于所有三个变量参数，将[0,0]单元设置为零以确保回溯完成，即，p[0,0]＝i[0,0]＝d[0,0]＝0。

得分参数在一定程度上是任意的，并可经调整以实现计算的性能。对于DNA的得分参数设置的一个实例(Huang，第3章：生物序列比较和比对(Bio-Sequence Comparison andAlignment)，Curr Top Comp Mol Biol.丛书，马萨诸塞州剑桥市：麻省理工学院出版社(The MIT Press)，2002年)将为：

MATCH_BONUS：10

MISMATCH_PENALTY：-20

INSERTION_PENALTY：-40

PENING_PENALTY：-10

DELETION_PENALTY：-5

以上空位罚分(INSERTION_PENALTY、OPENING_PENALTY)之间的关系有助于限制空位开放的数目，即，支持通过设置高于空位开放成本的空位插入罚分来归并空位。当然，MISMATCH_PENALTY、MATCH_BONUS、INSERTION_PENALTY、OPENING_PENALTY和DELETION_PENALTY之间可能存在替代关系。

一旦完成比对，可以组装比对后的序列以产生可与参考(即，遗传标准)相比以识别变异的序列。变异可以提供关于疾病、疾病期、复发等的洞察。在氨基酸比对的情况下，可以比较组装后的氨基酸序列与标准以确定关于蛋白质的进化信息或关于蛋白质的功能信息。然而，疾病比较的此标准方法是费时的，因为许多变异不一定与疾病相关。例如，当遗传标准来自具有与样本不同的血统的人群时，许多所谓的变异是归因于像毛色、肤色等的事物的差别。

发明内容

本发明提供算法和方法，该算法和方法的实施将线性的局部序列比对方法(例如，史密斯-沃特曼-后藤方法)转换成多维比对算法，该算法和方法提高并行计算、提高速度、提高精确度并且能够贯穿全基因组比对读数。本发明的算法提供(如史密斯-沃特曼算法中的)序列信息的“回顾”类型分析，然而，与已知线性方法对比，本发明的回顾贯穿包含多个通路和多个节点的多维空间而进行，以便提供对复杂和冗长序列读数的更精确比对，同时实现更低的总错配率、删除率和插入率。

在实践中，通过比对序列读数与跨越分歧点的一系列有向非循环序列来实施本发明，其考虑比对中的所有或几乎所有可能的序列变异，包括插入、删除和替代。通常表示为有向非循环图(DAG)的此类构建物可易于从可用的序列数据库来组装，可用的序列数据库包含“接受”参考序列和变异识别格式(VCF)条目。当结合DAG或其它有向构建物时，所公开的算法因此提供针对序列比对的多维方法，其大大改进比对精确度并提供通过传统算法不可能实现的序列分辨率。该技术可以与任何序列信息一起使用，然而实际上，如本文中所论述，该技术最适用于比对核酸序列和氨基酸序列。

本发明另外提供使用参考序列构建物(例如，表示基因组的每个基因座处的已知变异的DAG)在特定基因座处进行特定碱基响应的方法。因为在比对期间对序列读数与DAG进行比对，所以可以排除比较关于参考基因组的突变与已知突变的表的后续步骤。使用所公开的方法，需要做的仅仅是识别如位于DAG上表示的已知突变处的核酸读数并且响应该突变。替代地，当突变并非已知(即，未在参考序列构建物中表示)时，将找出比对并且将变异识别为新突变。该方法还使得有可能将例如特定疾病风险或疾病进展的另外的信息与并入到参考序列构建物中的已知突变相关联。此外，除了能够在比对期间找出所有基因相关的结果之外，所公开的方法还减少进行比对所需的计算资源，同时允许与多个参考序列的同步比较。

本发明另外包含用于构建表示生物体的序列内的位置处的已知变异的有向非循环图数据结构(DAG)的方法。DAG可以包含数千个位置处的多个序列，并且可以包含每个位置处的多个变异，包含删除、插入、平移、倒置和单核苷酸多态性(SNP)。还可能给DAG中的每个变异标记相关诊断信息，例如“乳腺癌”，由此减少识别针对提供样本的患者的风险所需的步骤。在一些实施例中，将对变异评分、加权或使其与其它变异相关以反映该变异作为疾病标志的发生率。

本发明另外包含用于执行本发明的方法的系统。在一个实施例中，系统包括处理器和存储装置的分布式网络，能够比较多个序列(即，核酸序列、氨基酸序列)与表示基因组或基因组区中观察到的变异的参考序列构建物(例如，DAG)。该系统另外能够使用高效比对算法来比对核酸读数以产生连续序列。因为参考序列构建物压缩大量冗余信息，并且因为比对算法如此高效，所以可以使用市售资源在全基因组上标记和组装该读数。该系统包括多个处理器，多个处理器同时执行多个读数与参考序列构建物之间的多个比较。可以累计比较数据并提供给医疗服务人员。因为该比较在计算上易处理，所以分析序列读数将不再表示在NGS测序与患者遗传风险的有意义的探讨之间形成瓶颈。

附图说明

图1描绘表示参考序列中的遗传变异的有向非循环图(DAG)的构建物。图1(A)示出开始参考序列和删除的添加。图1(B)示出插入和SNP的添加，因此得出用于比对的最终DAG；

图2描绘表示为有向非循环图的三个变异识别格式(VCF)条目；

图3(A)示出将核酸序列读数与解释插入情况的构建物以及参考序列进行比对的图形表示；

图3(B)示出用以识别核酸序列读数“ATCGAA”的适当位置的矩阵和回溯；

图4描绘用于并行处理的关联计算模型；

图5描绘用于并行计算的体系构建物。

具体实施方式

本发明包含用于比对序列(例如，核酸序列、氨基酸序列)与参考序列构建物的方法、用于构建参考序列构建物的方法、以及使用比对方法和构建物以产生比对和组装的系统。参考序列构建物可以是如下文所述的有向非循环图(DAG)，然而，只要是针对比对格式化的构建物，参考序列可以是反映物种内的不同生物体的序列中的遗传变异性的任何表示法。遗传变异性也可以是在生物体内的不同组织或细胞之间。一般来说，参考序列构建物将包括取样序列之间相同的部分和不同的部分。因此，该构建物可以被认为是具有包括相同序列的位置(即，根据某种规范排序)以及包括替代序列的一些位置，从而反映遗传变异性。本申请案另外公开用于基于核酸读数与构建物中的位置的比对来识别疾病或基因型的方法。该方法广泛适用于遗传测序和突变筛查领域。

参考序列构建物

与使用用以比对的单个参考序列和基因型核酸读数的现有技术序列比对方法不同，本发明使用可以考虑物种、人群内的遗传序列中的变异性、或甚至单个生物体中的不同细胞之间的遗传序列中的变异性的构建物。遗传变异的表示可以呈现为有向非循环图(DAG)(上文所论述)、行列比对矩阵或deBruijn图，并且只要比对算法的参数设置适当(下文论述)，这些构建物就可以用于本发明的比对方法。

在本发明的优选实施例中，该构建物是有向非循环图(DAG)，即，具有方向且具有非循环路径。(也就是说，一条序列路径无法多次穿过参考构建物上的一个位置。)在DAG中，序列中的遗传变异表示为替代节点。该节点可以是保守序列的一部分、或基因或简单地是核酸。贯穿构建物的不同可能路径表示已知的遗传变异。DAG可以是为生物体的全基因组而构建，或DAG可以是仅为基因组的一部分(例如，染色体，或更小一段遗传信息)而建构。在一些实施例中，DAG表示多于1000个核酸，例如多于10,000个核酸，例如，多于100,000个核酸，例如，多于1,000,000个核酸。DAG可以表示物种(例如，人类)或所选择的人群(例如，患有乳腺癌的妇女)，或甚至更小的亚群，例如在同一个体的不同肿瘤细胞之间的遗传变异。

图1中示出DAG构建物的简单实例。如图1(A)中所示，DAG以参考序列开始，该参考序列在图1(A)中示出为序列编号1：CATAGTACCTAGGTCTTGGAGCTAGTC。在实践中，参考序列通常长得多，并且可以是全基因组。序列通常存储为FASTA或FASTQ文件。(FASTQ已经成为产生自下一代定序器的序列数据的默认格式)。在一些实施例中，参考序列可以是标准参考，例如GRCh37。如所属领域的技术人员所认识到，序列中的每个字母(或符号)实际上对应于核苷酸(例如，脱氧核糖核苷酸或核糖核苷酸)或氨基酸(例如，组氨酸、亮氨酸、赖氨酸，等)。

在下一步骤，向参考序列添加变异，如图1(A)的底部图像中所示。如图1(A)中所示，该变异是从参考中删除图中的细线之间的序列“AG”，即，序列编号2。此删除以图形方式表示为将参考序列分成在删除之前和之后的节点，并在该节点之间插入两个字符串。节点之间的一个路径表示参考序列，而另一个路径表示删除。

在实践中，通过应用例如可在1000个基因组的专题网站中查询到的变异识别格式(VCF)文件中的条目将变异调用到DAG。因为每个VCF文件针对特定的参考基因组，所以不难识别字符串应位于何处。事实上，VCF文件中的每个条目可以被认为是与对参考组合以形成离散图，如图2中所显示。应注意，图2中的VCF条目不与图1的VCF条目相对应。

移到图1(B)，添加对应于特定位置处的插入“GG”的第二VCF条目以产生扩展的DAG，即，包含序列编号3和序列编号4。接着，可以添加第三VCF条目扩展DAG以在参考序列中较早地说明SNP，即，包含序列编号5到8。因此，在三个步骤中，已经创建DAG，可以针对该DAG比对核酸读数(如下文所论述)。

在实践中，DAG在计算机存储器(硬盘、闪存、云存储器等)中表示为节点集S，其中每个节点由字符串、父节点集和位置界定。字符串是节点的“内容”，即，序列；父节点界定节点相对于图中的其它节点的位置；并且节点的位置与系统中的某种规范排序有关，例如，参考基因组。虽然并非绝对必需相对于参考序列界定该图，但是它确实使输出数据的操作更为简单。当然，对S的进一步约束是它不能包含环路。

将此DAG方法推断至较大构建物，有可能构建这样的DAG：其并入表示参考的指定区的遗传序列中的已知变异的数千VCF条目。然而，由于DAG变得更庞大，因此计算确实花费更长时间，并且对于许多应用使用仅可表示序列的一部分(例如，染色体)的较小DAG。在其它实施例中，可以通过减小由DAG覆盖的人群的规模将DAG制作为更小，例如，从表示乳腺癌中的变异的DAG变为表示三阴性乳腺癌中的变异的DAG。替代地，可以使用基于易于识别的遗传标志物自定义的更长DAG，易于识别的遗传标志物将通常导致DAG的大部分在样本之间保持一致。例如，相比于解释同一序列上已知的人类的所有变异的DAG，针对以来自非洲血统女性的VCF条目创建的DAG，比对来自非洲血统女性的核酸读数集将更快速。还应认识到，本发明的DAG是动态构建物，这是因为它们可随时间推移而改变以并入最新识别的突变。另外，其中将比对结果以递归方式添加到DAG的算法也是可能的。

在字符串与DAG比对的实例中，可以调整空位罚分以使空位插入甚至更昂贵，因此更倾向于对序列的比对而不是在总体序列中开放新的空位。当然，通过DAG的改进(上文所论述)，因为在DAG中考虑了突变，所以应更进一步减少空位的发生。

比对算法

在一个实施例中，使用算法来比对序列读数与有向非循环图(DAG)。与背景技术中表达的算法对比，该比对算法通过识别关于在DAG(例如，参考序列构建物)上的位置处含有的每个序列的最大得分来识别C_i,j的最大值。实际上，通过在先前位置处“向后”看，有可能跨多个可能的路径识别最优比对。

如上文所论述，对读数(也称为“字符串”)和有向非循环图(DAG)执行本发明的算法。出于定义该算法的目的，假设S为要比对的字符串，假设D为将针对S比对的有向非循环图。以从1开始的索引对字符串S的单元加括号。因此，如果S是字符串ATCGAA，那么S[1]＝A、S[4]＝G等。

对于DAG，节点的序列的每个字母将表示为单独的单元d。d的前趋定义为：

(i)如果d不是其节点的序列的首字母，那么其节点中在d之前的字母是其(唯一)前趋；

(ii)如果d是其节点的序列的首字母，那么作为d的节点的父节点的任何节点的序列的最后一个字母是d的前趋。

所有前趋集继而表示为P[d]。

为了找出“最佳”比对，算法寻求M[j,d]的值，即S的第一j个单元与在d之前(且包含d)的DAG的部分的最佳比对的得分。此步骤类似于在背景技术章节的方程式1中找出H_i,j。具体来说，确定M[j,d]包括找出a、i、e和0的最大值，如下文所定义：

M[j,d]＝max{a,i,e,0} (6)

其中

对于P[d]中的p*，e＝max{M[j,p*]+DELETE_PENALTY}

i＝M[j-1,d]+INSERT_PENALTY

若S[j]＝d，则对于P[d]中的p*，a＝max{M[j-1,p*]+MATCH_SCORE}；

若S[j]≠d，则对于P[d]中的p*，a＝max{M[j-1,p*]+MISMATCH_PENALTY}

如上文所描述，e是S的第一j个字符与DAG的直到(但是不包含)d的部分比对的最高值加上另外的DELETE_PENALTY。因此，如果d不是节点的序列的首字母，那么仅存在一个前趋p，并且S的第一j个字符与DAG(直到且包含p)的比对得分相当于M[j,p]+DELETE_PENALTY。在其中d是节点的序列的首字母的实例中，可以存在多个可能的前趋，并且因为DELETE_PENALTY是恒定的，所以求[M[j,p*]+DELETE_PENALTY]的最大值就相当于选择与S的第一j个字符比对而具有最高比对得分的前趋。

在方程式(6)中，i是字符串S的第一j-l个字符与(直到且包含d的)DAG的比对加上INSERT_PENALTY，这类似于SW中的插入变量参数的定义(见方程式1)。

另外，a是S的第一j个字符与DAG的直到(但是不包含)d的部分比对的最高值加上或MATCH_SCORE(如果S的第j个字符与字符d相同)或MISMATCH_PENALTY(如果S的第j个字符与字符d不同)。如同e一样，这意味着如果d不是节点的序列的首字母，那么仅存在一个前趋，即，p。这意味着a是S的第一j-l个字符与DAG(直到且包含p)的比对得分，即，M[j-l,p]，取决于d与S的第j个字符是否匹配，再加上MISMATCH_PENALTY或MATCH_SCORE。在其中d是节点的序列的首字母的实例中，可以存在多个可能的前趋。在此情况下，求{M[j,p*]+MISMATCH_PENALTY或MATCH_SCORE}的最大值就相当于选择与S的第一j-l个字符比对而具有最高比对得分(即，候选变量参数的最高值)并取决于d与S的第j个字符是否匹配而加上MISMATCH_PENALTY或MATCH_SCORE的前趋。

此外，如在背景技术中所论述的SW算法中一样，可以调整罚分(例如DELETE_PENALTY、INSERT_PENALTY、MATCH_SCORE和MISMATCH_PENALTY)以促进与更少空位的比对等。

如以上方程式中所描述，该算法不仅通过计算该单元的插入、删除和匹配得分来找出每个读数的最大值，而且还向后看(逆着DAG的方向)至DAG上的任何先前节点以找出最大得分。因此，该算法能够详细研究贯穿DAG的含有已知突变的不同路径。因为图是有向的，所以逆着图的方向移动的回溯遵循朝向图的起点的优选变异序列，并且最大比对得分识别高度确定性的最可能比对。虽然以上方程式表示为“最大”值，但是“最大”预期涵盖任何形式的优化，包含(例如)转换所有方程式上的正负号并求解最小值。

在图3中例示所公开的算法的实施方案，其中对序列“ATCGAA”与表示序列编号10：TTGGATATGGG和已知插入情况序列编号11：TTGGATCGAATTATGGG的参考序列的DAG进行比对，其中插入是带下划线的。图3(A)示出与DAG相比的读数的图形表示，而图3(B)示出与该比较相对应的实际矩阵。类似于背景技术中所论述的史密斯-沃特曼技术，本发明的算法识别最高得分并且执行回溯以识别读数的适当位置。图3(A)和(B)还突出本发明产生字符串与构建物的实际匹配，而已知方法(例如，SW)将已更可能对字符串与参考的错误部分进行比对，或由于未产生包含于比对中的足够高的比对得分而拒绝该字符串。在其中序列读数包含未包含在DAG中的变异的实例中，将通过空位、插入等报告已比对序列。

并行计算的可能性

已经针对大规模并行计算调适并明显地修改了史密斯-沃特曼-后藤算法的序列版本。例如，美国专利公开案第2012/0239706号中描述了使用关联大规模并行计算(SWAMP)的被称作史密斯-沃特曼的ASC模型，该公开案以全文引用的方式并入本文中。SWAMP(和其它并行处理系统)的并行计算的部分源于沿任何反对角的值彼此独立的事实。因此，可以并行完成沿指定反对角的所有细胞以分配计算资源。以上递归方程中示出的数据相依性限制了可达到的并行性的水平，但是使用波前方法仍将加速此适用算法。沃兹尼亚克(Wozniak)对Sun Ultra SPARC实施的波前方法(生物科学中的计算应用(CABIOS)，13(2)：145-150，1997)使用专用SIMD类视频指令。沃兹尼亚克使用了SIMD寄存器来存储与次对角线并行的值，从而报告了优于相同机器上的传统实施方案的双倍加速。在沃兹尼亚克的实例之后，并行化码的类似方式将使用流SIMD扩展(SSE)设置用于x86体系构建物。由英特尔(Intel)设计，向量类操作一次对少量值(通常四个、八个或十六个)完成单一操作/指令。许多AMD和英特尔芯片支持各种型式的SSE，并且英特尔已经继续使用用于其调制解调器芯片组的高级矢量扩展(AVX)来开发此技术。

在其它实施方案中，罗盖(Rognes)和塞贝格(Seeberg)(生物信息学(Bioinformatics)(英格兰牛津)16(8)：699-706，2000)使用具有SSE的旧版-MMX SIMD指令的英特尔奔腾处理器用于其实施方案。出于罗盖和塞贝格的研究而针对并行比对研发的方法(生物信息学(Bioinformatics)16(8)：699-706，2000)不使用波前方法(罗盖，核酸研究(Nuc Acids Res)，29(7)：1647-52，2001；塞伯(Saebo)等人，核酸研究(Nuc Acids Res)，33(增刊2)：W535-W539，2005)。替代地，他们比对与查询序列并行的SIMD寄存器，从而使用预先计算出的查询专用得分矩阵一次计算八个值。此方法的另外的细节可以在U.S.7,917,302中找到，其以引用的方式并入本文中。罗盖和塞贝格布局SIMD寄存器的方式-北邻域相依性会去除达三分之一的由SSE并行“向量”计算获得的潜在加速。为了克服这点，它们结合了SWAT类优化。使用较大仿射空位罚分，北邻域大部分时间将为零。如果这是正确的，那么程序可以跳过计算北邻域的值，法勒(Farrar)称其为“惰性F评估”(生物信息学(Bioinformatics)，23(2)：156-161，2007)。罗盖和塞贝格能够减少方程式1的计算数目以通过在该值低于一定阈值时跳过该值来加速他们的算法。在(罗盖和塞贝格，生物信息学(Bioinformatics)16(8)：699-706，2000)中经由MMX/SSE指令和SWAT类扩展使用8路矢量报告了六倍加速。

在法勒(Farrar)完成的SSE研究(生物信息学(Bioinformatics)23(2)：156-161，2007)中，使用线条访问模式或交叉访问模式以排列与查询寄存器并行的SIMD寄存器。这样做避免了任何重叠相依性。此外，结合SWAT类优化(法勒，生物信息学(Bioinformatics)23(2)：156-161，2007)实现的加速是沃兹尼亚克(CABIOS 13(2)：145-150，1997)以及罗盖和塞贝格(生物信息学(Bioinformatics)(英格兰牛津)16(8)：699-706，2000)实施方案的2到8倍。块替代矩阵和高效且有独创性的内循环(具有移到内循环外的北(F)条件性)是重要的优化。用于处理的十六个8位单元的交叉存储器模式访问也改进了存储器访问时间，从而促进总体加速。

法勒(序列分析(Sequence Analysis)，2008)扩展了其对索尼(Sony)、东芝(Toshiba)和IBM制造的Cell处理器的研究。这种Cell处理器具有一个主核心和八个小核心。Cell宽带引擎是若干更多史密斯-沃特曼实施方案的开发平台，该实施方案包含斯卡科夫斯基(Szalkowski)等人的SWPS3(BMC研究笔记(BMC Res Notes)1(107)，2008)和韦拉万(Wirawan)等人的CBESW(BMC生物信息学(BMC Bioinformatics)9(377)2008)，均使用法勒的分段方法。鲁德尼茨基(Rudnicki)等人(基金信息(Fund Inform.)96，181-194，2009)使用了PS3以开发使用了对多个数据库序列进行并行计算的方法。

罗盖(BMC生物信息学(BMC Bioinformatics)12(221)，2011)还研发了并行处理多个数据库序列的被称为SWIPE的多线程方法。重点是对“普通CPU”使用SIMD方法。并行使用多个数据库序列来使用粗粒并行性分解研究的此调查类似于刘(Liu)等人(BMC研究笔记(BMC Res Notes)2(73)，2009)以及罗戈夫斯基(Ligowski)和鲁德尼茨基(高性能计算生物学八周年国际研讨会(Eight Annual International Workshop on High PerformanceComputational Biology)，罗马，2009)在CUDASW中描述的基于图形处理器单元(GPU)的工具。已经存在使用刘等人(BMC研究笔记(BMC Res Notes)3(93)，2010)和罗戈夫斯基等人(GPU计算，Gems Emerald版(GPU Computing Gems，Emerald Edition)，摩尔根考夫曼(Morgan考夫曼)，155-157，2011)的CUDASW++2.0进行GPU研究的其它实施方案。

在其它变化形式中，可以使用小规模的向量并行计算(8路、16路或32路并行性)以经由并行比对多个序列的GPU实施方案实现可用计算。对于计算的理论峰值加速是m的因素，这是最佳的。使用96个处理单元的ClearSpeed实施方案的96倍加速证实了该理论加速。

并行计算模型

用来开发和扩展史密斯-沃特曼序列比对的主要并行模型是关联计算(ASC)(波特(Potter)等人，计算机(Computer)，27(11)：19-25，1994)。本文中描述了史密斯-沃特曼算法的高效并行型式。本章节中详细描述此模型和另一个模型。

这里定义一些相关词汇。来自计算机体系构建物的费林分类法的两个关注的术语是MIMD和SIMD这两种不同模型的并行计算。将归类为多指令多数据(MIMD)模型的计算机集群用作概念证明以克服超大型比对中的存储器限制。第8节了描述MIMD模型的使用。还描述了被称为ASC的扩展的数据并行的、单指令多数据(SIMD)模型。

多指令多数据(MIMD)

多数据多指令模型或MIMD模型描述目前可用的大多数并行系统，并且包含目前流行的计算机集群。MIMD处理器具有完备的中央处理单元(CPU)，每个CPU具有其自身的本地存储器(奎因(Quinn)，并行计算(Parallel Computing)：理论与实践(Theory andPractice)，第二版，纽约：麦格劳-希尔出版社(McGraw-Hill)，1994)。与SIMD模型对比，MIMD处理器中的每一个存储并异步地执行其自身的程序。MIMD处理器经由允许其通信的网络连接，但是所使用的网络可大不相同，范围为以太网、Myrinet和机器(集群节点)之间的InfiniBand连接。通信往往采用比SIMD宽松得多的通信构建物，超出单个单元。沿着该网络通过个别处理器在个别处理器所执行的个别程序的控制下异步地移动数据。通常，由支持消息传递的若干不同的并行语言之一来处理通信。用于此的非常普遍的库被称为消息传递接口(MPI)。“SIMD类”型式的通信是可能的，但是数据的移动将是异步的。通过MIMD的并行计算通常需要大量的通信和频繁的同步，除非由处理器执行的各种任务是高度独立的(即，所谓的“易并行”或“乐于并行”的问题)。第8节中提出的研究使用经由InfiniBand连接的AMD Opteron集群。

与SIMD不同，消息传递所需的最坏情况下的时间难以预测或不可能预测。通常，使用往往由试验确定的平均情况下的估计来确定MIMD软件的消息传递执行时间，而不是通过对于SIMD来说典型的最坏情况下的理论评估来确定。由于对于MIMD软件最坏的情况往往是非常糟糕的且很少发生，因此平均情况下的估计更加有用。因此，MIMD在特定问题上所需的通信时间会比SIMD更多且通常比SIMD明显更多。这使得MIMD编程(尤其是在使用了消息传递时)的重要目标是使所需的处理器间通信的数目减到最少并使处理器通信之间的时间量增到最大。即使在单卡加速水平下(例如，使用图形处理器或GPU)也是如此。

数据并行编程也是MIMD编程的重要技术，但在这里，所有任务对不同的数据执行相同的操作，并且仅在各个关键点同步。MIMO系统的大部分算法以单程序多数据(SPMD)编程范例编写。每个处理器有其自身的相同程序的拷贝，执行特定于该处理器或其本地数据上的核心的代码部分。SPMD范例的普及源于这一事实：写入大量将同时跨不同处理器执行但仍然能够协作解决单个问题的不同程序是相当困难的。针对存储器密集型而不是计算密集型的问题所使用的另一种方法是使用在第8节中提出的研究创建如使用JumboMem完成的虚拟存储器服务器。这将在其底层实现中使用MPI。

单指令多数据(SIMD)

SIMD模型由被称作PE的多个简单的算术处理单元构成。每个PE具有其自身的本地存储器，PE可以自该本地存储器提取和存储程序，但是PE并没有编译或执行程序的能力。如本文所使用，术语“并行存储器”是指在计算系统中统称的本地存储器。例如，并行存储器可以是SIMD计算机系统的本地存储器(例如，PE的本地存储器)的集合、MIMD计算机系统中的处理器的本地存储器(例如，中央处理单元的本地存储器)的集合等。通过被称作控制单元(或前端)的处理器来处理程序的编译和执行(奎因(Quinn)，并行计算(ParallelComputing)：理论与实践(Theory and Practice)，第二版，纽约：麦格劳-希尔出版社(McGraw-Hill)，1994)。控制单元通常通过总线连接到所有PE。

所有有效PE步调一致地同步执行从控制单元接收到的程序指令。在第79页，“在任何时间单元中，单一操作与各自操作不同数据的多个处理单元上的执行状态相同”(奎因(Quinn)，并行计算(Parallel Computing)：理论与实践(Theory and Practice)，第二版，纽约：麦格劳-希尔出版社(McGraw-Hill)，1994)。虽然通过所有有效PE并行地同时执行相同指令，但是可允许一些PE跳过任何特定指令(贝克(Baker)，SIMD和MASC：CS 6/73301课程笔记：并行和分布式计算-PowerPoint幻灯片，(2004)2004)。这通常使用“如果-否则”分支结构来实现，其中一些PE执行如果指令，而其余的PE执行否则部分。这种模型对于具有可同时出现的最多少量如果-否则分支结构(例如，图像处理和矩阵运算)的本质上“数据并行的”问题是理想的。

控制单元可以将数据广播到所有有效PE，并且控制单元还可以从使用控制单元与PE之间的连接(通常是总线)的特定PE获取数据值。另外，该组PE通过例如线性阵列、2D网格或超立方体的互连网络连接，互连网络提供PE之间的并行数据移动。通过PE以同步并行型式贯穿此网络移动数据，该PE步调一致地执行包含数据移动的指令。正是控制单元将指令广播到PE。具体来说，如今SIMD网络不使用由多数并行计算机所使用的消息传递范例。如此的一个重要优点是，SIMD网络通信极其高效，并且可以通过控制该特定通信的算法的最坏情况下的时间来确定用于通信所需的最大时间。

本章节的剩余部分专门描述扩展的SIMD ASC模型。对于此论述，ASC是算法设计和开发的中心内容。

关联计算模型

关联计算(ASC)模型是扩展的SIMD，基于由固特异航空航天(GoodyearAerospace)的肯尼斯·巴彻(Kenneth Batcher)博士设计的STARAN关联SIMD计算机及其被海军舰船大量使用的后续ASPRO。

由肯特州立大学计算机科学学院开发，ASC是用于关联计算的算法模型(波特等人，计算机(Computer)，27(11)：19-25，1994)(波特，关联计算(Associative Computing)：大规模并行计算机的编程范例(A Programming Paradigm for Massively ParallelComputers)，普莱南出版公司(Plenum Publishing)，1992)。ASC模型产生于对固特异航空航天构建的关联处理器STARAN和MPP的研究。虽然目前在硬件中不受支持，但是目前有效模拟和设计用于此模型的计算机均取得了一定研究成果。

作为扩展的SIMD模型，ASC使用同步数据并行编程，从而避免了多任务和异步点对点通信路由两者。由于在任何时间仅执行一个任务，其中在所有有效处理单元(PE)上步调一致地执行此任务的多个实例，因此多任务是不必要的。如同SIMD编程器，ASC避免了包括负载均衡、同步、动态任务调度的问题、必须在MPI和其它MIND集群范例中明确处理的问题。

图4示出ASC计算机的概念模型。存在又称为指令流(IS)的单个控制单元和各自具有其自身的本地存储器的多个处理单元(PE)。控制单元和PE阵列通过广播/简化网络连接，并且PE通过PE数据互连网络连接在一起。

如图4中所见，PE可以利用位于其自身的本地存储器中的数据。数据保持在适当的位置，并且响应(有效)PE并行处理其本地数据。对词“关联”的引用涉及使用搜索以通过内容而不是存储地址来定位数据。ASC模型不采用关联存储器，替代地，其是关联处理器，其中一般循环是搜索-处理-检索。在(波特等人，计算机(Computer)，27(11)：19-25，1994)中可获得该模型的概述。

由于ASC数据结构的自然表格结构，该算法的表格性质适用于使用ASC的计算。SWAMP充分利用用于步调一致地转换北邻域和西北邻域的数据的跨PE互连网络的极为高效的通信，以及用于搜索和用于跨并行计算的最大值的快速恒定时间关联功能。

归因于ASC模型所需的额外硬件，在恒定时间执行关联操作(吉恩(Jin)等人，第15届国际并行与分布式处理会议(15th International Parallel and DistributedProcessing Symposium)(IPDPS'01)研讨会，旧金山，第193页，2001)。这些操作可通过任何SIMD类机器有效地执行(但不太迅速)，并已经成功调适以在若干SIMD硬件平台上有效地运行(Yuan等人，并行与分布式计算系统(PDCS)(Parallel and Distributed ComputingSystems(PDCS))，剑桥，MA，2009；塔汉(Trahan)等人，并行与分布式计算杂志(JPDC)(J.ofParallel and Distributed Computing(JPDC))，2009)。因此可以有效地在与SIMD密切相关的其它系统(包含向量机)上实施SWAMP和其它ASC算法，这就是该模型用作范例的原因。

控制单元提取程序指令和对程序指令进行解码并将控制信号广播到PE。PE在控制单元的指示下使用其自身的本地数据执行这些指令。所有PE以步调一致方式执行指令，其中指令之间隐式同步。ASC具有若干相关的高速全局操作：关联搜索、最大值/最小值搜索以及应答器选择/检测。这些在以下章节中进行描述。

关联功能

下文论述与SWAMP算法相关的功能。关联搜索

ASC算法中的基本操作是关联搜索。关联搜索同时定位其本地数据匹配指定搜索关键字的PE。具有匹配的数据的那些PE被称为应答器，而具有非匹配的数据那些PE被称为无应答器。在执行搜索之后，该算法随后可以限制进一步的处理以通过禁用无应答器而仅影响应答器(或反之亦然)。执行另外的搜索可进一步细化该组应答器。关联搜索由SWAMP+大量用于选择哪些PE在对角线内的并行行为中有效。

最大值/最小值搜索

除简单搜索之外，当每个PE使用标准比较运算符(等于、小于等)比较其本地数据与搜索关键字时，关联计算机还可以执行全局搜索，其中将来自整个PE阵列的数据组合在一起以确定该组应答器。最常见类型的全局搜索是最大值/最小值搜索，其中应答器是其数据为跨整个PE阵列的最大值或最小值的那些PE。最大值由SWAMP+用于其处理以追踪到目前为止所计算的最高值的每一对角线。最大值搜索的使用频繁发生，逻辑并行动作中发生一次，每次比对中发生m+n次。

应答器选择/检测

关联搜索可以产生多个应答器，并且关联算法可以三种不同的模式处理这些应答器：并行、顺序或单一选择。并行应答处理同时对每个应答器执行同一组操作。顺序应答处理个别地选择每个应答器，从而实现针对每个应答器的一组不同的操作。单一应答器选择(又称为挑选一个(pickOne))选择一个任意的所选应答器来进行处理。除多个应答器之外，关联搜索还可能产生无应答器。为了处理这种情况，ASC模型可以检测是否存在对搜索的任何应答，并且在这种情况下执行单独的一组动作(被称为任意应答器(anyResponder))。在SWAMP中，基于上述关联搜索，并行选择并处理含有待比对的字符的多个应答器。如果使用最大值/最小值搜索时存在具有完全相同的最大值的多个值，且当存在这样的多个值时，出现单一应答器选择。

PE互连网络

大多数关联处理器包含某种类型的PE互连网络以允许阵列内的并行数据移动。ASC模型本身并不指定任何特定的互连网络，而且事实上，许多适用的关联算法不要求指定特定的互连网络。通常，关联处理器实施简单的网络，例如ID线性阵列或2D网格。这些网络易于实施并且允许以同步的方式迅速传送数据。例如，ID线性阵列足以在SWAMP算法中用于PE之间的明确通信。

并行计算系统

图5中示出通用的并行处理体系结构。虽然每个组件示出为具有直接连接，但是应理解，各个单元可以是地理上分离的但是经由例如因特网的网络连接。虽然可能有混合配置，但是并行计算机中的主存储器通常或在单个地址空间中的所有处理单元之间共享，或呈分布式，即，每个处理单元具有其自身的本地地址空间。(分布式存储器是指这样的事实：存储器以逻辑方式分布，但通常意味着它也以物理方式分布。)分布式共享存储器和存储器虚拟化结合这两种方法，其中处理元件具有其自身的本地存储器以及对非本地处理器上的存储器的访问权。对本地存储器的访问通常比对非本地存储器的访问更快。

其中可以相等时延和带宽访问主存储器的每个单元的计算机体系结构被称为均匀存储器访问(UMA)系统。通常，这只能通过共享的存储器系统来实现，其中该存储器并非以物理方式分布。不具有此属性的系统被称为非均匀存储器访问(NUMA)体系结构。分布式存储器系统具有非均匀存储器访问。

可以多种方式在硬件中实施处理器-处理器和处理器-存储器通信，包含经由共享的(或多端口的或多路复用的)存储器、纵横开关、共享的总线或无数拓扑的互连网络(包含星形、环形、树形、超立方体、充足超立方体(在一个节点处具有超过一个处理器的超立方体))、或n维网格实施。

基于互连网络的并行计算机必须合并路由以实现并非直接连接的的节点之间的消息传递。用于处理器之间的通信的介质很可能在大型多处理器机器中分层。此类资源在市面上可购买用于专用用途，或者可以经由例如亚马逊的云计算的“云”访问这些资源。

计算机通常包含经由总线耦合到存储器的处理器。存储器可以包含RAM或ROM，并且优选地包含至少一个有形的非暂时性介质，其存储可执行以致使系统执行本文中所描述的功能的指令。如所属领域的技术人员在必要时将认识到或最适合于本发明的方法的性能，本发明的系统包含经由总线彼此通信的一个或多个处理器(例如，中央处理单元(CPU)、图形处理单元(GPU)等)、计算机可读存储装置(例如，主存储器、静态存储器等)、或其组合。

处理器可以是所属领域中已知的任何合适的处理器，例如由英特尔(加利福尼亚州圣克拉拉)以商标XEON E7出售的处理器，或由AMD公司(加利福尼亚州桑尼维尔)以商标OPTERON 6200出售的处理器。

存储器可以指计算机可读存储装置，并且可以包含在其上存储一个或多个指令集(例如，实施本文中存在的任何方法或功能的软件)、数据(例如，体现任何有形物理对象，例如在患者的染色体中发现的遗传序列)或这两者的任何机器可读介质。虽然在示例性实施例中计算机可读存储装置可以是单个介质，但是术语“计算机可读存储装置”应被理解为包含存储一个或多个指令集或数据的单个介质或多个介质(例如，集中式或分布式数据库、和/或相关联的高速缓冲存储器和服务器)。术语“计算机可读存储装置”应相应地被理解为包含(但不限于)固态存储器(例如，订户身份模块(SIM)卡、安全数码卡(SD卡)、微型SD卡、或固态驱动器(SSD))、光学和磁介质、以及任何其它有形存储介质。优选地，计算机可读存储装置包含有形的非暂时性介质。此类非暂时性介质不包括(例如)短暂波和信号。“非临时性存储器”应被解释为本身不包括计算机可读传输介质，例如信号。

根据本发明的输入/输出装置可以包含视频显示单元(例如，液晶显示器(LCD)或阴极射线管(CRT)显示器)、字母数字输入装置(例如，键盘)、光标控制装置(例如，鼠标或触控板)、磁盘驱动器单元、信号生成装置(例如，扬声器)、触摸屏、加速计、麦克风、蜂窝无线电频率天线、以及网络接口装置，网络接口装置可以是例如网络接口卡(NIC)、Wi-Fi卡或蜂窝调制解调器。

样本采集和制备

本发明包含用于产生对应于从生物样本回收的核酸的序列(例如，核酸序列、氨基酸序列)的方法。在一些实施例中，所得信息可以用来识别从受试者获得的核酸材料中呈现的突变。在一些实施例中，从受试者获得样本，即，核酸(例如，DNA或RNA)，处理(裂解、扩增和/或纯化)核酸并使用下文描述的方法对核酸测序。在许多实施例中，测序的结果不是线性核酸序列，而是数千或数百万个别短核酸读数的集合，其必须重新组装成针对受试者的序列。一旦比对了读数以产生序列，就可以将已比对序列与参考序列相比较以识别可以(例如)指示疾病的突变。在其它实施例中，可以基于读数与参考序列构建物(即，如上文所描述的有向非循环图(“DAG”))的比对识别受试者的特定突变。

出于任何上述目的，可以对生物样本应用各种方法。生物样本可以例如包括血液、全血、血浆、泪液、乳头吸液、血清、粪便、尿液、唾液、循环细胞、组织、活检样本、毛囊的样本或含有患者的生物材料的其它样本。在基于此类样本进行测试时的一个问题是，在大多数情况下，仅微小量的含有所关注的突变的DNA或RNA可以存在于样本中。例如腮抹试或血液样本等非侵入性样本尤其如此，其中突变异核酸存在极少量。在一些实施例中，核酸片段可以是天然的短，即，样本中的相关核酸的随机剪切可产生短的片段。在其它实施例中，为便于处理，或因为测序技术只能序列读取少于1000个碱基，例如少于500个碱基，例如少于200个碱基，例如少于100个碱基，例如少于50个碱基，所以核酸被特意分段。尽管本文中描述的方法可以用于比对不同长度的序列，但是在一些实施例中，大多数这些多个核酸读数将由测序方法产生，并且包括少于1000个碱基，例如少于500个碱基，例如少于200个碱基，例如少于100个碱基，例如少于50个碱基。

可以通过所属领域中已知的方法获取核酸。通常，可以通过各种技术从生物样本提取核酸，例如马尼亚迪斯(Maniatis)等人的分子克隆实验指南(Molecular Cloning:ALaboratory Manual)(冷泉港，纽约，第280-281页，(1982))所描述的那些技术，这些内容以全文引用的方式并入本文中。

可能有必要首先制备样本的提取物，然后执行进一步的步骤，即，微分沉淀、管柱层析、用有机溶剂提取等，以获取足够纯的核酸制备物。可以使用所属领域的标准技术通过例如细胞的化学或机械裂解来制备提取物。随后可以例如通过过滤和/或离心和/或用离液序列高的盐(例如异硫氰酸胍或尿素)或用有机溶剂(例如苯酚和/或HCCl₃)进一步处理提取物以使任何污染和潜在的干扰蛋白质改性。在一些实施例中，样本可以包括从受试者样本(例如，血液样本)收集的RNA，例如，mRNA。用于RNA提取的一般方法在所属领域中是众所周知的，并且在分子生物学的标准教科书(包含奥斯贝(Ausubel)等人的分子生物学实验室指南(Current Protocols of Molecular Biology)，约翰·威利父子出版公司(JohnWiley and Sons)(1997))中揭示。例如，在鲁普(Rupp)和洛克(Locker)的实验室投资(LabInvest.)56:A67(1987)和安德烈(De Andres)等人的生物技术(BioTechniques)18:42044(1995)中揭示了从石蜡包埋组织提取RNA的方法。这些参考文献中的每一个的内容以全文引用的方式并入本文中。具体来说，根据制造商的说明，可使用来自商用制造商(例如凯杰公司(Qiagen))的纯化试剂盒、缓冲液组和蛋白酶执行RNA分离。例如，可以用凯杰的RNeasy微型柱来分离从细胞培养的总RNA。其它市售RNA分离试剂盒包括MASTERPURE完整DNA和RNA纯化试剂盒(EPICENTRE，威斯康星州麦迪逊)和石蜡块RNA分离试剂盒(Ambion公司制造)。可以使用RNA Stat-60(Tel-Test)从组织样本分离总RNA。可以例如通过氯化铯密度梯度离心分离从肿瘤制备的RNA。

分析型测序

可以通过所属领域中已知的任何方法测序。DNA测序技术包含使用标记的终止子或引物以及板或毛细管中的凝胶分离的传统双脱氧测序反应(桑格方法)、使用可逆终止标记的核苷酸的使用合成法测序、焦磷酸测序、454测序、对标记的寡核苷酸探针的等位基因特异性杂交、对标记的克隆库使用等位基因特异性杂交的使用合成法测序(后接使用连接法测序)、在聚合步骤期间标记的核苷酸的掺入的实时监控、聚合酶克隆测序和SOLiD测序。分离分子的测序最近已通过使用聚合酶或连接酶的序列扩展或单次扩展反应以及通过使用探针库的单个或序列微差杂交得以证明。在测序之前可能另外有利的是扩增样本中的部分或全部核酸。在一些实施例中，使用所属领域中已知的聚合酶链反应(PCR)技术扩增核酸。

可以在所提供的本发明的方法中使用的测序技术的一个实例是Illumina测序(例如，MiSeq^TM平台)，它是基于聚合酶的使用合成法测序，可被用来扩增DNA或RNA。用于DNA的Illumina测序是基于使用折回PCR和锚定引物扩增固体表面上的DNA。基因组DNA被分成片断，并且在片段的5'和3'末端添加衔接子。附着到流动细胞通道的表面的DNA片段被扩展和桥式扩增。片段成为双链，并且双链分子改性。固相扩增接着改性的多个周期可以在流动细胞的每个通道中形成数百万簇的相同模板的单链DNA分子的大约1,000个拷贝。使用引物、DNA聚合酶和四荧光团标记的可逆终止子核苷酸来执行序列测序。在掺入核苷酸之后，使用激光来激发荧光团，并且捕捉图像及记录第一碱基的标识。从每个掺入的碱基中去除3'终止子和荧光团，并且重复掺入、检测和识别的步骤。当利用Illumina测序来检测RNA时，应用相同的方法，除分离并扩增RNA片段以确定样本的RNA表达之外。在用定序器询问序列之后，在例如FASTQ文件的数据文件中输出序列，该数据文件是用于存储生物序列和质量得分的基于文本的格式(见上文的论述)。

可以在所提供的本发明的方法中使用的DNA测序技术的另一实例是由生命技术公司(Life Technologies)提供的IonTorrent^TM测序。见美国专利申请第2009/0026082号、第2009/0127589号、第2010/0035252号、第2010/0137143号、第2010/0188073号、第2010/0197507号、第2010/0282617号、第2010/0300559号、第2010/0300895号、第2010/0301398号以及第2010/0304982号，这些申请案中的每一个的内容以全文引用的方式并入本文中。在Ion Torrent^TM测序中，DNA被剪切成大约300到800个碱基对的片段，并且这些片段是平端的。然后将寡核苷酸衔接子连接到片段的末端。衔接子充当用于片段的扩增和测序的引物。片段可附着到表面，并以使得这些片段可个别地解析的解析度附着。一个或多个核苷酸的添加释放了质子(H⁺)，在测序仪器中检测并记录其信号。信号强度与掺入的核苷酸的数目成正比。离子激流数据也可以作为FASTQ文件输出。

可以在所提供的本发明的方法中使用的DNA和RNA测序技术的另一实例是454^TM测序(罗氏公司(Roche))(马古利斯M(Margulies M)等人，2005，自然(Nature)，437，376-380)。454^TM测序是使用合成法测序技术，该技术还利用焦磷酸测序。DNA的454^TM测序包括两个步骤。在第一步骤中，DNA被剪切成大约300到800个碱基对的片段，并且这些片段是平端的。然后将寡核苷酸衔接子连接到片段的末端。衔接子充当用于片段的扩增和测序的引物。可以使用例如含有5'生物素标签的衔接子将这些片段附着到DNA捕获珠粒，例如，包覆抗生蛋白链菌素的珠粒。附着到珠粒的片段是在油水乳液的液滴内扩增的PCR。其结果是在每个珠粒上克隆扩增的DNA片段的多个拷贝。在第二步骤中，在孔(微微升大小)中捕获珠粒。并行对每个DNA片段执行焦磷酸测序。一个或多个核苷酸的添加产生了通过测序仪器中的CCD相机记录的光信号。信号强度与掺入的核苷酸的数目成正比。焦磷酸测序利用在核苷酸添加物上释放的焦磷酸(PPi)。在存在腺苷5'磷酰硫酸的情况下PPi由ATP硫酸化酶转化为ATP。荧光素酶使用ATP将萤光素转化为氧化萤光素，并且该反应产生被检测和分析的光。在另一实施例中，使用焦磷酸测序来测量基因表达。RNA的焦磷酸测序的应用类似于DNA的焦磷酸测序，并通过将部分rRNA基因序列的应用附着到微观珠粒中然后将附着物放置到个别孔来完成。随后扩增附着的部分rRNA序列以便确定基因谱。莎朗·马什(Sharon Marsh)，分子生物学方法中的

方案(

Protocols in Methods in Molecular Biology)，第373卷，15-23(2007)。

可以在所提供的本发明的方法中使用的DNA和RNA检测技术的另一实例是SOLiD^TM技术(应用生物系统公司(Applied Biosystems))。SOLiD^TM技术系统是可以用来运行DNA和RNA两者的大规模并行下一代测序的基于连接的测序技术。在DNA的SOLiD^TM测序中，基因组DNA被剪切成片段，并且在片段的5'和3'末端附着衔接子以产生片断库。替代地，可通过将衔接子连接到片段的5'和3'末端、使片段环化、消化环化的片段以产生内部衔接子、以及将衔接子附着到所得片段的5'和3'末端以产生配对库来引入内部衔接子。接着，在含有珠粒、引物、模板和PCR成分的微反应器中制备克隆珠粒群。在PCR之后，使模板改性并使珠粒富集以分离珠粒与扩展的模板。选定的珠粒上的模板经历3'修饰以允许接合到玻璃载片。可以通过序列杂交和连接具有由特定萤光团标识的确定中心的碱基(或碱基对)的部分随机寡核苷酸来确定序列。在记录颜色之后，分裂并去除连接的寡核苷酸，然后重复该过程。

在其它实施例中，使用SOLiD^TM基因表达系列分析(SAGE)来测量基因表达。基因表达系列分析(SAGE)是允许大量基因转录物的同步和定量分析而不需要为每个转录物提供个别杂交探针的方法。首先，产生含有足够的信息来唯一地识别转录物的短序列标签(约10-14bp)，只要该标签是从每个转录物内的唯一位置获得。然后，将许多转录物连接在一起以形成可经过测序的长系列分子，同时显示多个标签的标识。可以通过确定个别标签的丰度并且识别对应于每个标签的基因来定量评估任何转录物群的表达模式。更多细节请参考例如威尔克斯库(Velculescu)等人的科学(Science)270:484 487(1995)和威尔克斯库等人的细胞(Cell)88:243 51(1997，其中的每一个的内容以全文引用的方式并入本文中)。

可以在所提供的本发明的方法中使用的另一测序技术包含(例如)赫利克斯(Helicos)真实单分子测序(tSMS)(哈里斯T.D.(Harris T.D.)等人(2008)科学(Science)320:106-109)。在tSMS技术中，将DNA样本分裂成大约100到200个核苷酸链，并将多聚腺苷酸序列添加到每个DNA链的3'末端。每条链通过添加荧光标记的腺苷核苷酸来标记。然后将DNA链与流动细胞杂交，该流动细胞含有数百万被固定到流动细胞表面的寡-T的捕获部位。模板可以具有约1亿个模板/cm²的密度。然后将流动细胞载入至仪器(例如，HeliScope.TM.定序器)中，并且激光照射流动细胞的表面，从而显示每个模板的位置。CCD相机可以在流动细胞表面上映射模板的位置。然后分裂并洗掉模板荧光标记。通过引入DNA聚合酶和荧光标记的核苷酸开始测序反应。寡-T的核酸用作引物。聚合酶以有向模板的方式将标记的核苷酸掺入到引物。去除聚合酶和未掺入的核苷酸。通过对流动细胞表面成像来检测具有有向掺入的荧光标记的核苷酸的模板。成像之后，分裂步骤去除荧光标记，并且用其它荧光标记的核苷酸重复该过程直至实现所需读取长度。在每个核苷酸添加步骤收集序列信息。以下申请案中示出了tSMS的进一步描述：例如，在拉皮迪(Lapidus)等人(美国专利第7,169,560号)、拉皮迪等人(美国专利申请第2009/0191565号)、奎克(Quake)等人(美国专利第6,818,395号)、哈里斯(美国专利第7,282,337号)、奎克等人(美国专利申请第2002/0164629号)以及布拉斯拉维斯基(Braslavsky)等人，PNAS(USA)，100:3960-3964(2003)，这些参考文献中的每一个的内容以全文引用的方式并入本文中。

可以在所提供的本发明的方法中使用的测序技术的另一实例包含用以对DNA和RNA两者进行测序的太平洋生物科学(Pacific Biosciences)的单分子实时(SMRT)技术。在SMRT中，四种DNA碱基中的每一个可附着到四种不同萤光染料之一。这些染料是磷酸相连的。单个DNA聚合酶可与单分子模板单链DNA一起固定在零模式波导(ZMW)底部。ZMW可为能够相对于可在ZMW外迅速扩散(以微秒为单位)的萤光核苷酸背景观测单核苷酸藉由DNA聚合酶掺入的限制构建物。核苷酸掺入到生长链中耗时若干毫秒。在此期间，荧光标记被激发并产生荧光信号，并且荧光标签裂开。检测染料的对应萤光可指示掺入了哪种碱基。重复该过程。为了对RNA测序，用ZMW中的逆转录酶替换DNA聚合酶，并且相应地遵循该过程。

可以在所提供的本发明的方法中使用的测序技术的另一实例是纳米孔测序(瑞里G V(Soni G V)和米勒(Meller)，临床化学(AClin Chem)53:1996-2001)(2007)。纳米孔是直径为1纳米的数量级的小孔。纳米孔浸入在导电流体中以及跨纳米孔施加电势会因通过纳米孔的离子传导而导致轻微的电流。流动的电流量对纳米孔的大小敏感。随着DNA分子通过纳米孔，DNA分子上的每个核苷酸会不同的程度地阻碍纳米孔。因此，随着DNA分子通过纳米孔而通过纳米孔的电流的变化表示DNA序列的读取。

可以在所提供的本发明的方法中使用的测序技术的另一实例包括使用化学敏感场效应晶体管(chemFET)阵列以对DNA测序(例如，如在美国专利申请公开案第20090026082号中所描述的)。在该技术的一个实例中，DNA分子可以放入反应室中，并且模板分子可以与结合到聚合酶的测序引物杂交。通过chemFET，可以通过电流的变化检测到一个或多个三磷酸酯掺入到测序引物的3'末端处的新核酸链中。阵列可以具有多个chemFET传感器。在另一实例中，可以将单个核酸附着到珠粒，可以在珠粒上扩增核酸，并且可以将个别珠粒转移到chemFET阵列上的个别反应室，其中每个反应室具有一个chemFET传感器，并且可以对核酸测序。

可以在所提供的本发明的方法中使用的测序技术的另一实例包括使用电子显微镜(蒙德里安那基斯E.N.(Moudrianakis E.N.)和比尔M.(Beer M.)，美国国家科学院院刊(Proc Natl Acad Sci USA)，1965年3月；53:564-71)。在该技术的一个实例中，使用金属标记来标记个别DNA分子，这些金属标记是使用电子显微镜可识别的。然后在平坦表面上拉伸这些分子和并使用电子显微镜成像以测量序列。

另外的检测方法可以利用结合微阵列用于后续的荧光或无荧光检测、使用质谱方法的条码质量检测、发射的无线电波检测、已比对条码的散射光的检测、使用定量PCR或数字PCR方法的荧光检测。比较核酸杂交阵列是用于检测患者的样本DNA中的拷贝数变异的技术。使用例如不同的荧光团以不同方式标记样本DNA和参考DNA，然后将样本DNA和参考DNA与众多探针杂交。随后测量样本和参考的荧光强度，随后使用荧光强度比来计算拷贝数变异异。在希纳维M(Shinawi M)、张SW(Cheung SW)的阵列CGH及其临床应用(The array CGH and its clinical applications)，今日药物发现(Drug Discovery Today)13(17-18)：760-70中更详细论述了比较基因组杂交阵列的方法。微阵列检测可以不直接产生FASTQ文件，但是可使用程序将由微阵列定序器产生的数据转化为FASTQ格式或类似的格式。

检测DNA分子、RNA分子和拷贝数的另一种方法是荧光原位杂交法(FISH)，原位杂交方案(伊恩·达比(Ian Darby)编，2000年)。FISH是一种分子细胞遗传学技术，其检测特定的染色体重排，例如DNA序列的突变和拷贝数变化。DNA分子经过化学改性并分离成两链。随后用该改性的DNA链培育单链探针。信号链探针取决于目标序列部分而选择，并具有与互补序列部分的高亲和性。探针可以包含重复序列探针、全染色体探针和基因座特定探针。在培育时，将组合的探针和DNA链杂交。然后在显微镜下将结果可视化并进行定量以评估任何变异。

在另一实施例中，使用基于MassARRAY^TM的基因表达谱方法来测量基因表达。在西格诺公司(Sequenom,Inc.)(加利福尼亚州圣地亚哥)根据RNA和逆转录的分离开发的基于MassARRAY^TM的基因表达谱方法中，获取的cDNA被掺入合成DNA分子(竞争者)，其匹配所有位置中的目标cDNA区(除单一碱基外)并用作内标。cDNA/竞争者混合物经过PCR扩增并经受后PCR虾碱性磷酸酶(SAP)的酶处理，这使得剩下的核苷酸去磷酸化。在碱性磷酸酶失活后，来自竞争者和cDNA的PCR产物经历引物扩展，这产生针对竞争者-cDNA-派生PCR产物的不同质量信号。纯化后，在芯片阵列上这些产物，该芯片阵列预加载有用基质辅助激光解吸电离飞行时间质谱(MALDI-TOF MS)分析进行分析所需的成分。存在于反应中的cDNA随后通过分析所产生的质谱的峰面积比来对反应中存在的cDNA进行定量。进一步的细节请参考，例如，丁(Ding)和坎托(Cantor)，美国国家科学院院刊(Proc.Natl.Acad.Sci.USA)100:3059 3064(2003)。

其它基于PCR的技术包括：例如，差异显示(Liang和帕迪(Pardee)，科学(Science)257:967971(1992))；扩增片段长度多态性(iAFLP)(川本(Kawamoto)等人，基因组研究(Genome Res.)12:1305 1312(1999))；BeadArray^TM技术(Illumina公司，加利福尼亚州圣地亚哥；奥利芬特(Oliphant)等人，疾病标记发现(生物学技术补充)(Discovery of Markersfor Disease(Supplement to Biotechniques))，2002年6月；弗格森(Ferguson)等人，分析化学(Analytical Chemistry)72:5618(2000))；用于基因表达检测的珠阵(BADGE)(BeadsArray for Detection of Gene Expression (BADGE))，在基因表达的快速检测中使用市售的Luminex 100 LabMAP系统和多个彩色编码微球(路明克斯公司(Luminex Corp.)，德克萨斯州奥斯汀)(Yang等人，基因组研究(Genome Res.)11:1888 1898(2001))；以及高覆盖率表达谱(HiCEP)分析(福村(Fukumura)等人，核酸研究(Nucl.Acids.Res.)31(16)e94(2003))。这些文献中的每一个的内容以全文引用的方式并入本文中。

在某些实施例中，还可以识别或使用微阵列技术确认基因表达的差异，微阵列技术包含例如尼龙膜阵列、微芯片阵列和玻璃载片阵列，这些阵列例如可从昂飞公司(Affymetrix)(加利福尼亚州圣克拉拉市)购得。通常，RNA样本被分离，并经由逆转录转化为标记的cDNA。随后使用来自所关注的细胞或组织的特定DNA探针将标记的cDNA杂交到或者尼龙膜、微芯片上或者玻璃载片上。然后检测杂交的cDNA并进行量化，可以将所得的基因表达数据与对照物进行比较以用于分析。标记、杂交和检测的方法取决于微阵列支撑物是尼龙膜、微芯片还是玻璃载片而不同。尼龙膜阵列通常用P-dNTP标记的探针进行杂交。玻璃载片阵列通常涉及用两个不同的荧光标记的核苷酸来标记。在耶特曼(Yeatman)等人(美国专利申请案第2006/0195269号)中示出了制作微阵列和确定基因产物表达(例如，RNA或蛋白质)的方法，该申请案的内容以全文引用的方式并入本文中。

在一些实施例中，可单独使用或结合其它方法(例如，免疫测定或RNA测量测定)使用质谱(MS)分析，以确定本文所公开的一个或多个生物标志物在生物样本中的存在和/或数量。在一些实施例中，MS分析包含基质辅助激光解吸/电离(MALDI)飞行时间(TOF)MS分析，例如，直接点MALDI-TOF或液相层析MALDI-TOF质谱分析。在一些实施例中，MS分析包括电子喷雾离子化(ESI)MS，例如液相层析(LC)ESI-MS。可使用市售的光谱仪来完成质量分析。使用包含MALDI-TOF MS和ESI-MS的MS分析来检测生物样本中生物标记肽的存在和数量的方法在所属领域中是已知的。进一步的指导请参考例如美国专利第6,925,389号、第6,989,100号和第6,890,763号，这些专利中的每一个以全文引用的方式并入本文中。

可以使用相关领域的技术人员已知的多种技术来确定在本发明的方法、序列构建物和系统中使用的蛋白质序列。例如，可以通过以质谱分析蛋白质或蛋白质的一部分或使用埃德曼降解产生氨基酸序列和氨基酸序列读数。质谱可以包含：例如，基质辅助激光解吸/电离(MALDI)飞行时间(TOF)MS分析，例如直接点MALDI-TOF或液相层析MALDI-TOF质谱分析；电子喷雾离子化(ESI)MS，例如液相层析(LC)ESI-MS；或其它技术，例如MS-MS。可以使用商用仪器(例如，型号49X精确蛋白质/肽定序器((应用生物系统公司(AppliedBiosystems)/生命技术公司(Life Technologies))执行埃德曼降解分析。测序后的氨基酸序列(即，多肽，即，蛋白质)的长度可以为至少10个氨基酸，例如，长度至少为20个氨基酸，例如，长度至少为50个氨基酸。

参考文献并入

在本发明通篇中已经参考并且引用了其它文献，例如专利、专利申请案、专利公开案、杂志、书籍、论文、网络内容。所有此类文献在此出于所有目的以全文引用的方式并入本文中。

等效物

根据包含对本文中引用的科学和专利文献的参考的本文档的完整内容，所属领域的技术人员将显而易见除本文示出且描述的之外的本发明的各种修改以及其许多进一步实施例。本文中的标的物含有重要信息、范例和指南，其可适于本发明在其各种实施例及其等效物中的实践。

Claims

1.一种用于比对多个序列读数的系统，其包括处理器和非暂时性存储器，其中所述存储器包括指令，所述指令当被执行时致使所述处理器：

获取作为成串符号的多个序列读数；

比较对应于多个序列读数中的序列读数的每串符号与包括多个节点的参考有向非循环图(DAG)中的多个位置，其中每个节点以包含表示序列的一个或多个符号的串，界定所述节点相对于所述参考DAG中的其它节点的位置的父节点集，以及位置的形式存储在存储器中，其中遗传变异由所述参考DAG中的多个位置处的替代节点表示，

对每个序列读数与所述参考DAG中的所述多个位置中的每一个位置之间的重叠进行评分，其中更高得分对应于更大重叠量；

其中，所述评分进一步包括通过识别所述参考DAG中的先前位置在每个序列读数与所述参考DAG的一部分之间识别比对，其中节点的每串符号的前趋通过以下方式识别：(i)如果符号不是其节点的串符号的首个符号，那么在符号前面的符号是其唯一的前趋；或(ii)如果符号是其节点的串符号的首个符号，那么其节点的父节点的任何节点的串符号的最后一个符号是前趋；

基于所述评分将每个序列读数指派到所述参考DAG上的位置；以及

将对应于每个已比对序列读数的位置的文件写入到存储器。

2.根据权利要求1所述的系统，其中所述序列读数是核酸序列读数。

3.根据权利要求1所述的系统，其中所述序列读数是氨基酸序列读数。

4.根据权利要求1所述的系统，其中所述存储器另外包括指令，所述指令当被执行时致使所述处理器：

基于指派到每个序列读数的位置将多个已比对序列读数组装为组装序列；以及

将对应于所述组装序列的文件写入到存储器。

5.根据权利要求4所述的系统，其中所述组装序列对应于生物体的遗传序列。

6.根据权利要求5所述的系统，其中所述遗传序列涵盖染色体或基因组。

7.根据权利要求1所述的系统，其进一步包括基于一个或多个已比对序列读数的位置为生物体指派基因型。

8.根据权利要求7所述的系统，其进一步包括使所述指派基因型与疾病风险相关联。

9.根据权利要求8所述的系统，其中所述疾病是癌症。

10.根据权利要求1所述的系统，其中所述参考DAG包括存储在非暂时性计算机可读介质上的位置和符号的数据库。

11.根据权利要求1所述的系统，其中所述系统包括多个处理器，并且其中每个处理器经配置以比较所述多个序列读数的一部分与所述参考DAG并对其进行评分。

12.根据权利要求1所述的系统，其中所述存储器另外包括指令，所述指令致使第一处理器指示第二处理器进行获取、比较、评分、识别和指派，且随后发送指派到每个序列读数的位置到所述第一处理器。

13.根据权利要求1所述的系统，其中所述参考DAG上的每个位置对应于参考序列中的一个核酸。

14.根据权利要求1所述的系统，其中所述DAG上的每个位置对应于参考序列中的一个氨基酸。

15.根据权利要求1所述的系统，其中所述DAG上的每个位置对应于参考序列中的一个基因。

16.根据权利要求1所述的系统，其中所述参考DAG具有方向。

17.根据权利要求16所述的系统，其中在所述方向上贯穿所述参考DAG的路径表示生物体的基因组。

18.根据权利要求16所述的系统，其中在所述方向上贯穿所述参考DAG的路径表示生物体的染色体。

19.根据权利要求16所述的系统，其中在所述方向上贯穿所述参考DAG的路径表示蛋白质。

20.根据权利要求16所述的系统，其中在所述方向上贯穿所述参考DAG的路径是非循环的。

21.根据权利要求1所述的系统，其中所述DAG中的所述替代节点表示在所述位置处的已知序列变异。

22.根据权利要求1所述的系统，其中所述替代节点对应于核酸序列，因碱基插入或碱基删除而彼此不同。

23.根据权利要求1所述的系统，其中所述替代节点对应于多肽，因氨基酸插入或氨基酸删除而彼此不同。

24.根据权利要求1所述的系统，其中所述参考DAG包括多于1,000个符号。

25.根据权利要求24所述的系统，其中所述参考DAG包括多于1,000,000个符号。

26.根据权利要求1所述的系统，其中所述多个序列读数包括多于1,000个序列读数。

27.根据权利要求1所述的系统，其中所述多个序列读数的至少一部分的长度多于100个符号。

28.根据权利要求1所述的系统，其中所述多个序列读数是下一代测序方法的输出文件。