LinearFold:RNA二级结构的线性时间预测

来源:谷歌翻译: https://www.biorxiv.org/content/10.1101/263509v2.full

摘要

快速而准确地预测RNA序列的二级结构在许多应用(例如药物设计)中很有用。最新的预测器有一个基本的局限性:它们的运行时间随输入序列的长度呈三次方缩放,这对于较长的RNA来说很慢,并且限制了二级结构预测在全基因组应用中的使用。为了解决这个瓶颈,我们设计了第一个线性时间算法来解决这个问题。可以与热力学评分和机器学习评分功能一起使用。像以前的工作一样,我们的算法基于动态编程(DP),但有两个关键区别:(a)我们以从左到右而不是自下而上的方式递增地处理序列,并且(b)由于这种增量处理,我们可以进一步采用波束搜索修剪来确保实践中的线性运行时间(以精确搜索为代价)。尽管我们的搜索是近似的,但令人惊讶的是,它在具有已知结构的各种序列数据库上甚至可以提供更高的总体准确性。更有趣的是,它导致对该数据库中最长序列家族(16S和23S核糖体RNA)的预测更加准确,并且提高了远程碱基对(相距500+个核苷酸)的准确性。

作者贡献:LH根据DH的建议构想了这个想法。LH,DD和KZ设计了算法。LH和DD在Python中实现了一个原型。DD和KZ在C ++中实现了快速版本。DHM和DH对算法进行了监督测试。DD进行了测试并绘制了数字。DD,LH和DHM撰写了手稿;卫生署对此进行了修改。

核糖核酸(RNA)参与许多细胞过程。尽管许多RNA编码蛋白质(信使RNA,mRNA),但非编码RNA(ncRNA)具有固有功能,不会翻译成蛋白质(1)。ncRNA的序列催化的反应(23),调节基因表达(4 –  6),提供一种用于蛋白质(位点识别78)和蛋白质(的投放服务9)。最近发现并鉴定了各种类型的长非编码RNA,即长于200 nt的 ncRNA (10),在确定其功能和作用机制方面提出了新的机遇和挑战。此外,RNA既是遗传物质又是功能分子的双重性质导致了RNA世界假说,即RNA是生命的第一个分子(11),并且这种双重性质也已被用于开发体外方法来进化功能。序列(12)。最后,RNA是一种重要的药物靶标和剂(13 – ⇓ ⇓ ⇓  18)。

预测RNA序列的二级结构,定义为所有规范碱基对的集合(A–U。G–C,G–U)。是一个重要的和具有挑战性的问题(1920)。知道结构揭示了关于RNA的功能,这是在许多应用中,从非编码RNA的检测(有用的关键信息21 –  23),以寡核苷酸的设计用于消息的敲低(2425)。鉴于基因组数据的大量增加(每年约10 21个碱基对),能够快速确定结构很有用(26),并给出了具有实验确定的结构的序列的一小部分。实验测定可以提供一种能够提高预测RNA二级结构(的准确性信息27),以及这些分析现在可用于全转录物和在体内28 –  30)。最近的研究侧重于预测的准确度提高(31 。⇓ 。⇓ 。⇓ – 35),但没有对预测的速度足够重视。

图。1。

RNA二级结构和我们工作的高层次构想。左上方:大肠杆菌 tRNA Gly的二级结构;右上:对应的圆图;中央:对应的点括号格式。下图:我们工作的示意图。简而言之,我们的算法从左到右扫描序列,并将每个核苷酸标记为“。”(未配对),“(”(与未来的核苷酸配对)或“)”(与先前的核苷酸配对) 。

重要性声明

快速而准确地预测RNA二级结构(规范碱基对的集合)是一个重要的问题,因为RNA结构揭示了有关其功能的关键信息。现有方法对于相对较短的RNA可以达到合理的准确性,但是它们的运行时间几乎与序列长度成三次方,对于较长的RNA来说太慢了。我们开发了用于RNA二级结构预测的第一个线性时间算法。出人意料的是,我们的算法不仅运行速度快得多,而且还提高了结构已知的各种RNA序列的整体准确性,这些改进对长RNA家族(如16S和23S核糖体RNA)具有重要意义。更有趣的是,它对于远程碱基对也更加准确。

