CN105849555B - 用于错误校正的序列读数迭代聚类 - Google Patents

Abstract

示例性实施方式提供了用于错误校正的序列读数的迭代聚类的方法和系统。示例性实施方式的方面包括接收序列读数的集和相关的质量值；将序列读数基于序列相似性分组为初始簇的集；生成各初始簇的簇共有区；基于与序列读数相关的质量值和簇共有区迭代改进聚类；并且生成并输出各簇的最终簇共有区。

Description

用于错误校正的序列读数迭代聚类

相关申请的交叉引用

本申请要求2013年12月18日提交的题为“用于从混合群体中生成共有序列的方法”的美国临时专利申请系列号61/917,777，和2014年7月24日提交的题为“用于错误校正的序列读数迭代聚类”的美国临时专利申请系列号62/028,741的权益，两者转让给本申请的受让人，并且通过引用纳入本文。

发明背景

生物分子序列确定的进步，尤其是针对核酸和蛋白质样品，已经彻底改变了细胞和分子生物学领域。由自动化测序系统发展促进，现在能够对样品核酸的混合群进行测序。然而，序列信息的质量必须得到仔细监控，并且可被与生物分子本身或使用的测序系统相关的许多因素削弱，包括生物分子的组成(例如，核酸分子的碱基组成)、实验和系统杂音、观察的信号强度的变化、和反应效率的差异。如此，必须采用方法以分析和改善来自这类测序技术的数据的质量。

除了影响生成的序列读数的总体准确性以外，这些因素可能使碱基判定为真变异或者错判(例如，序列读数中的插入、删除或错配错误)的设计复杂化。例如，当序列读数具有在同源染色体之间不同的碱基判定时，能够确定不同的碱基判定是同源物之间的真变异或者仅仅是测序错误是重要的。另外，个体中的病毒群体可能在群体中的个体病毒基因组之间存在许多变异，尤其是可高度突变的病毒，如HIV。能够鉴定具有不同来源(例如，不同染色体或基因组来源)的测序读数是能够精确表征核酸混合群体的关键。对于生成100％精确的读数的理论测序平台而言，读数可简单地用简单字符串比对算法与另一读数比较。读数之间的任何差异表明真变异，因此表明不同来源。然而，任何现实原始测序数据可能含有错误，因此简单字符串匹配算法方法是不足的。当对转录组进行测序时尤为如此。

转录组是所有RNA分子的集合，包括一种细胞或细胞群体生成的mRNA、rRNA、tRNA和其他非编码RNA。因为该数据包括细胞中的所有mRNA转录本，转录组反映了在任何给定时间上正处于主动表达的基因。目前，有两种推导转录组的一般方法。一种方法将序列读数映射到转录组正受到研究的生物体或紧密相关的物种参照基因组上。另一种方法是转录组从头组装，其使用软件来从短序列读数直接推导出转录本。

然而，市售的基因组比对器不能对转录组测序中的全长度长序列读数进行错误校正。例如，在RS II设备上产生的读数平均为5-6kb，并且常规生成长达20kb的读数。对于这种长读数能力，可对全长mRNA转录本进行测序，例如，在转化为cDNA之后。这可有助于研究人员鉴定难以使用短读数测序技术重构的剪接模式。然而，公共可得的序列比对器，例如GMAP，和功能性注释工具几乎都需要具有接近100％准确性的读数。PacBio设备从具有使其难以直接应用这些序列比对工具的错误特征的单一模板分子生成读数。然而，在测序插入(转录本)远短于聚合酶阅读长度的情况中，可生成高度准确的共有序列：通过对单一分子的冗余测序，cDNA模板的长度与系统中聚合酶测序引擎的处理能力的组合可产生充足的冗余以实现这些分析工具所需的准确性。然而，这仅可应用于较短的转录本，而较长的转录本仍然在它们实现适合生物分析的准确性水平之前需要额外处理。

目前，存在用于在转录组测序中对长读数(例如，cDNA长读数)进行错误校正公开的2种工具，PacBioToCA和LSC。两种工具都使用短读数(例如，短读数)，并且按照以下一般方案：对于各长读数，将短读数与长读数比对仿佛其是基因组“支架”，并且基于短读数比对生成最佳共有区。这种一般方案有几个缺陷：(1)由于短读数仅为50-10bp，它们可能非特异性映射并且引入更多的错误；(2)所有现有的短读数技术携带它们自身的系统错误，其可能使校正产生偏差；(3)没有利用相同的转录本通常由多个长读数表示的事实，其在来自太平洋生物科学公司(Pacific Biosciences)的长读数的情况没有系统系统偏差；(4)没有使用来自长读数的质量值(QV)；和(5)该方案需要2种不同的测序系统。

需要一种解决转录组测序中错误的问题的算法，并且优选设计为处理从头合成转录组，即没有参照基因组的算法。

发明内容

示例性的实施方式一般涉及用于分析来自核酸的混合群体的序列数据、用于将各序列读数分配到特定来源、和用于最终鉴定来自序列信息的一个或多个生物分子目标序列的一个或多个共有序列的方法。本文提供的方法不仅可应用于几乎没有错误的序列数据，也可应用于具有较高频率的插入、删除和/或错配错误的序列数据。因此，本发明还涉及进行这些方法的系统。

参照以下详细说明和附图将更好地理解本发明和各种具体方法及实施方式，其中，在各种具体方面和实施方式中描述了本发明。提供这些是为清楚起见，并且不用被认为限制本发明。本发明及其方面可应用于多种类型的本文未具体公开的方法、装置和系统。在某些方面中，示例性实施方式提供了用于错误校正的序列读数的迭代聚类的方法和系统，其通过在至少一个处理器上进行的至少一个软件组件来进行。在某些实施方式中，这类方法包括接收序列读数的集和相关的质量值；将序列读数基于序列相似性分组为初始簇的集；生成各初始簇的簇共有区；基于与序列读数相关的质量值和簇共有区迭代改进聚类；并且生成并输出各簇的最终簇共有区。

