生物信息标识Link to Publisher's site
生物信息学. 2015年6月15日;31(12):1913–1919年。
2015年1月31日在线发布。 内政部:10.1093/生物信息学/btv053
PMCID:PMC4765884号
PMID编号:25638815

星码:基于全对搜索的序列聚类

摘要

动机:测序技术的不断增长为计算生物学提供了新的应用和挑战。在许多这样的应用程序中,排序错误需要纠正。当测序从一个未知的参照物(如随机的DNA条形码)读取时,这一点尤其重要。在这种情况下,可以通过对所有条形码进行成对比较来进行纠错,这是一个计算复杂的问题。

结果:在这里,我们解决了这个难题,并描述了一个精确的算法来确定在给定的Levenshtein距离内的序列对。为了纠错或减少冗余,匹配对然后被合并成相似序列的簇。星码的效率归功于poucet搜索,这是一种在trie节点上执行的Needleman-Wunsch算法的新实现。在随机条码匹配方面,星码在速度和精度上都优于序列聚类算法。

可用性和实施:C源代码在http://github.com/gui11aume/starcode.

联系人: moc.liamg@noilif.emualliug

1简介

所有测序技术都有一定程度的不精确性。例如,Illumina平台(马古利斯et公司 艾尔。2005年)有1-2%的错误率由替换组成(多姆et公司 艾尔。2008年;中村et公司 艾尔。2011年)PacBio平台有15%的错误率,包括插入和删除(开斋节et公司 艾尔。2009年). 这类技术的巨大吞吐量最近为开发高效的纠错算法创造了额外的需求。

测序错误可以通过将读数与参考基因组进行比较来发现。然而,这种参考并不总是可用的。当策略是从一个未知的或随机的序列中选取一个错误来纠正。例如,当出现这种情况时,使用条形码来跟踪(阿克塔尔2013年;舍佩尔斯et公司 艾尔,2008年). 排序错误将产生错误(不存在)的条形码,必须将其删除。

序列聚类可以看作是图上的一个群体检测问题,其中节点代表序列,边表示相关序列之间的匹配。该过程包括一个匹配阶段(计算量最大的阶段),在该阶段构建图,并在聚类阶段识别社区。

在这里,我们描述了一种序列聚类算法,称为“星码”,参考随机条码的簇,通常是星形的。星码是基于所有对搜索的,即在图的构造阶段,在给定的Levenshtein距离下的所有序列对都被识别出来。匹配通过无损过滤进行,然后对前缀trie的分支进行彻底搜索。该算法的新颖之处在于poucet策略,它利用了按字母顺序排序的序列的冗余性,避免了不必要的重新计算,提高了速度。

在这篇文章中,我们展示和基准星代码。我们表明,在真实的生物数据集上,starcode比现有的序列聚类软件快一个数量级。即使星码是为纠错而设计的,我们也证明了它也可以用于其他问题。作为一个例子,我们用它来识别细菌基因组和蛋白质-RNA相互作用实验中丰富的基序。

2种方法

2.1使用tries的不精确字符串匹配

星码匹配方法是基于Needleman-Wunsch(NW)算法的一种改进(奈德曼和温什,1970年). 在原始算法中(图1a) 两个序列之间的Levenshtein距离是通过在一个术语(编辑矩阵),其中n是各自的序列长度。这种动态规划方法的复杂性是O().

保存图片、插图等的外部文件名为btv053f1p.jpg

NW序列对比。()GTTGCA与GATCCA的比较。利用NW动态规划算法初始化编辑矩阵的边距,从左到右、从上到下计算单元。E[j],坐标项(,j)计算为最小值⁡(E[−1个j]+1E[j−1]+1E[−1个j−1]+Δ(j)),其中Δ(j)=0如果与第一个符号的序列相同j从第二个符号Δ(j)=1否则。两个序列之间的Levenshtein距离是右下角单元格的值。(b)确定GTTGCA和GATCCA是否为2-匹配的低复杂度算法。外部单元格中的值是在初始化期间设置的。动态规划算法如上所述,不同之处在于,如果对角线单元格(粗体边框)的值大于2,则会中止该算法。初始化单元中的值可能与原始NW方案(箭头)不同,但计算单元中的值仍然相同。空单元格的值从不计算,这有助于降低复杂性