虽然有两个主要的方法模拟 RNA二级结构,即经典热力学方法(3637)以及最近基于机器学习的方法(3839),所有这些努力使用几乎相同的动态规划(DP)算法(4041),以找到最好的计分的结构。然而,该算法,从计算语言学(借用4243),具有一个运行时间øÑ 3),该秤立方体与序列长度ñ。对于长RNA来说这很慢(n > 1,000),并且在实践中,许多研究人员求助于在整个序列中的短区域上运行该算法,这不可避免地忽略了跨片段的碱基对(44)。计算和实验研究表明,天然RNA序列末端之间的碱基配对是可预期的。

在本文中,我们设计了第一个线性时间RNA二级结构预测算法LinearFold,其灵感来自于我们先前在线性时间自然语言解析中所做的工作(45)。虽然经典的O3)动态编程是自下而上的,为每个跨度解决了最佳子结构,但我们的算法却是从左到右,以点括号格式(未配对的“。”,以“ (”或“。”结尾)。虽然此天真的版本在指数时间O(3 n)中运行,但我们使用了一种从计算语言学(46)借用的有效合并方法,该方法将运行时间减少到On3)。在此从左到右算法的基础上,我们进一步应用波束搜索,这是一种流行的启发式方法,用于修剪搜索空间(45),该算法仅在每个核苷酸处保留得分最高的前b个状态,从而得出O(n)时间近似搜索算法。即使我们的搜索不是精确的,根据经验,在合理的波束大小(例如b = 100)下,它也接近精确搜索,实际上比精确搜索具有更好的预测准确性。

我们的算法可以与热力学模型和机器学习模型一起使用。特别是,我们使用了Vienna RNAfold(37)的热力学自由能模型实现了LinearFold的两个版本,LinearFold-V和CONTRAfold(38)的机器学习模型实现了LinearFold-C(见图1(下))。我们在结构完善的RNA序列的不同数据集上评估我们的系统,结果表明,尽管LinearFold效率更高,其所有族群的平均准确度更高,而LinearFold令人惊讶相对于最长的16S和23S核糖体RNA家族的精确搜索方法,其精确性要高得多。更有趣的是,LinearFold在相距超过500个核苷酸的远程碱基对上也更准确,这对当前模型而言是众所周知的挑战性问题(47)。

结果

LinearFold的效率和可伸缩性

为了证明LinearFold的效率和可伸缩性,我们将其运行时间与基线系统中使用的常规立方时间预测算法CONTRAfold和Vienna RNAfold进行了比较。图2显示了两个数据集的结果:(a)ArchiveII数据集(48),具有已知结构的各种RNA序列集(请参见“方法”部分和表SI 1的详细信息),以及(b)(采样的)子集RNAcentral(49)的文献,这是来自许多数据库的全面(元)ncRNA序列集。尽管ArchiveII集包含长度为3,000或更短的序列,但RNAcentral集具有很多更长的序列,最长的是244,296 nt(来自NONCODE数据库的Homo Sapiens成绩单NONHSAT168677.1(50))。我们使用的机器运行的Linux是3.40GHz Intel Xeon CPU和32G内存。所有程序均使用GCC 4.9.0编译的C / C ++编写。

表SI 1。

ArchiveII数据集的详细信息以及CONTRAfold MFE,LinearFold-C,Vienna RNAfold和LinearFold-V的预测准确性。统计显着性由*(0.01≤标记p <0.05)和**(p <0.01)。

图2。

答:在ArchiveII数据集上的运行时间比较:LinearFold-C(光束大小为100和50)与两个基线CONTRAfold MFE和Vienna RNAfold(LinearFold-VLinearFold-C具有相同的运行时间)比较。B:RNAcentral数据集上的运行时间比较(对数对数)。C:内存使用情况比较(RNAcentral set,log-log)。

图2 A确认LinearFold的运行时间与序列长度成线性比例,而两个基线系统则超线性地变化,其经验运行时间为O2.6),由曲线拟合确定。图2 B在更长的序列(对数尺度)上证实了这一事实,对于约10,000 nt的序列(例如,HIV基因组),LinearFold(默认光束大小为b = 100)仅花费7秒,而基线需要4分钟。对于长度为32,753的序列,我们的LinearFold仅需23秒,而CONTRAfold需2个小时,RNAfold需1.7个小时。这清楚地显示了我们的线性时间预测算法在很长的ncRNA上的优势。