在另一个方面中，迭代改进聚类还包括：使用质量值计算属于各簇的各序列读数的概率；将个体序列读数从一个簇重分配至具有最高计算概率的另一个簇；并且合并高度相似的簇。

在一个实施方式中，输入序列读数包括长度至少0.5kb至长度1、2、3、4、5、7或10kb的全长的长读数，并且使用簇共有区和非全长读数来生成最终簇共有区，其可用于提供序列数据的全覆盖率以提供更高水平的共有区。

附图的一些方面的简要说明

图1是显示用于实施使用用于转录组测序数据的错误校正的测序读数的迭代聚类的进程的计算机系统的一个实施方式的图。

图2是显示按照示例性实施方式用于错误校正的序列读数的迭代聚类的方法的某些方面的流程图。

图3是显示来自已经比对以产生成对比对的相同同种型的2个读数的示例性部分的图。

图4是显示示例性相似性图像的图。

图5是显示一个用于区分比对的读数之间的真同种型差异与序列错误的实施方式的图。

图6是显示初始分配至错误簇的序列读数的示例的图，其中相同填充模式的序列读数来自相同同种型。

图7是显示分别针对各簇生成的示例性簇共有区C1、C2、C3和C4的图。

图8是显示将序列读数从一个簇重分配至具有最高的成员计算概率的簇的图。

图9是显示从孤儿产生新簇的示例的图。

图10是显示2个簇合并的图。

发明详述

本发明的多个实施方式和组分采用在多个技术领域熟悉的信号和数据分析技术。为了清楚地说明，本文不提供已知分析技术的详细内容。这些技术描述于多个可及的参考文献中，如：R.B.Ash，《真实分析和概率》(Real Analysis and Probability)，学术出版社(Academic Press)，纽约，1972；D.T.Bertsekas和J.N.Tsitsiklis，《概率介绍》(Introduction to Probability)，2002；K.L.Chung，《固定转移概率的马尔科夫链》(Markov Chains with Stationary Transition Probabilities)，1967；W.B.Davenport和W.L Root，《随机信号和噪音理论介绍》(An Introduction to the Theory of RandomSignals and Noise)，麦格劳-希尔公司(McGraw-Hill)，纽约，1958；S.M.Kay，《统计学处理基础》(Fundamentals of Statistical Processing)，第1-2卷，(精装-1998)；MonsoonH.Hayes，《统计学数据信号处理和建模》(Statistical Digital Signal Processing andModeling)，1996；R.M.Gray和L.D.Davisson的《统计学信号处理介绍》(Introduction toStatistical Signal Processing)；Steven M.Kay的《现代光谱估计：理论及应用》(ModernSpectral Estimation：Theory and Application)/书和光盘(Prentice-Hall信号处理丛书)(精装-1988年1月)；Steven M.Kay的《现代光谱估计：理论及应用》(Modern SpectralEstimation：Theory and Application)(平装-1999年3月)；Burkhard Buttkus的《应用地球物理中的光谱分析和过滤理论》(Spectral Analysis and Filter Theory in AppliedGeophysics)(精装-2000年5月11日)；Donald B.Percival和Andrew T.Walden的《物理应用的光谱分析》(平装-1993年6月25日)；J.L.Starck和F.Murtagh的《天文图像和数据分析》(Astronomical Image and Data Analysis)(天文学和天体物理学图书馆)(精装-2006年9月25日)；Daniel S.Sem的《蛋白质组学中的光谱技术》(Spectral Techniques InProteomics)(精装-2007年3月30日)；Dhammika Amaratunga和Javier Cabrera的《DNA微阵列和蛋白质阵列数据的探索和分析》(Exploration and Analysis of DNA Microarrayand Protein Array Data)(概率和统计学Wiley丛书)(精装-2003年10月21日)。

转录组分析的长读数错误校正与基因组组装的错误校正不同。两者都可归结为聚类问题。在基因组组装中，仅存在与染色体一样多的“簇”；各染色体彼此非常不同。与整个染色体尺寸比较，共享的重复区域非常小，并且只要存在跨重复的连续长读数，相对容易决定其起点。

相反，对于转录组分析，存在与转录本一样多的簇。在真核生物中，基因可具有许多不同的剪接形式。在一个极端示例中，转录本的一个同种型有额外的20bp外显子，而其他同种型则没有。对于许多生物学问题而言，能够将2种同种型区分开是重要的。这种详细差异的水平很少在基因组规模上发现，因此，现有的方法，例如生成高质量(>99.999％准确)从头组装的分级基因组组装过程(HGAP)不能直接应用于转录组问题(HGAP描述于2013年7月12日提交的美国专利申请13/941,442)。

“准种问题”是一般转录组聚类问题的具体应用。像转录组测序那样，簇的总量是未知的并且必须迭代地“猜测”簇和簇共有区。对于HIV基因组而言，该问题更简单，因为HIV基因组是已知的并且目前可着眼于预期的突变数量。Zogardi等，(2010)“HIV准物种的可靠估计和下一代测序数据的错误校正(Error correction of next-generation sequencingdata and reliable estimation of HIV quasispecies)”Nucl.Acids.Res.doi：10.1093/nar/gkq655中提供了准物种问题的其他信息，其通过引用全文纳入本文用于所有目的。

按照示例性实施方式，提供了解决转录组测序中的错误的问题的算法。然而，与使用“种子读数(seed read)”来比对较短读数以生成高度精确的预组装读数的HGAP概念不同，示例性实施方式的算法采用簇共有区。

示例性实施方式一般涉及从混合群体中生成共有序列。更具体地，示例性实施方式提供了基于主要使用长读数数据的同种型迭代聚类对读数进行错误校正的方法和系统。迭代计算各输入序列读数属于各簇的概率，然后将序列重分配至具有更高成员概率的簇。另外，该进程合并高度相似的簇。按照示例性实施方式，迭代同种型水平聚类去除了转录本冗余并且改进了转录组共有区准确性，全部都不需要参照基因组。