在许多情况下,感兴趣的信息是找出序列是否是τ-匹配(即它们的距离小于或等于固定的阈值τ)。在这种情况下,复杂性可以降低到O(τ最小值⁡(n)). 它不是计算编辑矩阵的所有项,而是初始化为,如所示图1只计算对角线周围的项。如果一个对角项的值大于τ,则该过程将停止,因为序列不是τ-匹配的。

此方法可用于根据前缀树(也称为trie)匹配序列(科努肯,1995年). 编辑矩阵的项按行更新,而深度优先搜索遍历trie(图2). 每次访问节点时,计算一行,每次搜索回溯时,都会删除一行。如果对角项超过阈值τ,则所有下游序列的Levenshtein距离也必然大于τ。因此,不再在该路径中发现更多的命中,深度优先搜索将回溯到父节点。当进程停止时,该搜索路径上的每个尾节点(对应于数据库的一个序列)都是查询的τ-匹配。这种方法是有效的,因为它消除了大面积的搜索空间,而且查询与数据库每个前缀的NW比较只计算一次。

保存图片、插图等的外部文件名为btv053f2p.jpg

NW算法在尝试。索引的每个序列都是trie中的一条路径。查询GTTGCA写在矩阵的顶部,矩阵初始化如所示图1b、 trie通过深度优先搜索(路径以箭头结尾)遍历。在每个深度处,添加到路径的节点被写入编辑矩阵的左侧,并计算该行。从1到4的检查点(带圆圈的数字)显示搜索过程中编辑矩阵的状态。因此,叶对应于一个标记为-3的查询。发现命中后,搜索路径返回到标记为2的节点,编辑矩阵的最后一行被删除。然后搜索路径转到标记为4的节点,在这种情况下,新计算的对角线单元格超过阈值(带圆圈)。即使此节点有子节点,也不会访问它们(交叉),因为没有要发现的2-match

2.2 poucet搜索算法

可以进一步改进搜索策略。如果两个连续查询共享一个长度前缀k,直到k对于两个查询,编辑矩阵的第行将完全相同。因此,计算中间产物可以存储在trie的节点中,以便下一个trie搜索可以从深度开始k. 然而,在一些编辑节点中存储矩阵行遇到了困难。的确,在k第二行,对角线右侧的术语依赖于两个查询之间不共享的字符。这个问题的解决方法是在每个节点中存储行和列项的组合,它们形成一个角度形状,看起来像一个水平翻转的L(图3). 使用这种结构,计算的中间产物存储在一个节点的深度k只依赖第一个k查询的字符。

保存图片、插图等的外部文件名为btv053f3p.jpg

Poucet搜索算法。该算法的原理与上所示的相同图2不同的是,编辑矩阵不是按行更新的,而是沿着水平翻转的L进行的。随着深度优先搜索的进行,这些值存储在trie的节点中。由于翻转L的垂直部分中的值对于节点的每个子节点都是相同的,因此只计算一次(箭头)。当搜索路径访问节点时,将计算灰色单元格中的值。将中间层存储在节点中允许下一个查询在深度处重新启动k如果它共享一个共同的长度前缀k使用当前查询

最大限度地利用连续查询之间的字母顺序进行排序,从而最大限度地利用连续输入的前缀。在童话故事“Le Petit Poucet”中,主人公用白色的鹅卵石为哥哥们寻找回家的路,这让人想起一个较小的查询(按字母顺序)为下一个查询铺路的方式。因此,我们称这种搜索算法为“poucet”。

2.3无损过滤

当查询没有匹配项时,最好省略trie搜索。为此,starcode使用了一种类似于吴和曼伯(1992). 查询最初分为τ + 1段。假设所有段的长度至少为τ,那么数据库中的每个τ-匹配将至少包含一个查询段的逐字副本。实际上,查询和匹配之间最多有τ个版本,分布在τ + 1区域中,因此至少有一个片段是未修改的。由于前面片段中可能存在插入和删除,共享片段可能会从查询的原始位置向上移动到左侧(所有插入)或右侧(所有删除)的τ核苷酸。