此外,LinearFold在内存使用方面也具有优势,从而可以在极长的序列上实现更好的可伸缩性。基线三次时间算法需要随序列长度按平方缩放的存储空间,因为从直观上讲,它们需要为整个序列的每个子串[i,j]找出最佳得分或最小自由能结构。这意味着如果序列长度加倍,则需要4倍的内存。另外,由于设计上的缺陷,CONTRAfold MFE或Vienna RNAfold均不能以大于32,767 nt的任何序列运行。另一方面,LinearFold不仅占用线性时间,而且还使用线性内存,而无需使用大小为O的二维表(2)。结果,LinearFold能够在不到3分钟的时间内处理RNAcentral中最长的序列(244,296 nt)。实际上,在我们的32GB内存机器上,LinearFold甚至可以扩展到10,000,000 nt的序列。

LinearFold的精度

接下来,我们比较LinearFold和两个基线系统的预测精度,报告ArchiveII数据集中每个RNA家族的正预测值(PPV;已知结构中预测对的分数)和敏感性(预测已知对预测的分数)。在先前的工作之后,我们还使用配对的双尾t检验检验了统计显着性(51)。图3表明,与所有系列的平均值相比,LinearFold-C的CONPfold上的 PPV和灵敏度提高了+2.1%/ + 1.4%(绝对值)。这是令人惊讶的,因为LinearFold使用少量的运行时即可生成更准确的结构。分别是LinearFold-C在三个家族中,PPV /敏感性的PPV /敏感性显着提高:I组内含子,16S和23S核糖体RNA,最后两个是该数据集中最长的家族。16S rRNA的平均长度为1548 nt,PPV /灵敏度绝对提高+3.89%/ + 3.08%,而23S rRNA的平均长度为2927 nt,绝对值为+14.00%/ + 9.98%。更令人惊讶的是,LinearFold的精度提高在更长的序列上更为明显。LinearFold-V也提高了精度与维也纳RNAfold相比,差异较小(PPV /灵敏度的绝对值整体提高了+0.2%/ + 0.2%)。单独地,对I族内含子和16S rRNA这两个家族的PPV /敏感性的改善是显着的。

图3。

在ArchiveII数据集上的PPV和灵敏度(按族),将LinearFold与相应的基线CONTRAfold MFE和Vienna RNAfold进行比较。每列代表一个族的准确度,在该族的所有序列中取平均值。表中报告了所有家庭的平均准确度。统计显着性标记为*(0.01 < p <0.05)和** {p <0.01)。有关准确度的详细信息,请参见表SI 1。有关PPV /灵敏度指标和重要性测试方法的详细信息,请参见“方法”部分。

光束尺寸的影响

上面我们使用b = 100作为默认光束大小。现在我们研究不同光束大小的影响。我们首先研究搜索质量的影响。由于我们的LinearFold算法使用近似搜索而不是精确搜索,因此我们将精确搜索自由能和返回自由能之间的差用作搜索质量的度量-它们越接近,我们的搜索质量就越好。从图4 A中我们可以看到,当光束尺寸增加时,搜索质量越来越近,而Linear-Fold使用默认的光束尺寸实现了类似的模型成本/自由能。同样,图4 B在每个光束尺寸下绘制预测的对数。它表明ViennaRNA模型倾向于高估,而CONTRAfold模型则低估。同样,尽管与两种模型中的精确预测相比,我们的算法总是预测不足,但是当使用默认波束大小时,我们预测的碱基对数与每种模型几乎相同。当b = 100 时,平均差为0.29(CONTRAfold)/0.01(ViennaRNA)对。

图4。

光束大小的影响。AD说明了光束尺寸增加时不同变量的趋势。答:内部成本,即两个版本的LinearFold的平均自由能/模型成本;B:将这些方法预测的对数(按序列平均),与地面真相进行比较;C:PPV和灵敏度随光束尺寸的增加而变化;D:改变光束大小时,PPV-灵敏度的权衡;EF:ArchiveII数据集中的PPV和对对距离的灵敏度,将LinearFold-C与CONTRAfold MFE 进行了比较。每个点代表特定长度范围内所有碱基对的总体PPV /灵敏度。