计算机实施

图1是显示用于实施用于错误校正的序列数据迭代聚类的进程的计算机系统的一个实施方式的图。在具体实施方式中，本发明可整体或部分体现在固定介质上记录的软件。计算机100可以是具有至少一个处理器102(例如，CPU等)、存储器103、输入/输出(I/O)104和数据储存库106的任意电子装置。CPU102、存储器103、I/O 104以及数据存储库106可通过系统总线，或者使用任意类型的通信连接来连接。虽然未显示，计算机100也可包括用于有线和/或无线通信的网络接口。在一个实施方式中，计算机100可包括个人计算机(例如，台式机、笔记本、平板等)、服务器、客户端计算机，或可穿戴装置。在另一个实施方式中，计算机100可包括任意类型的用于与远程数据应用相互作用的信息电器，并且能够包括这类装置如互联网功能电视、手机等。

处理器控制计算机100的运行并且可从存储器103和/或数据存储库106读取信息(例如，指令和/或数据)并且相应执行指令以执行示例性实施方式。术语“处理器102”往往包括一个处理器、多个处理器，或者一个或多个多核处理器。

I/O 104可包括任意类型的输入装置，如键盘、鼠标、麦克风等，以及任意类型的输出装置，例如，监视器和打印机。在计算机100包括服务器的一个实施方式中，输出装置可耦合至本地客户端计算机。

存储器103可包括任意类型的静态或动态存储器，包括闪存、DRAM、SRAM等。存储器103可存储程序和数据，包括序列比对器/重叠器110、簇共有区算法111、迭代簇错误校正(ICE)组件112，和平滑组件114(例如，Quiver)。这些组件/算法可用于本文所述的转录组序列组装进程。

数据存储库106可存储几个数据库，包括其存储序列读数116、读数质量值(下文中QV)118、最大团120、簇122、簇共有区124、概率126，和最终簇共有区128的一个或多个数据库。造转录组测序实施方式中，序列读数116包括同种型序列读数，其可包括全长序列读数(下文“全长读数”)116-1和非全长序列读数(下文“非全长读数”)116-2。同样，在该实施方式中，簇122可包括同种型水平簇。

在一个实施方式中，数据存储库106可位于计算机100内。在另一个实施方式中，数据存储库106可通过网络端口或外部装置连接至计算机100。数据存储库106可包括分离的服务器或任意类型的存储装置(例如，盘型光学或磁性介质、固态动态或静态存储器等)。数据储存库106可任选地包括多个辅助存储装置，例如，用于分开存储输入序列(例如，序列读数)、序列信息、计算结果和/或其他信息。计算机100可在此后使用该信息来指导服务器或客户端逻辑，如本领域所理解的那样，以体现本发明的方面。

操作中，操作者可通过显示屏(未显示)上呈现的用户界面与计算机100相互作用以序列读数116和各种软件程序所需的其他参数。一旦援用，包括序列比对器/重叠器110、簇共有区算法111、ICE组件112和平滑组件114的存储器103中的程序由处理器102执行以实施本发明的方法。

序列比对器/重叠器110从数据存储库106中读取选择的序列读数116并且在序列读数116上进行序列比对以鉴定相似的区域，其可以是结构或功能或其他序列读数116之间的关系的结果。在一个实施方式中，全长读数116-1一般是高准确性读数，例如，至少约98％或99％准确，并且可以是来自提供这种高质量读数的测序技术的原始读数，或者可以是构建自较低质量的测序读数数据的预组装的高质量读数，如本文他处所述。比对的序列117在序列比对期间由序列比对器/重叠器110生成。在某些实施方式中，序列比对器/重叠器110以C、C++、Java、C#、F#、Python、Perl、Haskell、Scala、Lisp、Python/C混合式和本领域已知的其他语言执行。

ICE组件112通过基于相似性和最大团120将序列读数116分成初始簇的集来生成类似序列读数的簇122。簇共有区算法111生成各簇的簇共有区124。ICE组件112然后通过基于簇共有区124和与序列读数相关的质量值的迭代来迭代改进聚类，其包括基于概率126将序列读数116从一个簇重分配至另一个簇，并且合并基本相似的簇。然后，平滑组件114可按照示例性实施方式生成各簇122的最终簇共有区128，如下文进一步所述。

该处理的输出可包括最终簇共有区128的列表，其各自代表簇的“共有区”。在一个实施方式中，各簇122可代表单一、独特的转录本。因此，在一个实施方式中，本发明可提供使用全长读数116-1从混合群体中鉴定独特的全长转录本的集的方法和系统。

在一个实施方式中，处理的结果还可任选地包括质量信息、技术信息(例如，峰特征、预期的错误率)、替代(例如，第二或第三好)共有区确定、置信标准等。在产生初始簇、生成簇共有区、迭代聚类和生成最终簇共有区的进程期间和之后，这一处理的过程和/或结果可保存到存储器103和数据存储库106和/或通过I/O 104输出用于在显示装置上显示和/或保存到其他存储装置(例如，CD、DVD、蓝光、闪存卡等)，或打印。

图2是显示按照示例性实施方式用于错误校正的序列读数的迭代聚类的进程的某些方面的流程图。在一个实施方式中，可使用该进程来校正转录组测序期间长读数中的错误。该进程可通过序列比对器/重叠器110、簇共有区算法111、ICE组件112、和平滑组件114的组合来进行(图1)，其虽然显示为分开的组件，各自的功能可合并成较少或较大数量的软件算法/组件。

可通过接收一组序列读数116和相关的读数质量值118(块200)来开始该进程。序列读数116优选包括但不限于一组全长的长读数116-1。读数质量值(QV)118是由测序机器生成的对每个位置碱基判定准确性的估计。

