Pasta : 超大数据多序列比对
摘要:在这篇文章中,我们介绍一个新的、高可扩展性的算法:Pasta,用来估计大规模的序列比对。论文用了真实的和模拟的数据进行测试,初测试结果比SATe好一些,比其他方法有显著提高。(其实整个程序就是基于SATe开发的)
1、简介
多序列比对是许多生物信息分析的基本步骤,包括预测RNA和蛋白质的结构和功能,还有评估进化关系。但是只有一小部分的多序列比对方法可以分析超过10000条序列的大规模数据集。在大规模数据集上的多序列比对方法性能评估显示,一些多序列方法在用传统的比对标准评估时可以产生高准确率的比对。然而,专注于对核苷酸数据库进化关系的评估发现,只有一小部分的多序列比对方法在面对10000条以上的核苷酸序列数据库时,可以提供好的比对,从而产生高准确率的进化树。然而这些研究依靠的标准测试是最多28000条序列,因此关于大规模数据集的序列比对准确性及其对进化树构建的准确性所知甚少。目前对于超过100000核苷酸序列的系统进化关系分析有两个组织研究,iPTOL项目以及Thousand Transcriptome 项目。
在这篇文章中,展示了PASTA,practical alignments using sate
and transitivity (利用SATe和传递性的使用序列比对方法), 用来处理极大规模的多序列比对。PASTA首先对序列和树的估计是使用了一个非常简单的基于隐马尔科夫技术(HMM),然后根据进化树对序列进行了二次比对。如果需要可以根据新的序列比对重新构造一棵进化树,然后算法可以迭代。
PASTA的准确性和可扩展性的关键是采用了新的技术用来估计序列比对和进化树。如同SATe一样,PASTA利用了centroid edge 数据库分解技术把序列集合分成一个个子集合,然后在子集合上使用了MAFFT-linsi进行序列比对,但是PASTA和SATE在合并这些子集合比对成为一个完整比对集合的时候使用了不同的方法。SATE用了Opal来合并子集合的比对,PASTA只利用Opal合并一对对邻近子集合以及重叠的比对,然后视每组比对为等价关系,然后利用传递性来合并这些大规模比对。这样的再比对(re-alignment)方法是高度并行的,同时可以扩展到大规模数据集上。同时,这个再比对过程对于PASTA来说是一个微不足道的部分,然而对于SATE来说是一个主要占用运行时间的部分(10000条序列占44%的运行时间,50000条序列战78%的时间)。这样PASTA在大数据处理上相对SATE来说有显著的提升。同时,PASTA 比SATE产生了更高准确率的比对和进化树。将在收集到的数据计算PASTA的速度以及准确率,这个数据集是200,000条RNAsim数据集,在这个数据集上的比对时间小于24小时,程序运行在12核的机器上。
PASTA
PASTA对数据的序列集合S利用了迭代分治的策略,同时使用了下面的输入参数:一个初始树,子集合的大小k(默认200),子集合比对技术(默认是MAFFT-linsi),子集合合并技术(OPAL),停止规则(默认是3轮迭代)。
第一轮迭代由初始树开始,上一轮的迭代输出作为下一轮迭代的输入,每一轮的迭代包含六个步骤:
(1)分解输入集合S成为若干个子集合,S1,S2,S3,…Sm,子集合序列最多k条。
(2)计算生成树T*,来连接子集合S1……Sm。
(3)利用子集合比对技术来比对每个子集合。
(4)将生成树T*的每条边上的端点合并。
(5)利用传递闭包的连续性合并在第4步中的重叠和兼容序列。
(6)利用FastTree-2在整个多序列比对上计算最大似然树(ML tree)
初始树的构建方法: PASTA可以由任意合理的初始树开始。没有初始树的话就从S中任意抽取100条序列,然后再从100条序列中抽取子集合X,然后在这个子集合上计算SATE比对A;这部分叫做“主比对”。我们然后用HMMER在A上来计算隐马尔科夫模型,然后依次将S-X 与A做比对,至此就创建了整个数据集的比对。 我们然后利用FastTree-2在这些集合上构建最大似然树(ML tree)。(这样就有了初始比对和初始树)
步骤1: 利用SATE的质心分解技术在现有的指导树(guide tree,我觉得应该就是上面构造的初始树)上将序列集合分解成不相交的集合,S1,S2…Sm,每个子集合最多k条序列。如果树最多只有k个叶节点,那我们就返回序列集合(就是说数据集比较小,整个数据集都不到k条序列);否则的话,我们就找到树的边缘将集合分裂为大小大致相同的两个子集合,然后在每个子集合上重复分裂。
步骤2:我们对子集合构造生成树T*。(下文是如何构建生成树)对于集合Si上的两个叶子节点之间的路径上的点全都标成Si。然后假如还有节点i有被标注,就找到这个节点最近标注节点,直到所有节点都被标注。然后把端点上的标注是相同的边都删除,这样就生成了S1,S2,…Sm上生成树。
步骤3:对每个子集合Si 计算多序列比对(Type 1 sub-alignment)。
步骤4:生成树T*的每个节点都是一个子集合,每个子集合都由Type1 比对。对于生成树T*上的每条边(v,w),用特定的对齐合并技术合并v 和w,这产生了新的比对集合,每个包含最多2k的序列,称为(Type 2 sub alignmnets),要求Type2比对不能改变Type1的比对,
步骤5:根据传递性原理计算合并。(没看懂)
步骤6:运行FastTree-2估计出最大似然树。
为了加速计算,掩盖那些插入空格超过99.9%的位点。注意到PASTA的比对合并步骤在引进新的同源上是比较保守的,这样PASTA就趋向于产生许多插入多个空格的空隙。掩盖这些插入空格多的空隙对树的估计步骤是没有损害的,但是对运行时间会有很大的影响。
运行时间 :PASTA把比对以压缩的格式保存在内存,压缩的格式是保存原始的序列,以及序列中每个字符在最终序列里的位置(‘ACCA’,[1,3,5,6],‘A-C-CA’)这种格式减低了内存的需求,同时减小了传递合并的运行时间。
只要每对传递合并是正确执行,那最终的输出多序列比对就不取决于是从生成树的那条边开始执行,但是执行顺序会影响运行时间。但是如果合并sub-alignments用的是逆序的质心分解,那么运行时间就是有界的:
定理1: 给定m个Type1 alignments 和 m-1个Type2 alignments 那么算法用来计算传递合并花费的时间是O(Kmlog(m)+Lm),其中K 是任何序列的最大长度,L是输出的比对长度。任意的输入时间复杂度是O(Km^2+Lm)
实验设置
数据集: 为了评估一个中型数据集的准确率,用1000-taxon 模拟数据库。为了评估大型数据库的性能,用RNAsim,一个拥有1000000条序列的模拟RNA的数据库,二次抽样生成10000,50000 100000,200000 条序列。http://www.cs.utexas.edu/users/phylo/software/pasta/.对于10000的样例,创造10个不同的副本,但是对于其他样例由于运行时间的需求,只进行一个复制。最后用16s biological数据库Gutell lab ,以及1000-taxon(数据已经下载)。 这些数据库包括16S.3,6323条序列,16S.T,7350条序列,16S.B.ALL ,27643条序列。生物数据库的参考比对是根据二级结构,参考树的计算是使用了RAxML。
方法: 比较了PASTA ,SATE, Muscle, MAFFT-Profile ,ClustalW 。 PASTA结果是基于默认设置的:三轮迭代,子集合大小为200,MAFFT-linsi用于子集合,Opal 用于合并序列。初始树由之前的section2 描述的方法构建。MAFFT-Profile是MAFFT的一个可以添加新序列到已经存在的主干比对的新版本。实验提供给MAFFT-Profile与PASTA的初始树相同的主比对。运行SATE 的时候与PASTA执行完全一致的初始树,在SATE上同样是执行三轮迭代。由于在大规模数据集上OPAL的时间花销太大了,当序列条数超过了5000条时候,用Muscle来合并SATE的比对。最后用FastTree-2来计算最大似然树。
标准:测试了比对的准确率,树的错误,运行时间。
比对的准确率使用的是FastSp ,根据两种不同的度量:SP—score (参考序列在估计比对中同源序列的百分比) 以及modeler score(同源序列占正确的估计比对的百分比)然后平均一下,获得一个标准。注意SP -score是SP-FN 错误率的补,同时modeler-score是SP-FP错误的补,这样预测比对的准确率同时受到false positive 和false negative的影响。
除此之外,还统计了在估计比对中列恢复正确(the number of columns
that are recovered entirely correctly)的数量(TCscore)。
对树的评估的标准是bipartition distance , 也叫 Robinson - Foulds(RF)rate , 但是当参考树无法构建(参考树的计算是使用了RAxML)的时候这个度量是不合适的。因此用False Negative(FN )rate ,这个是统计估计树中在真实树上分错的边树。值得注意的是,当估计树和参考树是二分的时候FN值和RF值是一样的。
计算平台:所有的方法执行在Lonestar Linux 集群(TACC), 每台机器是24G内存12核。