这些观察结果是具有100%灵敏度的过滤方法的基础。更精确地说,这些片段的定义如下:去掉序列的第一个τ核苷酸,序列的其余部分被划分为大小不超过1的τ + 1片段(为了一致性,较长的片段总是在3'范围内)。每次将一个序列添加到trie中,它都会被分区,并且它的段被添加到τ + 1不同的索引中。第一个片段被添加到第一个索引中,第二个片段添加到第二个索引中,等等。在搜索之前,按相同的方式对查询进行分区,并在索引中查找其片段。如果没有找到匹配项,则该查询在当前数据库中没有τ-match,因此省略trie搜索。如果至少要找到一个trie,则至少要找到一个trie。

如上所述,在查询和τ-匹配之间共享的片段可以被发现上移到τ核苷酸上。因此,根据图4,确保不遗漏匹配:在τ + 1索引中查找最右边的片段,在τth索引中查找第二个最右边的片段和被一个核苷酸移位的相邻片段,依此类推,直到在第一个索引中查找到第一个片段及其移位最多τ个核苷酸的相邻片段。

保存图片、插图等的外部文件名为btv053f4p.jpg

用̴长度为̴的例子说明̴长度为̴的序列。删除查询的最后一个τ核苷酸,其余的被分成四个连续的片段。每个序列都是针对编号为I–IV的不同索引进行查询的。例如,针对索引I查询的唯一片段是GTTG,而针对索引II查询的片段是GCAA、CAAT和AATA。如果在适当的索引中找到任何一个片段,则执行trie搜索,否则将忽略它,因为不可能存在τ-匹配。不管结果如何,标记为I–IV的段都会添加到相应的索引中(即每个索引中只添加一个片段)

2.4寻找和建造

同时建立了搜索空间的动态搜索方法。trie在单词序列中相互匹配之前插入。如果A和B是相互τ-匹配,则当B在trie中时将查询A,反之亦然。不管怎样,匹配A-B都会被发现。这保证了每个τ-匹配都被发现,同时保持trie尽可能“薄”,从而减少搜索时间。算法1和算法2所示的伪码概括了整个匹配过程。

算法1星码算法-

1号文件:定义:τ

第二章:变量: 种子,开始 = 0,高度,顺序,特里亚,最后一个序列,k

三:容器: 击打,卵石

第四章:阅读序列文件

第五章:heght ← 决定最大序列长度

第六章:衬垫序列高度

第七章:分类按字母顺序排列

第八章:k ← 计算过滤器段长度

第九章:tre ← 创造空的高度高度

十:插入的根节点特里亚在里面卵石深度0

十一:对所有人序列

十二:se ← 得到下一个序列

十三:如果至少一个k-的mer顺序在筛选器索引中然后

14日:seed ← 长度当前序列和下一序列之间共享前缀的

十五:strt ← 长度共享前缀顺序最后一个序列

十六:清楚的击打

17日:清楚的 卵石深度>strt

18日:对所有人 卵石在深处开始

19日:node ← 得到下一个节点来自卵石

二十:呼叫 波塞特(顺序,节点,种子,击打,卵石)

21日:结束

22日:过程 击打链接匹配到顺序

23日:stse ← se

24日:结束if

25日:插入 顺序进入特里亚

26日:插入 顺序k-mers进入过滤指数

27日:结束

算法2 Poucet搜索算法-

1号文件:程序 波塞特(查询,节点,种子,击打,卵石):

第二章:计算 节点-NW◃后的特定列图1

三::对所有人 小孩中的节点节点

第四章:计算 小孩-NW◃后的特定行图1

第五章:5:计算使用行和列的中心值◃图1

保存图片、插图等的外部文件名为btv053f5p.jpg

可扩展性。()运行时间的对数与要聚类的序列数的对数。(b)运行时间是聚类距离的函数。(c)运行时间和输入序列的长度。(d)不同并行线程数的相对性能提高

第六章:6:如果中心值>τ 然后◃超出不匹配。

保存图片、插图等的外部文件名为btv053f6p.jpg

已知聚类结构的人工数据集的基准测试结果(见正文)。()准确度由识别出的簇的数量来衡量。Starcode可以识别正确数量的簇,而seed和cd-hit则可以识别每个真阳性大约40个误报。绘制出第一条平分线,并显示出完美的结果。(b)通过识别对的数量来测量精度。Slidesort识别的配对比starcode少5-10%。水平线表示比率为1。(c)不同工具的运行时间。随着星团的大小,星码的运行时间增加,但仍具有竞争力。(d)不同工具的内存使用情况。星码的内存使用量随着星团大小的增加而减少。

第七章:持续与下一个小孩

第八章:结束 如果

第九章:如果 节点深度 = 高度 然后◃发现命中。

十:节约 节点顺序击打

十一:持续与下一个小孩

十二:结束 如果

十三:如果 节点深度≤seed 然后

14日:节约 节点在里面卵石在当前深度

十五:结束 如果

十六:呼叫波塞特(查询,小孩,种子,击打,卵石)

17日:结束 对于

18日:结束 程序

2.5并行化

查询在连续的块中进行排序和分区。然后匹配步骤分两个阶段进行。在构建阶段,根据上述算法从每个块的序列中构建一个不同的trie。在第二种方法中,所有序列块都将根据所有其他尝试进行查询。如果查询被分区到N第一阶段包括N寻找和建设工作,而第二个工作包括N(N−1)/2查询作业。在每个阶段中,作业之间不存在依赖关系,因此可以有效地并行化匹配算法N大于独立线程数。

2.6聚类

星码的默认聚类算法是为了纠正排序错误而设计的。此方法使用消息传递(麦凯,2002年)识别和计算“典型”序列(在聚类术语中也称为质心)。默认情况下,每个序列将其读取计数转移到最接近的τ-匹配,前提是后者的计数至少多出5倍。如果不满足条件,则不进行传输。如果这些序列中有数个相等的相合数τ相等,则它们之间的相合数相等。从读取计数最少的序列开始,递归地重复该过程。进程结束时读计数为正的序列被认为是规范的。集群由所有将其读取计数传输到同一规范序列的序列组成(将丢弃将其读取计数传输到不同规范序列的序列)。请注意,簇的半径可以大于用于匹配的最大距离。

由于没有一种测序技术的错误率高于20%,因此,从测序错误中出现的序列的读取计数通常是标准序列的5倍或更低。否则,序列更有可能是不相关的,或者两者都是从同一个规范序列派生的。可以使用命令行选项修改此行为簇合比允许更灵活或更严格的聚类,例如将唯一的输入序列聚类在一起,簇合比必须设置为1。

对于其他序列聚类问题,starcode实现了一种称为“球体聚类”的多用途算法。在球面聚类中,序列按出现频率排序。从最频繁的序列开始,每个序列都成为规范序列并声明其所有的τ-匹配,从而形成半径τ的簇(因此得名)。它们只能立即属于一个序列。

7.2基准条件

所有测试都是在一个16核双处理器Intel Xeon E5-2687W v2系统上进行的,该系统具有256 GB DDR3-RAM,频率为1866兆赫。在所有软件中,命令行参数都被等效地设置为以单核模式运行,允许长度为50的输入序列最多有三个不匹配。表1和22分别总结仿真和真实数据集中使用的执行选项。

表1。

仿真基准测试中使用的软件执行选项

软件命令行选项
星码-1.0星码 -d3
幻灯片-2幻灯片演示2 -d -t E -c DNA
Cd-hit-est-4.6.1cd命中测试 -n 9 -c 0.9 -米 0 -r 0
种子-1.4.1种子 –不匹配

表2。

实际数据基准测试中使用的软件执行选项

软件命令行选项
星码-1.0星码 -d3
幻灯片-2幻灯片演示2 -d -美国 -t E -c DNA
Cd-hit-est-4.6.1cd命中测试 -n 8 -c 0.94 -米 0
种子-1.4.1种子 不匹配 –轮班
彩虹-2.0.3彩虹 集群 -米

3结果

3.1介绍和基本性能

Starcode是一个通用的DNA序列聚类工具,主要用于纠错。错误被假定为不匹配、插入或删除(这里提供的实现匹配最多有8个错误的序列)。输入序列可以是单端或成对端读,上限为1024个核苷酸(配对端为512个)。序列可以是可变长度的,它们可以被修剪和过滤以保证质量。与starcode兼容的文件格式有raw sequence、raw sequence with count、FASTA或FASTQ(在这种情况下,starcode忽略质量)。Starcode要么返回聚类结果的详细信息,即规范序列、簇大小及其组成序列的完整列表。或者,只打印规范序列,这对于从输入文件中过滤出冗余序列非常有用。默认情况下,聚类是在假设实验误差(测序错误、聚合酶链反应突变等)发生偏差的情况下进行的,对于其他聚类问题也可以使用更通用的算法(示例见第3.3节)。

我们展示了星码在伪随机序列数据集上的基本性能和可伸缩性(图5). 标准配置是一组长度为40的1 000 000序列,运行在一个线程上,最大Levenshtein距离为3。在每个测试中,只有一个参数被修改,而其他参数保持不变。由于聚类步骤不需要额外的内存分配,而且比所有对搜索都快得多,第3.1节和第3.2节中给出的性能结果同时适用于消息传递和球体聚类算法。

图5a显示了星码的运行时间作为输入序列数的函数n. 在双对数尺度下,趋势是一条斜率为1.5的直线,表明星码的运行时间复杂度低于二次型(全对搜索的朴素实现)。请注意,此数据集的序列不匹配,有关更真实数据集的性能评估,请参阅第3.2节。图5b表明,运行时间随用于聚类的最大Levenshtein距离呈指数增长。原因是,随着最大距离的增加,trie呈指数级散开,搜索在更大深度退出。作为序列长度的函数,运行时间先增加,然后急剧下降并保持在较低水平(图5c) 一。超过阈值长度后,过滤算法开始变得有效,并且大多数查询都在不搜索trie的情况下得到解决。最后,我们展示了随着图5d、 搜索算法是完全并行的,相对性能线性增加到12个线程。随后观察到的弯曲有两个原因:第一个是输入读取和聚类步骤很短但不是并行的;第二个是由于硬件限制,即没有足够的内存带宽来满足增加的并行内存访问需求。

3.2基准

我们根据序列聚类算法slidesort对starcode进行了基准测试(清水和筑田,2011年),种子(et公司 艾尔。2011年),彩虹(et公司 艾尔。2012年)以及cd播放(et公司 艾尔。2012年). 尽管slidesort是一种全对搜索算法,但它也被包含在基准测试中,因为序列比较是序列聚类问题中计算量最大的步骤。Rainbow只在成对的结束读取上运行,而其他工具只在单个读取上运行,因此所有工具不能在同一个数据集中运行。

序列聚类算法的性能对数据集中的簇的大小非常敏感,这在许多应用中是先验的。因此我们在已知的聚类分析工具上建立了一个精度测试工具。我们生成了四个数据集,这些数据集由1到1000个簇组成,共有100万50个mer。每一个星团由100个相同质心序列的重复序列组成,再加上通过合并三个误差(包括最多一个索引)从质心导出的卫星。每个集群的卫星数量从999 900到900。Rainbow没有在这个基准测试中测试,因为生成的数据是单次读取。此外,评估slidesort的准确性是有问题的,因为每个数据集中的3个匹配的数量是未知的(同一个集群中的成对卫星可能是3个匹配或不匹配)。出于这个原因,我们只比较了starcode找到的对的数量与slidesort找到的对的数量。测试结果总结如下图6.

虽然starcode在所有四个数据集上都实现了完美的聚类,但是通过seed和cd-hit实现的聚类是不完整的。这两个工具在所有数据集上识别每个真集群大约40个假集群(图6a) 一。我们还观察到slidesort在所有数据集上发现3个匹配比starcode少5-10%(图6b) 一。我们对这个结果感到惊讶,因为slidesort声称是一个精确的算法。然而,当我们在较小的数据集上运行额外的测试时,情况显然不是这样的,在这些数据集上,朴素的成对比较是可行的(更多信息可从starcode存储库获得,参见http://github.com/gui11aume/starcode/tree/master/misc).

不同工具的运行时间与集群大小的函数如所示图6c、 slidesort和starcode的运行时间分别呈现线性和次线性趋势。无论集群大小,Seed和cd hit大约在恒定时间内运行。尽管如此,星码的性能仍然具有竞争力,即使是对于100万个序列的簇。内存使用情况如所示图6c、 slidesort和cd-hit实现了最小的内存占用,与其他工具相比,最大的差异是一个数量级。请注意,与slidesort的比较并不完全公平,因为它没有在内存中保存集群所需的完整图形。星码的内存使用对于1000大小的集群是最高的,但是随着集群大小的增加,它会降低并低于seed的内存使用量。总之,starcode是唯一一个以增加内存占用为代价实现这些数据集精确性的工具。考虑到输出的精确性,starcode在运行时间方面保持了竞争力。

人工数据的性能并不总是与实验数据集的性能一致。典型的实验带来了额外的困难。例如,簇的大小可能不均匀,并且读取可能包含几乎恒定的区域,这些区域通常会降低基于滤波器的算法的性能。我们对序列聚类算法进行了基准测试,以解决将数千名记者集成在并行(TRIP)条形码中的问题(阿克塔尔et公司 艾尔。2013年). 简单地说,TRIP的原理是用随机条形码标记报告者转录本,并测量RNA中条形码的丰度作为基因表达的代理。因为不存在与条形码匹配的未知条形码。

中总结了用于基准测试的数据集的基本属性表3. 数据集1(SRR950457)已进行预处理,以提取条形码并删除读取的常量部分。文件中只有15到17个核苷酸之间的条形码。数据集2(PRJEB7686)由原始Illumina单次读取组成。这些数据集因读取大小、总读取计数和经验簇大小而不同。根据starcode的输出,数据集1的最大簇包含约70 000个序列,而数据集2包含4个超过100万个序列的簇。数据集3(SRR950477)已经被包括在成对端聚类模式下的starcode与rainbow的基准测试。

表3。

用于基准测试的生物数据集摘要

数据集读取计数读取长度类型
SRR9504576 542 30916± 1单身
PRJEB7686型127 675 53750单身
SRR9504772 460 226100 + 100成对端

所有数据集都是Illumina reads。

星码、种子、slidesort、rainbow和cd-hit的运行时间总结如下表4. 我们调整了第一个数据集的距离阈值,以补偿减少的序列长度。starcode和slidesort都是用选项'-d 2'并且cd hit的标识设置为'-c 0.85’. 由于最小序列长度的限制,我们无法在数据集1上运行seed。在所有数据集上,Starcode比其他工具都快得多。Seed和cd-hit排在第二位,在数据集1和2上的运行时间分别是前者的35倍和20倍。Rainbow在聚类成对结束读取的工作中慢了近一个数量级。我们没有记录过去10天的确切运行时间,因为这比星码的运行时间高出几个数量级。

表4。

软件在三个生物数据集上的运行时间(秒)

软件SRR950457PRJEB7686型SRR950477
星码5289844
种子60 374
幻灯片4055>10天
彩虹306
Cd命中测试170512 591

在过去的10天里,没有记录确切的运行时间。破折号表示软件不能用于此数据集。

相同数据集上不同工具的内存占用如所示表5. 这些值表示在上述数据集上运行期间内存使用的峰值。在短读(数据集1)中,starcode利用trie压缩的优势优于其他工具。在数据集3上,starcode的内存使用量比rainbow大得多。Starcode和cd-hit在数据集2上使用了相似的内存量。两者所需内存是slidesort的两倍,slidesort的优点是在所有对比较期间不存储完整的图形。

表5。

内存使用量(GB)

软件SRR950457PRJEB7686型SR04R9577
星码0.6530.95.2
种子53.9
幻灯片1.3013.9
彩虹0.5
Cd命中测试0.8028.5

3.3识别富集序列基序

作为一种序列聚类算法,星码还可以用于其他应用,如识别丰富的基序。序列模体被认为在DNA代谢中起着重要作用。关键调控因子,如转录因子、核小体和非编码rna具有序列偏好,将它们定位于它们的作用位点。识别这些序列是确定调节器和它们所涉及的机制的一种方法。然而,它们的序列并不完全相同,因此它们的序列并不完全相同。这个问题对于长基序(超过12-13个核苷酸)来说,由于组合标度的原因,计算上变得困难。

我们根据脑膜炎病原体做了一个测试脑膜炎奈瑟菌. 这种细菌的基因组中穿插着一个被称为DNA摄取序列的12-bp序列(史密斯et公司 艾尔。1999年). 我们从2.19 Mb基因组的两个方向提取了12个mers,得到了439万个12 mers,包含277万个独特序列。在Levenshtein距离为2的范围内用星码对12个mer进行聚类,我们鉴定出已知的DNA摄取序列N。脑膜炎(atgccgtcgaa)是最丰富的12个mer,精确命中1466次,不精确命中2096次。这一结果证明了星码可以用来识别细菌基因组中与生物学相关的基序。

为了在另一个应用程序上测试星码,我们使用了由rnacetre产生的RNA-蛋白质相互作用数据(et公司 艾尔。2009年). 哺乳动物剪接因子SRSF1已知能结合富含RNA-GA的基序,但对其识别的基序存在一些分歧(潘迪特et公司 艾尔。2013年). 对于rnacetre数据集中每个人类SRSF1的复制,我们用它们的秩替换微阵列信号,并从微阵列探针中提取10个mers。这10个mer被给予与其所属探针等级相等的分数,并使用最大Levenshtein距离为2的星码球形聚类法发现了丰富的基序。因此,最丰富的10百万富翁的分数就是这一距离内所有10百万富翁排名的总和。在6个重复中,最富集的10个mers是aggacgga、aggacgga、aggacgg、aggacgg、aggacgga和AGGATACAGG。除最后一次复制外,母题由AGGAC和GGA组成,间隔区长度可变。这表明SRSF1与RNA的结合可能涉及一个间隔序列,这解释了6-mers和7-mers的基序之间的不一致。

4讨论与结论

Starcode是一种基于所有对匹配的序列聚类算法。它可以达到很高的精度,并且在实验数据集上它可以比流行的启发式算法更快。按照设计,starcode是为在多核平台上处理高吞吐量的测序数据而定制的,这些平台有足够的内存。由于其优越的精度和更快的运行时间,它填补了现有软件的空白,允许充分利用中高端硬件。

有点令人惊讶的是,starcode在实验数据集上比竞争工具快得多,而seed和cd-hit在人工数据集上的速度更快。Starcode是从TRIP实验数据集基础上开发的,并选择poucet搜索来给出最佳的经验结果。我们通常从实验数据中推测出随机结构的好处。

星码的速度和精确性也使得它对于其他的聚类任务也很有用,比如在微生物基因组和实验数据中识别丰富的基序。这里,我们给出了两个此类应用的例子。首先,我们在N。脑膜炎. 在第二部分中,我们恢复了人类RNA结合蛋白SRSF1的基序,并注意到它似乎由一个连接子分开的两部分组成。这一假设与SRSF1通过两个连续的RNA识别基序结合RNA的事实是一致的,这两个序列已知能结合3-4个核苷酸(多布内et公司 艾尔。2013年). 包含插入和删除的Levenshtein距离比位置权重矩阵表示更有可能捕获双部分结合基序。使用聚类方法来解决这个问题是不寻常的,但它说明了基于距离的方法的潜在优势。

星码看起来比其他工具快的原因之一是它被设计成聚集相对相似的序列。当对相关序列进行聚类时,Levenshtein距离将不得不增加,从而导致运行时间成倍增长(图5b) 是的。然而,对于纠正测序引入的错误这一重要的实际案例,starcode说明了仍有空间开发比当前技术更快速和更精确的算法。

致谢

我们要感谢Maria Chatzou对初稿的宝贵反馈,感谢Heng Chang Chen的表演果蝇旅行实验。

基金

这项研究得到了加泰罗尼亚政府(经济和知识部)和西班牙经济和竞争力计划(Nacional BFU2012-37168)(Centro de Excelencia Severo Ochoa 2013-2017)的支持(SEV-2012-0208). P、 C.奖学金部分由西班牙经济和竞争力部提供支持【国家培训次级方案:博士生培训博士学位前研究金,2013年】。

利益冲突声明:无。

工具书类

  • Akhtar W.等人。(2013年)。数千名记者并行整合检测染色质位置效应.细胞 , 154公元914-927年。[公共医疗][谷歌学者]
  • 鲍娥等。(2011年)。种子:下一代序列的有效聚类.生物信息学 , 27公元2502-2509年。[PMC免费文章][公共医疗][谷歌学者]
  • Chong Z.等人。(2012年)。Rainbow:一个集成的工具,用于高效地聚类和组装RAD-seq-reads.生物信息学 , 28第2732-2737页。[公共医疗][谷歌学者]
  • Daubner G.M.等人。(2013年)。RRM-RNA识别:核磁共振或晶体学和新发现.货币。奥平。结构。生物。,23100-108年。[公共医疗][谷歌学者]
  • Dohm J.C.等人。(2008年)。高通量DNA测序中超短读取数据集的重大偏差.核酸研究。,36,e105。[PMC免费文章][公共医疗][谷歌学者]
  • Eid J.等人。(2009年)。单聚合酶分子的实时DNA测序.科学类 , 323第133-138页。[公共医疗][谷歌学者]
  • Fu L.等人。(2012年)。CD-HIT:加速下一代测序数据的聚类.生物信息学 , 28,3150–3152年。[PMC免费文章][公共医疗][谷歌学者]
  • MacKay D.J.C.(2002年)。信息论、推理与学习算法. 剑桥大学出版社,纽约。[谷歌学者]
  • Margulies M.等人。(2005年)。微加工高密度picolitre反应器的基因组测序.自然 , 437公元376-380年。[PMC免费文章][公共医疗][谷歌学者]
  • Nakamura K.等人。(2011年)。Illumina定序器的序列特定误差分布.核酸研究。,39,e90。[PMC免费文章][公共医疗][谷歌学者]
  • Neederman S.B.,Wunsch C.D.(1970年)。一种用于寻找两种蛋白质氨基酸序列相似性的通用方法.J、 分子生物学。,48,443-453年。[公共医疗][谷歌学者]
  • Pandit S.等人。(2013年)。全基因组分析揭示了SR蛋白在调控剪接中的合作与竞争.分子电池 , 50公元223-235年。[PMC免费文章][公共医疗][谷歌学者]
  • Ray D.等人。(2009年)。RNA结合蛋白RNA识别特异性的快速系统分析.纳特。生物技术。,27公元667-670年。[公共医疗][谷歌学者]
  • Schepers K.等人。(2008年)。用细胞条形码技术剖析T细胞谱系关系.实验医学学报,205公元2309-2318年。[PMC免费文章][公共医疗][谷歌学者]
  • Shimizu K.,Tsuda K.(2011年)。SlideSort:所有配对的相似性搜索,用于短阅读.生物信息学 , 27公元464-470年。[PMC免费文章][公共医疗][谷歌学者]
  • 史密斯H.O.等人。(1999年)。自然转化细菌的DNA摄取信号序列.微生物学研究。,150公元603-616年。[公共医疗][谷歌学者]
  • Ukkonen E.(1995年)。后缀树的在线构造.算法 , 14公元249-260年。[谷歌学者]
  • 吴S,Manber U.(1992年)。快速文本搜索:允许错误.公社。ACM公司 , 3583-91年。[谷歌学者]

文章来自生物信息学在这里提供牛津大学出版社