迭代聚类错误校正(ICE)组件112基于序列相似性(块202)将序列读数分成初始簇的集。簇共有区算法111生成各初始簇的暂时簇共有区(块204)。ICE组件112基于与序列读数相关的读数质量值118和簇共有区迭代改进聚类(块206)，如下文进一步详述。

在其他实施方式中，该进程还包括生成并输出各簇的最终簇共有区128(块208)。在一个实施方式中，最终簇共有区128可包括最终簇共有序列的列表，其各自代表簇的共有序列(并且因此在一个实施方式中代表转录本)。在一个实施方式中，一旦完成了迭代聚类进程，即可生成最终簇共有区128。在输入包括全长读数116-1的另一个实施方式中，可通过将非全长读数116-2输入最终平滑进程，其随后生成最终簇共有区128来生成最终簇共有区。如本领域所熟知，最终簇共有区128可保存到，例如，存储器103和/或数据存储库106，或送至I/O 104用于在监视器上显示和/或由打印机打印。

上述步骤的进一步详细描述见下。

序列读数

转录本同种型测序的一个目的是使用准确、未组装的、全长的长读数理解转录组复杂性。通过测序机器自动捕获并鉴定全长读数116-1，但是示例性实施方式通过迭代聚类提高准确性。

在示例性实施方式中，输入序列读数116包括例如转录本的全长的长读数116-1。然而，在另一个实施方式中，输入序列读数116可包括非全长读数116-2。序列读数116可任选地包括冗余序列信息，例如，其中相同的转录本经重复测序以生成包括转录本的多个拷贝的长序列读数。此外，与序列读数116相关的其他信息可包括相关的测序技术输出的特征(例如，踪迹特征(积分的每峰计数、峰的形状/高度/宽度、与相邻峰的距离、相邻峰的特征)、信噪比、功率噪音比、背景标准、信号强度、反应动力学等)等。

初始聚类

迭代聚类进程包括2个阶段。第一阶段包括基于序列相似性将序列读数116分成初始簇的集(块202)。在一个实施方式中，例如，使用初始聚类有助于确定那些序列读数116来自相同的转录本同种型。聚类的背景想法是对源自相同同种型的多个拷贝的多个序列读数的观察。例如，以下显示了源自相同同种型的转录本读数的3个拷贝：

其中，加粗的碱基表示主要是插入缺失标记(插入或缺失)的随机分布错误。对这类源自相同同种型的读数的聚类可能产生更高准确性的共有序列。

比对

图2显示了初始聚类进程的进一步详细说明(块202)。在一个实施方式中，初始聚类进程可通过由序列比对器/重叠器110比对序列读数116开始产生比对的读数。可使用许多已知的序列比对进程，例如，使用基本局部比对与连续细化(Basic Local Alignmentwith Successive Refinement)(BLASR)算法映射单一分子序列读数116，进一步描述于美国专利公开号20120330566，其通过引用全文纳入本文用于所有目的。

图3是显示来自已经使用序列比对器/重叠器110比对以产生比对的读数300的相同同种型的2个读数的示例性部分的图。在该实施例中，显示为“查询”的第一比对的读数的长度是1,675bp，并且显示为“目标”的第二比对的读数的长度是1,680bp。比对的读数300之间的比对(“nMatch”)是1.6kbp并且相似性百分比(“％sim”)是99.1677，其包括2个插入和11个缺失的插入缺失标记(表示为“*”)。

比对之后，下一个步骤是形成同种型簇。可使用参照基因组并且将读数与参照基因组比对并且确定位于特定基因座的读数代表同种型。然而，每个基因座存在许多尚未确定的替代性同种型。另外，该方法依赖于比对器并且需要以好的参照基因座开始，这限制了该方法应用于具有已存在的参照基因组的那些应用。

按照一个示例性实施方式，提供了不使用参照基因组鉴定同种型簇，并因此适用于不存在参照基因组的应用的方法和系统。

相似性图像

再次参考图2，在比对之后，使用比对的读数300来构建相似性图像(块202-2)。构建相似性图像，使得各序列读数116表示为图像中的节点，并且序列读数116之间的比对表示为节点之间的连接边缘，以显示2个序列读数具有比对命中(即，足够高的相似性百分比)。

图4是显示示例性相似性图像400的图。用于发现同种型聚类的算法采用成对比对，其中相似性图像400中的各节点402表示读数，并且连接节点对的边缘404表示存在成对比对，如图3所示，其中查询和目标读数由于其高的相似性百分比将在图像中表示为节点402的对并且通过边缘404连接。

最大团

一般而言，类似性图像进程导致形成多个相似性图像400。再次参考图2，之后在相似性图像中发现所有最大团(块202-3)。团(chique)是指包括节点的集的图像，其中对于每2个节点402存在连接两者的边缘。最大团是不含与其他团重叠的节点402的最大可能尺寸的团。最大团发现算法非确定性地将相似性图像400划分成不重叠的最大团。存在许多发现所有最大团的方法。在一个实施方式中，可运行最大团发现算法，如贪婪随机自适应检索法，其迭代构建随机化的、贪婪偏差的解决方案，其然后扩大到局部最优解决方案。参见例如，Abello等，On maximum clique problems in very large graphs(非常大图像中最大团问题的研究)，AT&T实验室研究技术报告(AT&T labs Research Technical Report)，1998，其通过引用全文纳入本文用于所有目的。

将相似性图像400划分成非重叠最大团需要比较序列读数116以检测同种型比对差异以确定序列读数116是否属于相同的团。一种检测同种型比对差异的方法是检测2个比对的读数之间的比对中的大间隙。例如，如果考虑2个比对的读数，其中一个相对另一个有大插入，则非常可能插入是额外的外显子，并且因此，可检测到同种型差异。然而，检测同种型比对差异随着比对中的间隙变得越来越小而变得越来越有问题。例如，2个比对的读数之间仅7个碱基的插入差异可能表示聚合物延伸。需要确定的是，这是真同种型差异还是序列错误。