图4 C绘制了PPV和灵敏度与光束大小b的关系图。LinearFold-C优于CONTRAfold MFE在两个PPV和灵敏度与b ≥75(尽管它会收敛到CONTRAfold MFE当b ⟶+∞)。和LinearFold-C的 PPV /灵敏度与稳定b ∈[100150]。图4 D显示了PPV和LinearFold灵敏度的权衡,以及光束大小的变化。它从PPV和灵敏度的增加开始,在b = 120 达到峰值,然后下降直至收敛到精确搜索。虽然峰值发生在光束大小为120时,但是我们可以看到性能LinearFold 当其光束大小为[100,150]时,它是一致的,因为两个PPV /灵敏度保持几乎相同。

远程碱基对的准确性提高

我们进一步评估了碱基配对核苷酸之间的距离对预测性能的影响。如图图4 E和F预测长途对时,LinearFold可以超越在两个PPV和灵敏度先前方法。与LinearFold对于长距离对无法准确预测的担心相反,即使在超过500 nt的配对距离下,LinearFold仍然优于以前的方法。支持信息中对LinearFold-V和Vienna RNAfold在PPV,灵敏度和预测质量方面的详细比较。

预测示例:I组内含子,16S和23S rRNA

我们比较了LinearFold-C和CONTRAfold MFE(图5将来自不同RNA家族的3个实例的预测二级结构可视化,它们分别是Group I内含子C. Mirabilis,16S rRNA A. Pyrophilus, 23S rRNA E. coli。通过正确预测更多碱基对并减少错误预测,圆形图显示了我们的性能提升。该可视化还展示了LinearFold对长距离碱基对的预测优于基线的预测,如I组Mirabilis(底部一半,对距离250),16S A. Pyrophilus(左侧,550),23S E. coli(左侧) ,600)。图SI 3显示了将LinearFold-V与Vienna RNAfold 进行比较的相应结果。我们还构建了一个Web演示,以可视化这三个系列中所有序列的结果。*

图SI 1。

该图对应于图4(c)(d),但带有ViennaRNA版本,在ArchiveII数据集上运行。左:随着光束尺寸的增加,PPV和灵敏度的趋势;右:根据光束大小,PPV和LinearFold的灵敏度。

图SI 2。

此图对应于在ArchiveII数据集上运行的带有ViennaRNA版本的Figure4 (c)(d)。与维也纳RNA折叠相比,它通过Mathews数据集中的对长度显示了LinearFold-V的PPV和灵敏度。每个条形图代表一个长度范围内所有碱基对的总体PPV /灵敏度。

图SI 3。

比较Vienna RNAfold和LinearFold-V的3个RNA序列的圆图(对应于图5)。

图5。

比较CONTRAfold MFE和LinearFold-C的3种RNA序列(从3种不同的RNA家族中选择)的圆图蓝色弧线表示正确预测的对,红色弧线表示错误预测的对,而浅灰色弧线表示基于对的未预测。每个图从顶部开始使用RNA序列的顺时针顺序。从3个不同的RNA家族中选取了这3个实例,与CONTRAfold MFE相比,LinearFold-C的性能得到了显着改善。

讨论区

RNA结构预测对于推断RNA功能非常重要,并具有许多应用,包括药物设计。现有的RNA二级结构预测算法的运行时间与序列长度成立方比例,这对于较长的非编码RNA而言太慢了。例如,这项工作中的基线系统CONTRAfold和Vienna RNAfold是最流行的两种预测软件,它们分别需要2个小时和1.7个小时的时间才能获得32,753 nt的序列。此外,现有算法还需要根据序列长度进行二次扩展的存储空间,因此,两个基线系统都无法在超过32,767 nt的序列上运行。实际上,RNAcentral数据集中最长的RNA序列为244,296 nt,是该限制的7倍。

我们使用动态编程和束搜索设计了第一个线性时间,线性内存预测算法LinearFold,并将该算法应用于机器学习模型和热力学模型。在图2中确认了时间和存储器上的线性。更有趣的是以下三个令人惊讶的发现:

  1. 1.尽管LinearFold与基线系统中的现有算法相比仅使用了少量的运行时和内存,但我们预测的结构在PPV和灵敏度以及机器学习模型和热力学模型上总体上都更加准确(请参见图3)。 。

  2. 2. 在更长的家族(例如16S和23S rRNA)上,LinearFold的准确性提高更为明显(见图35)。

  3. 3.更令人惊讶的是,LinearFold也比相距超过500 nt的远程碱基对的基线更准确(图4 E–F),众所周知,这对当前模型而言是一个具有挑战性的问题(47)。

  4. 4.尽管模型取决于波束大小b,但是碱基对的数量和预测的准确性对于波束大小的变化非常稳健(当b在100-200范围内时)。