按照示例性实施方式的一个方面，可通过促使来自包括插入的原始序列读数116的各碱基与估计每个位置的准确性并且显示各碱基是取代错误、插入错误或缺失错误的概率的质量值(QV)相关的事实改变来从序列错误中确定同种型差异。

图5是显示一个用于区分比对的读数300之间的同种型差异与序列错误的实施方式的图。在一个实施方式中，可使用差异阵列500来保持2个比对的读数300之间的位置差异的踪迹。在各碱基位置处有“+”的2个比对的读数300之上和之下的取代(S)、插入(I)和缺失(D)的行显示相关的QV 118表示足够可能出现错误的位置。差异阵列500中的各位置可包括值，例如，0或1，其中0值表示2个比对的读数300之间的差异是由于测序错误造成的(而不是真同种型差异)，并且1值表示2个比对的读数300之间的差异不能由测序错误解释。

然后确定在差异阵列中是否存在任何足够大的1值的区域，即从[I，J]中寻找大于或等于阈值长度T的范围，并且1值的区域的总和大于阈值百分比C的差异阵列中的区域：

J-I≥T且(D[I：J])的加和≥C*T

例如，估计阈值长度T设为10个碱基，并且阈值百分比C设为50％。将检索差异阵列500中长于10个碱基的其中超过50％的碱基具有1值的区域。如果无法发现这种区域，则2个比对的读数300可被认为来自相同同种型。在图5所示的示例中，差异阵列500中不存在这种区域，使得2个比对的读数300被确定为来自相同同种型并且因此置于同一团中。如果，另一方面，发现这种区域，则将确定2个比对的读数300来自不同的同种型并且因此不会放在同一团中。对于其他信息，参见Tseng和Tompa，用于在多个序列比对中定位极端保守的元件的算法(Algorithms for Locating Extremely Conserved Elements in MultipleSequence Alignments)，BMC Bioinformatics(2009)，其通过引用全文纳入本文用于所有目的。

注意根据定义，团需要各节点402互相连接。按照示例性实施方式，在最大团发现进程之后，术语“团”将置于更宽的术语“簇”之下，因为在最大团发现进程之后不需要或使用节点之间的边缘404。

在比对序列读数(块202-1)之后，构建相似性图像和最大团发现进程(块202-2到202-3)，如何分组序列读数116的问题一般可能仍然存在。即，在形成的第一组簇中可能存在模糊性。例如，在图4中，对于节点/读数402的对，可发现最大团发现进程对节点/读数402属于哪个团是模糊的。最大团发现仅仅是对各团成员的初始估计。因此，在该进程的阶段1结束时，一些序列读数116可能被分配至不正确的簇，并且一些应该在一起的序列读数116可能被分配至分开的簇。

图6是显示初始分配至错误簇112的序列读数的示例的图，其中相同填充特征的序列读数/节点来自相同同种型。如图所示，标记1-3的序列读数已被不正确地放在与序列读数4-5不同的簇中，其全部来自相同同种型。另外，序列读数10和11被不正确地与读数12分为一组，并且读数6已被不正确地与读数7-9的组分开。

再次参考图2，按照示例性实施方式，在初始聚类202之后进行的进程设计为解决初始簇122的模糊性。

簇共有区

在形成初始簇122之后，簇共有区算法111生成各初始簇的暂时簇共有区(块204)，其中各簇共有区124用于表示簇的所有成员的序列。簇共有区生成是本领域熟知的。例如，簇共有区算法111可基于使用有向非循环图来编码多个序列比对，例如，DAGCon(有向非循环图共有区)算法。考虑比对的读数300的集合，DAGCon取了一组成对的比对，其针对其他读数所比对的参照或主干/种子进行比对(基因组从头组装，最长的序列读数用作主干/种子)以生成有向非循环图，其中各条通过图的路径表示比对之一。该图然后简化并且确定最可能的通过图的路径，其为共有区。参见Chin等，来自长读数SMRT测序数据的非混合、精炼微生物基因组组装(Nonhybrid，finished microbial genome assemblies from long-readSMRT sequencing data)，Nature Methods(2013)，其通过引用纳入本文。

图7是显示分别针对各簇122生成的示例性簇共有区C1、C2、C3和C4，其中各输入序列读数116精确地属于一个簇122。

再次参考图2，在暂时簇共有区生成(块204)之后，援用错误校正(ICE)进程的迭代聚类的第二阶段。ICE的第二阶段由基于簇共有区124和读数质量值118迭代改进聚类开始(块206)。在该进程中，序列读数116从一个簇自动“再分配”到另一个簇，或者命名为“孤儿”并且用于生成新簇，如果序列读数被确定为不属于任意已有簇，并且合并高度相似的簇，如下所述。

图2显示了用于迭代改进聚类的进程(块206)的进一步详细说明。ICE组件112可通过使用质量值(QV)计算属于各簇(C)的各序列读数(S)的概率来开始迭代进程(块206-1)。这可通过将各簇122中的各序列读数116与各簇共有区C比对来完成。更具体地，各读数Si与各簇共有区Cu比对，其中“i”＝1至序列读数的总数，并且“u”＝1至簇共有区的总数(在图7所示的实施例中，i＝12并且u＝4)。

如果使用上述的检测同种型比对差异的进程，现有的序列读数(S)没有以足够高的相似性百分比与任意簇共有区(C)比对上(即，没有同种型命中)，则由于该序列具有差的概率而忽略该序列读数。在一个实施方式中，可使用线性时间算法来滤去具有较大插入缺失标记的比对。(参见，例如，用于在多重序列比对中定位极端保守的元件的算法(Algorithms tor locating extremely conserved elements in multiple sequencealignments)，Tseng和Tompa，BMC Bioinformatics，2009)。

如果现有的序列读数不与簇共有区中的一个或多个比对，则ICE组件112考虑现有读数的QV和簇共有区计算现有读数属于各簇的概率：

Pr(Si|Cu，QVs(Si))

如果QV不可得，则：

Pr(Si|Cu QVs(Si))＝(θ匹配)计数(匹配)(1/3θsub)计数(sub)(1/3θins)计数(ins)(1/3θdel)计数(del)，

其中θ分别是取代(sub)、插入(ins)和缺失(del)的匹配概率。

参考图7作为示例，当计算读数6的概率时，ICE组件112确定读数6属于簇共有区C3的概率大于读数6(S6)属于簇共有区C4的概率：

Pr(S6|C3)>Pr(S6|C4)

这可能是由于簇C3含有来自与S6同组的同种型。

概率计算的输出是对各读数序列属于各簇122计算的概率的列表。在一个实施方式中，计算的概率的数量是节点/序列读数的总数乘以簇的总数，一些概率具有“未知”的值。

再次参考图2，在计算概率之后，ICE组件112将个体序列读数从一个簇重分配至具有最高计算概率的另一个簇(块206-2)。

图8是显示将序列读数从一个簇重分配至具有最高的计算成员概率的簇的图。该实施例显示将读数6从簇C4重分配至C3。

应理解如果确定对现有序列读数有最高计算概率的簇就是该序列读数已经是成员的簇，则没有重分配。

再次参考图2，按照示例性实施方式的另一个方面，如果不存在比对(即任意序列读数与簇之间的概率未知)或如果线性时间算法排除了任意序列读数的所有比对，则序列读数可被视为孤儿，然后可从孤儿形成新簇(块206-3)。新簇可使用与上述初始阶段相同的过程从孤儿形成。

图9是显示从孤儿序列产生新簇的示例的图。在该实施例中，确定读数S12没有同种型命中。读数S12被称作孤儿并且产生含读数S12的新簇C6。

上述方法存在一个小问题，即当孤儿被分配至新簇，例如C6时，其仅有一个序列读数，该读数是其本身的代表。因此，具有一个序列读数的簇的计算概率将始终是1，这表示没有其他的簇会对该读数有更高的计算概率并且该读数不会重分配至另一个簇。仅具有一个序列读数的簇可被称为单现突变(singleton)并且没有成员以产生多样性，导致单现突变从来没有节点可重分配至的更好的簇，即便这样的一个簇可能存在。

按照一个实施方式，该问题可通过随机生成各孤儿节点的概率来解决。即，随机数生成器可用于生成预定范围，例如0-1的值(可能是其他范围)。如果随机概率小于预定阈值概率，例如，0.30，则孤儿被重分配至对孤儿的成员具有非零计算概率的簇之一。

再次参考图2，簇经处理以确定是否存在基本相同的簇，并且如果存在，将簇合并成新簇(块206-4)。基本相同的簇可在处理期间发生，由于大致最大团发现和迭代共有区判定进程。在一个实施方式中，基于它们的簇共有序列的相似性确定2个簇是否基本相同可能通过用户限定的参数，例如相似性百分比＝>95％来控制。

图10是显示2个簇合并的图。在该实施例中，来自图9的之前的簇C1至C4被确定为同种型命中并且具有大于99.5％的阈值相似性百分比。因此，在图10中，簇C1至C4已经合并成新簇和相应的簇共有区C7。

再次参考图2，每次簇的数量变化时，相应的簇共有区也可能变化。因此，ICE组件112更新各变化的簇的簇共有区并且更新所有序列读数的概率Pr(Si|Cu，QVs(Si))(块206-5)。这通过经线206-6判定簇共有区算法111(块204)来完成，并且因此产生错误校正进程的迭代聚类的第二阶段的“迭代环”，其第一步骤是重计算各序列读数的概率(块206-1)。

在一个实施方式中，簇共有区算法111可在每次发生变化时产生簇的簇共有区124。然而，在一个实施方式中，当基于簇大小判定簇共有区算法111即如果簇大小较大时，可任选地使用预定阈值来限制，如果特定簇中节点的数量大于预定阈值，判定簇共有区算法111可跳过块206-5。在某些实施方式中，簇共有区算法步骤是可并行的。

在一个实施方式中，新的额外序列读数可在第二阶段期间的任何时候通过将额外序列读数针对所有已有共有序列比对来增加至已有的簇组。如果已有簇C具有最高的在后概率并且比对不被拒绝，则新的序列分配至簇C。否则，该序列读数可如上述初始阶段那样被认为孤儿并且形成新簇。

一旦通过重分配序列读数和/或合并簇无法进一步改进簇，则用于错误校正进程的迭代聚类(块206)完成。

一旦用于错误校正进程的迭代聚类完成，则判定平滑组件114使共有区结果平滑(块208)。在一个实施方式中，平滑组件114可基于Quiver算法，如2013年7月12日提交的U.S.13/941,442中所述，其通过引用纳入本文。如上所述，ICE组件112判定簇共有区算法111生成各簇的簇共有区。在一个实施方式中，这些簇共有区用作全长读数116-1在迭代聚类进程期间所比对的“参照”。

按照示例性实施方式的一个方面，平滑(polishing)步骤的输入可包括簇共有区和非全长读数116-2，其然后比对至各簇共有区用作参照。在平滑期间，非全长读数116-2用于向序列读数施加全覆盖率以使用上述相同的“同种型命中”标准提供更高水平的共有区。在一个实施方式中，与全长输入序列不同，非全长读数116-2不必排他性地比对并且可属于多个簇。同样，使用线性时间算法来拒绝不利的比对。一旦非全长读数116-2与簇共有区比对，则平滑组件114生成各簇的最终簇共有区128(图1)。该进程的输出可包括最终簇共有区128的列表，其各自代表簇的“共有”序列。在一个实施方式中，各簇可用于表示单一、独特的转录本。