为什么我们的波束搜索算法即使是近似的,在准确性方面(尤其是在16S和23S rRNA中)也要优于精确的搜索基线?首先,当前的热力学和机器学习模型远非完美,因此,次优结构(就自由能或模型得分而言)完全有可能比最优结构更准确(就PPV /灵敏度而言)。例如,对于约400个核苷酸的序列,可以发现约80%正确的结构,其自由能在最佳结构的5%以内(52)。但是我们的算法如何系统地选择一个更准确的次优结构而看不到基本事实?我们怀疑这是因为波束搜索会在每个步骤中修剪得分较低的(子)结构,从而要求在每个前缀处对幸存的(子)结构进行高分。此额外约束可能会补偿模型的不准确性。

我们的算法有几个潜在的扩展。首先,可以扩展LinearFold来计算分区函数和碱基对的概率,因为该任务的经典方法McCaskill算法(53)与立方时间结构预测算法相似,本文的基线。其次,线性时间算法来计算碱基对的概率应促进假结的线性时间识别,通过在那些假结预测程序(一个线性时间一个替换立方时间麦卡斯基尔算法5455)。第三,线性时间为LinearFold与使用结构化预测方法的立方时间CONTRAfold相比,该方法还便于更轻松,更快地训练参数(56),并且我们预想了适合线性时间预测的再训练模型应该更加准确。

方法

我们的LinearFold方法分为四个步骤,从最幼稚但易于理解的穷举搜索版本(图6A)开始,然后逐步扩展到线性时间版本(图6D图结构的堆栈和波束搜索。

图6。

LinearFold的四步演示,仅查找最大对数即可,而不是实际的MFE模型。每个节点代表结构的预测前缀,并显示未配对开口的堆栈。每个箭头对应一个动作(推动,跳过和弹出);死角状态带有红色×。在C语言中,具有相同颜色的节点由图结构堆栈合并。有边界的节点代表地面真实路径。

线性时间预测的基本思想是从左到右递增预测,将每个核苷酸标记为未配对的“•”,打开“(”或关闭“)”。仅考虑无假结结构。

给定一个输入RNA序列x = 1 … Ñ -i其中 ∈{A,C,G,U},我们的算法旨在找到最好的结构Y = ÿ ý 1 … ý Ñ -1,其中i {“。”,“(”,“)”}的自由能最小(或模型成本最小):嵌入式影像这里ÿ(x)是该组所有可能的结构,即,{Y | y具有圆括号},c是成本函数(即自由能函数),w是模型(和参数)。

天真的穷举增量预测: O(3 n时间。通过从左到右详尽地预测y,我们遍历y(x)中所有可能的结构,并选择具有最小自由能或模型成本的结构。我们在步骤jj∈ {0,…,n })将每个状态形式化为三元组,s =⟨σ|。i, j⟩:y,其中σ| i是迄今为止不匹配的开口组成的堆栈,其中i是堆栈的顶部,这意味着i是最后一个不匹配的开口核苷酸。ÿ是直到j – i的对应点括号(子)序列。对于每个状态,它可以通过以下三个动作之一转换到后续状态:推入,将当前核苷酸j标记为左括号“(”,将其放在栈顶,跳过,标记j作为点“。”,则使堆栈保持不变,然后弹出,如果j匹配i并从堆栈中弹出i,则将j标记为右括号“)” 。参见图4 SI(a)用于演绎系统。该算法取O(3 n),以期详尽遍历所有可能的结构(请参见图6 A)。

图SI 4。

演绎系统比较了穷举搜索算法和带有图结构堆栈的动态规划。在这里,o表示字符串连接。如果(a,b)是允许的对之一(CG,AU,GU),我们说(a,b)“匹配” 。

通过全栈合并进行动态编程: 嵌入式影像 时间。现在,我们在此详尽的方法之上应用动态编程来利用共享计算。考虑一个简单的情况下,两种状态可以合并:如果在相同的步骤中两个状态Ĵ, ⟨σ,Ĵ ⟩:Y“和⟨σ,Ĵ ⟩:Y”共享完全相同的堆σ,但不同的点状括号中的字符串yy’,我们说这两个状态是“等效的”,我们可以将它们合并(并且仅在yy’之间保持更好的得分。)图6 B说明了这种合并。尽管我们合并以减少状态数,但这仍然是指数时间,因为每个步骤中可能会有成指数的不同堆栈。该算法需要嵌入式影像 时间。

通过图结构堆栈进行动态编程: O3时间。为了避免成倍地考虑多个状态,我们进一步将状态与不同的堆栈合并。考虑在同一步骤两种状态Ĵ,⟨σ 0 | Ĵ ⟩和⟨σ 1 | Ĵ ⟩共享最后不成对开口(即,堆栈顶部)。我们称这些状态为“暂时等效”,因为直到未配对的开口i关闭(因此从堆栈中弹出),它们都可以被视为完全相同。换句话说,我们可以代表两个堆栈σ 0和σ 1 | 作为… 在哪里…表示我们目前不在乎的历史部分。堆栈的这种分解被富田(46)称为“图结构堆栈”(GSS )。合并后,我们将新状态定义为⟨i  j⟩ 因此我们保持O2)状态。对于每个状态⟨i  j⟩ 弹出动作可能花费最坏的On)时间,因为⟨i  j⟩ 可以与步骤i中的每个⟨k ,i combine合并。因此,总时间复杂度为O3)。有关合并过程的示例,请参见图6 C;对于演绎系统,请参见图SI 4(b)。