在另一个实施方式中，最终簇共有区128可映射至基因组，其中去除冗余并且同种型塌缩，从而生成高质量全长同种型。

按照示例性实施方式，用于错误校正进程的序列读数的迭代聚类可具有许多应用。例如，ICE可用于全长cDNA测序、生物信息学分析和生物学应用。

全长cDNA测序的示例可包括，但不限于，构建全长转录本富集的cDNA文库；使用琼脂糖凝胶或BluePippin^TM系统进行尺寸选择；对全长达到10kb的转录本进行测序；和对各转录本进行单分子观察。

除了同种型水平聚类以生成高质量转录本共有序列，生物信息学分析的示例可包括，但不限于，鉴定推定的全长转录本；和检测人工嵌合体。

最后，生物学应用的示例可包括，但不限于，新转录本；替代性剪接；替代性聚腺苷酸化；保留的内含子；融合基因；和反义转录。

在一些实施方式中，该系统包括可操作地耦合至处理器的计算机可读介质，其储存由处理器执行的指令。指令可包括下述的一种或多种：接收序列读数的输入的指令(和，任选的，参照序列信息)、构建预组装的读数的指令、比对序列读数的指令、生成字符串图像的指令、生成图像的指令、鉴定字符串束的指令、确定主要重叠群的指令、确定相关重叠群的指令、校正读数的指令、生成共有序列的指令、生成单倍型序列的指令、计算/储存与方法的各步骤相关的信息的指令(例如，字符串图像中的边缘和节点，字符串图像中的重叠和分支点、主要和相关重叠群)、和记录该方法结果的指令。

在某些方面中，该方法是计算机实施的方法。在某些方面中，算法和/或结果(例如，生成的共有序列)在计算机可读介质上储存，和/或在屏幕或打印纸张上显示。在某些方面中，结果经进一步分析，例如，以鉴定遗传变异，以鉴定序列信息的一个或多个来源，以鉴定个体或物种之间保守的基因组区域，以确定2个个体之间的相关性，以提供个体诊断或预后，或以提供健康护理专业人员可用于确定患者的合适治疗策略的信息。

此外，本发明的功能性方面在计算机或其他逻辑处理系统或电路上实施，如本领域普通技术人员所理解，可使用任意合适的实施环境或编程语言如C、C++、Cobol、Pascal、Java、Java-script、HTML、XML、dHTML、汇编或机器代码编程、RTL等来实施或完成。

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

在一些方面中，本发明包括用于多倍体基因组字符串图像汇编的制品，其包括含有一个或多个程序的机器可读介质，该程序在执行时实施本文所述的本发明的步骤。

应理解上述说明是示例性的而非限制性的。对本领域的普通技术人员而言，显而易见的是，可以对本发明进行各种修改而不会偏离本发明的范围和精神。因此，本发明的范围不应参照以上的说明决定，而应参照所附权利要求及其等同方案的全部范围决定。在本发明中，引用多个参考文献、专利、专利申请和公开。除非另外说明，各自出于所有目的通过引用纳入。出于说明和公开可与本发明关联使用的试剂、方法和概念的目的引用本文所有的出版物。本文并不旨在理解为承认这些参考文献相对于本文所述的发明是现有技术。

Claims (27)

1.一种用于错误校正的序列读数迭代聚类的方法，所述方法在至少一个处理器上执行的至少一个软件组件上进行，包括：

(A)接收序列读数的集和相关的质量值，其中所述相关的质量值是对每个位置碱基判定准确性的估计；

(B)将所述序列读数基于序列相似性分成初始簇的集；

(C)生成各所述初始簇的簇共有区；

(D)基于与所述序列读数相关的质量值和所述簇共有区迭代改进所述聚类，其中迭代改进所述聚类包括：(i)使用所述质量值计算各序列读数属于各簇的概率；(ii)将个体序列读数从一个簇重分配至具有最高计算概率的另一个簇；和(iii)合并高度相似的簇；

(E)生成并输出各所述簇的最终簇共有区。

2.如权利要求1所述的方法，其中使用所述质量值计算各序列读数属于各簇的概率还包括：

将各所述簇中的各序列读数与各所述簇共有区比对；

对于没有以足够高的相似性百分比与任意簇共有区比对上的现有序列读数，由于具有差的概率而忽略所述序列读数；

对于与簇共有区(C)中的一个或多个比对的现有序列读数(S)，考虑所述现有序列读数的质量值(QV)和簇共有区来计算所述现有序列读数属于各所述簇的概率：

Pr(Si|Cu，QVs(Si))。

3.如权利要求2所述的方法，其中对于QV不可得，则计算：

Pr(Si|Cu QVs(Si))＝(θ匹配)计数(匹配)(1/3θsub)计数(sub)(1/3θins)计数(ins)(1/3θdel)计数(del)，

其中θ分别是取代(sub)、插入(ins)和缺失(del)的匹配概率。

4.如权利要求1所述的方法，还包括：

对于任意所述序列读数和任意所述簇之间没有比对，认为所述序列读数是孤儿，并且从所述孤儿形成新簇。

5.如权利要求4所述的方法，还包括：

对于仅具有一个序列读数的新簇，对于各孤儿节点随机生成随机概率；并且对于随机概率小于预定阈值概率，则将所述孤儿重分配至对所述孤儿的成员具有非零计算概率的簇之一。

6.如权利要求1所述的方法，其中接收的序列读数包括全长的长读数，并且其中所述生成和输出最终簇共有区还包括：

将非全长读数输入最终平滑进程，其生成最终簇共有区。

7.如权利要求4所述的方法，其中所述读数包括长度范围为0.5kb至1、2、3、5、10、15、或20kb的全长的长读数。

8.如权利要求1所述的方法，其中基于序列相似性将所述序列读数分成初始簇的集还包括：

比对所述序列读数以产生比对的读数；使用所述比对的读数构建相似性图像；并且使用所述相似性图像发现最大团。

9.如权利要求8所述的方法，其中发现最大团包括：非确定性地将所述相似性图像划分成非重叠最大团。

10.一种在计算机可读介质上存储的可执行软件产品，所述产品含有用于错误校正的序列读数迭代聚类的程序指令，所述程序指令在至少一个处理器上执行，包括：

(A)接收序列读数的集和相关的质量值，其中所述相关的质量值是对每个位置碱基判定准确性的估计；

(B)将所述序列读数基于序列相似性分成初始簇的集；

(C)生成各所述初始簇的簇共有区；

(D)基于与所述序列读数相关的质量值和所述簇共有区迭代改进所述聚类，其中迭代改进所述聚类包括：(i)使用所述质量值计算各序列读数属于各簇的概率；(ii)将个体序列读数从一个簇重分配至具有最高计算概率的另一个簇；和(iii)合并高度相似的簇；

(E)生成并输出各所述簇的最终簇共有区。

11.如权利要求10所述的可执行软件产品，其中使用所述质量值计算各序列读数属于各簇的概率还包括：

将各所述簇中的各序列读数与各所述簇共有区比对；

对于没有以足够高的相似性百分比与任意簇共有区比对上的现有序列读数，因具有差的概率而忽略所述序列读数；

对于与簇共有区(C)中的一个或多个比对的现有序列读数(S)，考虑所述现有序列读数的质量值(QV)和簇共有区来计算所述现有序列读数属于各所述簇的概率：

Pr(Si|Cu，QVs(Si))。

12.如权利要求11所述的可执行软件产品，其中对于QV不可得，则计算：

Pr(Si|Cu QVs(Si))＝(θ匹配)计数(匹配)(1/3θsub)计数(sub)(1/3θins)计数(ins)(1/3θdel)计数(del)，

其中θ分别是取代(sub)、插入(ins)和缺失(del)的匹配概率。

13.如权利要求10所述的可执行软件产品，还包括：

对于任意所述序列读数和任意所述簇之间没有比对，认为所述序列读数是孤儿，并且从所述孤儿形成新簇。

14.如权利要求13所述的可执行软件产品，还包括：

对于仅具有一个序列读数的新簇，对于各孤儿节点随机生成随机概率；并且对于随机概率小于预定阈值概率，则将所述孤儿重分配至对所述孤儿的成员具有非零计算概率的簇之一。

15.如权利要求10所述的可执行软件产品，其中接收的序列读数包括全长的长读数，并且其中所述生成和输出最终簇共有区还包括：

将非全长读数输入最终平滑进程，其生成最终簇共有区。

16.如权利要求13所述的可执行软件产品，其中所述读数包括长度范围为0.5kb至1、2、3、5、10、15、或20kb的全长的长读数。

17.如权利要求10所述的可执行软件产品，其中基于序列相似性将所述序列读数分成初始簇的集还包括：

比对所述序列读数以产生比对的读数；使用所述比对的读数构建相似性图像；并且使用所述相似性图像发现最大团。

18.如权利要求17所述的可执行软件产品，其中发现最大团包括：非确定性地将所述相似性图像划分成非重叠最大团。

19.一种用于错误校正的序列读数迭代聚类的系统，包含：

存储器；和

与所述存储器耦合的处理器，其设置成：

(A)接收序列读数的集和相关的质量值，其中所述相关的质量值是对每个位置碱基判定准确性的估计；

(B)将所述序列读数基于序列相似性分成初始簇的集；

(C)生成各所述初始簇的簇共有区；

(D)基于与所述序列读数相关的质量值和所述簇共有区迭代改进所述聚类，其中迭代改进所述聚类包括：(i)使用所述质量值计算各序列读数属于各簇的概率；(ii)将个体序列读数从一个簇重分配至具有最高计算概率的另一个簇；和(iii)合并高度相似的簇；

(E)生成并输出各所述簇的最终簇共有区。

20.如权利要求19所述的系统，其中使用所述质量值计算各序列读数属于各簇的概率还包括：

将各所述簇中的各序列读数与各所述簇共有区比对；

对于没有以足够高的相似性百分比与任意簇共有区比对上的现有序列读数，因具有差的概率而忽略所述序列读数；

对于与簇共有区(C)中的一个或多个比对的现有序列读数(S)，考虑所述现有序列读数的质量值(QV)和簇共有区来计算所述现有序列读数属于各所述簇的概率：

Pr(Si|Cu，QVs(Si))。

21.如权利要求20所述的系统，其中对于QV不可得，则计算：

Pr(Si|Cu QVs(Si))＝(θ匹配)计数(匹配)(1/3θsub)计数(sub)(1/3θins)计数(ins)(1/3θdel)计数(del)，

其中θ分别是取代(sub)、插入(ins)和缺失(del)的匹配概率。

22.如权利要求19所述的系统，还包括：

对于任意所述序列读数和任意所述簇之间没有比对，认为所述序列读数是孤儿，并且从所述孤儿形成新簇。

23.如权利要求22所述的系统，还包括：

对于仅具有一个序列读数的新簇，对于各孤儿节点随机生成随机概率；并且对于随机概率小于预定阈值概率，则将所述孤儿重分配至对所述孤儿的成员具有非零计算概率的簇之一。

24.如权利要求19所述的系统，其中所述接收的序列读数包括全长的长读数，并且其中所述生成和输出最终簇共有区还包括：

将非全长读数输入最终平滑进程，其生成最终簇共有区。

25.如权利要求22所述的系统，其中所述读数包括长度范围为0.5kb至1、2、3、5、10、15、或20kb的全长的长读数。

26.如权利要求19所述的系统，其中基于序列相似性将所述序列读数分成初始簇的集还包括：

比对所述序列读数以产生比对的读数；使用所述比对的读数构建相似性图像；并且使用所述相似性图像发现最大团。

27.如权利要求26所述的系统，其中发现最大团包括：非确定性地将所述相似性图像划分成非重叠最大团。