通过波束搜索进行动态编程: O(n) 时间。实际上,精确的搜索算法仍在O3)时间内运行。但是,与所有现有系统用于RNA结构预测的传统的自下而上的O3)搜索不同,这种从左到右的O3)搜索很容易“线性化” 。我们进一步采用波束搜索修剪(56)以将复杂度降低至线性时间。通常,我们只保留b个得分最高的状态⟨i ,j⟩每一步。这样,所有得分较低的状态都将被修剪,并且如果一个结构生存到最后,则它在每个步骤中都必须是前b个状态之一。修剪还意味着在弹出动作中,状态(i,j)可以与步骤i中最多b个状态(k,i)组合。因此,总时间复杂度为Onb 2)。但是,不是使用pop动作生成2个新状态,而是使用多维数据集修剪57)生成最佳b状态,这将花费Ob log b) 时间。因此,整个长度为n的序列的总运行时间为Onb log b),请参见图6 D中的波束搜索。

数据集,评估指标和重要性测试。我们选择ArchiveII数据集(48),它是3,000多种具有已知二级结构的RNA序列的集合。但是由于当前的CONTRAfold机器学习模型(v2.02)是在S处理数据集(58)上训练的,因此我们删除了出现在S处理数据集中的那些序列。我们使用的结果数据集包含9个科的2889个序列,平均长度为222.2 nt。由于比较分析中存在碱基对匹配的不确定性,因此我们认为如果碱基对也因一条链上的一个核苷酸滑移而被正确预测((48))。通常,如果一对(i,j)位于预测结构中,如果(i,j),(i – 1,j),(i + 1,j),(i,j – 1),(i,j + 1)中的一个是正确的)处于地面真相结构中。我们报告敏感性和PPV,嵌入式影像我们使用配对的双尾t检验来计算统计显着性,其I型错误率与先前的方法一致(51)。

LinearFold:RNA二级结构的线性时间预测》有1个想法

  1. 百度于2019年7月首次提出。该算法使得整序列、整基因组的 RNA 结构预测成为可能,也是 RNA 结构预测领域40年来第一次重大提速。这项工作发表于生物信息学顶级会议 ISMB 2019 和生物信息学权威杂志 Bioinformatics。

    针对此次新型冠状病毒的基因组(长达3万个碱基),采用该算法,27秒就可以预测其结构。相较于经典算法,现在只需不到半分钟就可以拿到病毒的结构资料,提升基因检测、疫苗研发等科研中心的工作效率,让病毒的研究及疫苗开发速度快速提升。

发表评论

电子邮件地址不会被公开。 必填项已用*标注