经典回顾01RNA二级结构预测的动态规划算法

Posted 2021-04-18 BioHack

经典回顾主要是记录我对经典英文文献的翻译，是我学习过程的一部分。文献翻译没有专门进行多次校正，所以还远远达不到“信达雅”的水平，但是也能勉强一看。如果有错误的地方，也欢迎读者（如果有读者的话）指出和探讨。今天翻译的文章是一篇于1981年发表在《Nucleic Acids Res》上的文章，是RNA二级结构的预测领域，非常经典的一篇文章。考虑到篇幅问题，今天的内容只翻译，不进行分析和梳理。

原文题目：Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information

原文作者：Zuker M, Stiegler P

译文作者：太昊

译文

利用热力学（thermodynamics）和辅助信息（auxiliary information）优化大型RNA序列的构象计算

摘要：本文展示了一种新的计算折叠RNA分子的方法，该方法利用已发表的堆积能（stacking energies）和失稳能（destabilizing energies）的数值寻找最小自由能构象。它基于应用数学中的动态规划算法（dynamic programming），比目前生物文献中出现的程序效率更高、速度更快，以及能够折叠更大的分子。通过对长度为459个核苷酸的免疫球蛋白γ1重链信使RNA片段的折叠，证明了它的效力（power）。我们超出了基本方法的要求，展示了如何将额外信息（additional information）合并到该算法中。额外信息包括在化学反应性（chemical reactivity）和酶敏感性（enzyme susceptibility）的数据。我们通过折叠大肠杆菌（Escherichia coli）16S核糖体RNA的两个大片段来说明这一点。

引言(Introduction)

一个RNA分子的核苷酸序列携带着它在三维空间中实际构象所需的信息。研究核苷酸链的二级结构或者折叠可能推导出该分子组织结构的初步草图。迄今为止，通过应用各种拓扑学和热力学规则去发现给定序列在能量上最支持的结构的方式，已经可以预测RNA的二级结构。Tinoco等人设计一种用于估计RNA二级结构的基本方法（1）。通过解链温度（melting temperature）以及在双链的寡聚核苷酸（oligoribonucleotides）和多聚核苷酸（polyribonucleotides）上的热力学数据，将稳定数（stability numbers）分配给预测结构。然后，几位作者通过实验完善了自由能（free energy）贡献的估算值(2-4)。稳定数让位于更精确的堆积能（stacking energies）和失稳能（destabilizing energies），从而改进了预测RNA二级结构的方法（5）。最近，Salser整理编辑了文献中获得的热力学数据（6）。

使用试图最小化自由能的计算方法预测RNA二级结构并不是什么新鲜事。Pipas和McMahon发表了利用热力学预测RNA二级结构的早期尝试（7）。这种早期技术能够轻松地处理相对较短的序列，例如转运RNA的序列，但是，由于它需要计算时间与2N成正比，因此在折叠长的RNA序列方面是低效的，其中N是序列中核苷酸的数目。Studnicka等人发表了这项技术的重大改进。他们定义了折叠RNA分子的精确拓扑规则，并且完全利用了可以获得的热动力学数据。然而，他们的算法太过复杂，并且为了获得一个解决方案，需要人为干预几个阶段。此外，与应用数学中的优化算法相比，该算法速度慢、效率低。最终，由于各种捷径和折衷方案，解决方案通常是次优选择（suboptimal），尤其是对于较长的链来说。

在数学方面，两篇相关的论文值得一提。Nussinov等人制定了一个杰出的动态规划算法，它可以在一个折叠的分子中获得最大数目的碱基配对（9）。该算法易于修改，从而用于为各种碱基配对类型分配不同的权重，但是它不能按照由Salser编制的热动力学标准来计算最优化的结构。堆积能和失稳能无法合并到此算法中。Waterman和Smith定义了一种效率较低的算法，但是考虑到了堆积能和失稳能，尽管他们的能量计算不是Studnicka等人使用的那些（8）。此外，他们的算法仅限于200个核苷酸以下的序列。

本文提出了上述两种方法的本质融合。我们使用与Studnicka等人相同的折叠规则和能量计算（8）。改进之处在于创建了一个专门为这个问题设计的动态规划算法。因此，它在微妙但重要的方面与来自Nussinov以及Waterman和Smith等人的算法不同（9-10）。这一结果产生了一个算法，该算法已经被成功地用于折叠多达600个核苷酸的序列，并且有可能处理长达800个核苷酸的序列。对于含有N个核苷酸的序列，它可以在与N3成正比的时间内计算出最优结构，与Studnicka等人的程序的N5相比（8），这是一个真正的改进。对于长达200个核苷酸的序列，计算时间几乎可以忽略不计。该程序不需要人工干预，且只需一个计算周期。最终，由于不采取任何节省时间或其他方式的折衷办法，因此根据所使用的规则，所生产的结构具有尽可能少的自由能。

我们在算法中加入了额外的信息，从而超越了通常的折叠规则。具体来说，当有关于某些核苷酸对化学修饰的反应性的信息时，或者当酶学研究表明在部分水解的条件下，哪个磷酸二酯键最有可能被切断时，我们可以将这些信息直接构建到算法中，并预测与这些数据相容的最优结构。即使是二级结构保守性的系统发育数据，或是RNA酶切检测（examination of RNA digests）获得的特定长程相互作用的证据，也可以被该算法使用。事实上，我们基本的假设是使用这种额外信息不仅是可取的，而且是必要的。无论使用多么强大的计算机程序，目前可用的基本折叠规则和热力学信息都不足以非常可靠地预测正确的二级结构。使用的额外信息越多，预测的折叠效果越好。这将在稍后关于E.coli 16S核糖体RNA折叠的讨论中加以说明。

方法（Methods）

A. 定义(Definitions)

在本节中，我们定义将从中选择最佳结构的一组结构。尽管我们的折叠规则与Studnicka等人使用的折叠规则相同（8），但是我们的方法是制定非常简单的规则，这些规则甚至允许一些不可能的结构。通过为它们分配高的能量，可以消除这种无效结构。结构的图形化表征（schematic representation），即在数学理论中的图（graphs），被引入。这些表征实际上与Nussinov等人采用的那些一致（9）。虽然对理论而言不是必需的，但它们有助于准确地阐明正在考虑的结构，并且是一种展示如何计算结构能量的好方法。

按照惯例，我们从分子的5'末端开始给RNA分子的核苷酸编号，用Si表示第i个位置的核苷酸（1<=i<=N），其中N是核苷酸总数。字母S单独表示整个分子，而Sij表示从Si到Sj包括的核苷酸，假设1<=i<j<=N。图1显示了在半圆上等距排列的RNA分子的未指定核苷酸。在这个几何图中，分子的N个核苷酸被称为顶点（vertices）。碱基间半圆的N-1个圆弧(arcs)被称为外部边（exterior edges），并且它们代表在连续核苷酸之间的磷酸二酯键。碱基配对由核苷酸之间的线段表示，即由两个顶点之间的半圆上的弦(chords)表示。如果一个弦连接两个可以碱基配对的核苷酸：G-C，A-U或G-U，则该弦被称作可接受的（admissible）。这些弦称为内部边（interior edges）。边和顶点的全部集合称为图（graph）。一个可接受的结构（admissible structure）被定义为一个具有仅含有互不接触或相交的可接受弦组成的图的结构。不允许弦接触等价于一个核苷酸可以与最多一个其他核苷酸发生碱基配对。禁止两个弦相交是更为深刻的。这个条件正是排除所有打结结构（knotted structures）所需的条件。使用Studnicka等人的命名法（8），我们的可接受结构都是“正交的”（orthodox），并且包括所有他们的“复杂高级结构”（complex hyperstructure）。

通过定义结构的自由能或等效地定义其图来完成描述。一个图的面（face）被定义为在所有侧面（sides）上都被边（edges）限制的的任何平面区域。结构的自由能不与健相关，但是与键之间的区域有关，如Nussinov等人所做的那样（9）。就图形表征而言，能量取决于图的面。具有单独内部边的面叫做发夹环（hairpin loop）。有两个内部边的面被分为三个组：如果内部边在两个侧面被单个外部边分开，则该面称为堆积区（stacking region）；如果它们在一侧被一个外部边隔开，但在另一侧被一个以上外部边隔开，则该面称作突起环（bulge loop）；否则，该面被称作内部环（interior loop）。最后，具有三个或更多内部边的面称为分叉环（bifurcation loop）。“发夹”（hairpin）一词是指一类结构或亚结构，该结构的面包含一系列连续的堆积区，突起环或内部环，并以发夹环结尾。图1给出一个结构及其对应的图的示例。所有不同类型的区域都有出现。从核苷酸6到17和核苷酸18到32形成的亚结构都是发夹结构。

如果F是一个面，我们让E(F)表示它对应的自由能。例如，如果F是在两个A-U碱基对之间的堆积区，则E(F) = -1.8 kcal/mole（6）。一个结构的能量等于它的面的能量的总和，当前的问题是找到具有最小自由能的结构。能量函数是一个强大的工具。正确使用它可以排除无效或不想要的结构。例如，如果F是少于四个外部边的发夹环，则设E(F) = ∞ 。这将排除选择在空间上（stearically）不可能的结构，包括少于三个碱基的发夹环。甚至是仅在一系列堆积区的内部允许G-U碱基配对的更复杂的规则，也可以通过正确使用能量函数来实现。另一方面，能量函数可以用来支持某些类型的结构。这将会在稍后讨论。我们使用Salser编制的能量数据（6）。迄今为止，这是有关堆积区的稳定效应（stabilizing effect）和各种环的失稳效应（destabilizing effects）的最佳信息。关于分叉环的失稳效应，实际上仍然是未知的。我们随后以两种不同的方式处理它们。在这个算法的一个版本中，它们都具有零能量，而另一个版本，他们被视为内部环。后一种算法当然更复杂。本文显示的所有结果都是使用第二种方法从计算机程序生成的。

B. 算法(The algorithm)

现在，我们将详细描述当分叉环被赋予零能量时如何获得二级结构的最小自由能。该算法简单但功能强大。没有折衷方案。不会忽略任何一种可能性，但是该算法从许多结构中选择了一种最小能量的结构，即使对于5S核糖体RNA如此小的分子，也是如此。主要的数学技巧是为给定的RNA序列的每个子序列S ij计算两种可能不同的能量。对于所有满足1<=i<j<=N的成对的i和j，令W(i,j)为子序列S ij形成的所有可能的可接受结构的最小自由能。此外，令V(i,j)为S i和S j互相配对时，由S ij形成的所有可能的可接受结构的最小自由能。如果S i和S j不能配对，那么V(i,j) = ∞。V(i,j)和W(i,j)的数字是递归计算出来的，首先是对所有五核苷酸子序列，然后用于所有连续的越来越大的S的子序列。五核苷酸序列是非常容易处理的。它们没有形成稳定的结构，所以如果j-i=4，则W(i,j)=0。如果在一个3个核苷酸发夹环的情况下，S i和S j是G-C或A-U碱基对，则V(i,j) = +8.4 或 +8.0 kcal/mole（6）。如果j-i = d > 4，则在不同的满足j'-i'< d的i'和j'配对情况下，依据V(i',j')和W(i',j')可以计算V(i,j)和W(i,j)。这些数据已经被计算出来了。设想在S ij上能量V(i,j)的可接受结构，假设S i和S j可以配对。我们用FH(i,j)表示包含在S i和S j之间的内部边的发夹环，FL(i,j,i',j')表示恰好包含两个内部边的面，其中一个在S i和S j之间而另一个在S i'和S j'之间（假设i<i'<j'<j）。FL表示的面是堆积区或者是其他突出环或内部环。如图2A所示，与S i和S j间的边相邻的面是三种可能的类型之一。它有一个，两个或两个以上的内部边。在第一种情况下，V(i,j)=E(FH(i,j))。在第二种情况下，对于满足i<i'<j'<j的一些成对的i'和j'，V(i,j)=E(FL(i,j,i',j')) + V(i',j')。在最后一种情况下，对于一些满足i+1<i'<j-2的i'，V(i,j)=W(i+1,i')+ W(i'+1,j-1)。在这里，能量分为了两个亚结构能量之和；故此，用“分叉”（bifurcation）表示。因此，V(i,j)可以通过这三种方式获得的能量中的最小值，所以我们可以写成V(i,j)=min{E1,E2,E3}，这里

E1 = E(FH(i, j))

E2 = min{E(FL(i,j,i',j')) + V(i',j')},i<i'<j'<j

E3 = min{W(i+1,i')+ W(i'+1,j-1)},i+1<i'<j-2

现在设想在S ij上能量W(i,j)的可接受结构。同样，三种可能性。如图2B所示，S i或S j（或都是）不参与到这个结构中，或他们彼此碱基配对，或他们都配对，但不是和彼此。第一种情况是不重要的。该结构有至少一个悬挂端（dangling end），并且W(i,j)=W(i+1,j)或者W(i,j-1)。在第二种情况下，W(i,j)=V(i,j)，已经被计算。最后一种情况称为开放的分叉结构因为这个结构分成了两个独立的部分，在S i和S j之间的没有连接。如果S i和S i'碱基配对，并且S j和S j'碱基配对，其中i<i'<j'<j，那么W(i,j)=W(i,i')+W(i'+1,j)=W(i,j'-1)+W(j',j)。因此，从这三种情况可获得的最小能量W(i,j)由下式给出：

W(i,j)=min{W(i+1,j),W(i,j-1),V(i,j),E4}，其中E4=min{W(i,i')+W(i'+1,j)}（i<i'<j-1）

启发式地（heuristically），此递归算法的工作方式是一次向序列添加一个核苷酸，然后观察每个步骤的最优结构。要计算的最后一个数字W(1,n)是所需的答案。它是一个在S上的可接受结构的最小能量。然而，实际上，在计算W(1,n)上所花费的劳动产生了更多的能量，因为在S的每个子序列上的可接受结构的最小能量也是已知的。剩下的就是结构的构造，这相当于识别对应的图的内部边。这是通过对矩阵W和V的回溯（traceback）实现的，很简单。

C. 计算技巧(Computation techniques)

这个算法已用Fortran语言编写好了。它的实现是在带有IBM TSS/370分时操作系统的IBM 3032处理器上。该程序的最有效版本将V和W矩阵的能量存储为同一正方形数组中的一半整数，以节省空间。能量计算是在一个小的子例程中完成的，该子例程有助于更改规则和创建程序的特殊版本。

能量函数起着比单纯定义热力学研究中的能量更重要的作用。我们已经表明它被用于执行某些拓扑和折叠规则，例如发夹环的最小尺寸和G-U碱基配对的允许位置。通过给相邻的面提供非常大的正能量，可以排除禁止的碱基配对。相似地，可以通过将非常大的负能量分配给相邻面来强制进行一些碱基配对。在这两个极端之间，可以通过使用较小的“加成”（bonus）能量来支持某些碱基配对。大肠杆菌16S核糖体RNA折叠的讨论中给出了一个例子，在那里我们展示了如何将某些系统发育数据纳入算法。能量函数也能够被用来将有关化学修饰和酶可及性(enzyme accessibility)的其他信息纳入算法中。在大肠杆菌16S核糖体RNA的例子中，有关化学反应性信息来自于单链鸟苷残基的乙氧二羟丁酮修饰(kethoxal modification of single stranded guanylic residues)（12,13）。酶的可及性数据来自蛋白质-RNA复合物（14,16）或30S核糖体亚基中16S核糖体RNA的拓扑学研究（P. Stiegler, P. Carbon and C. Ehresmann, 未发表结果）。任何可以被乙氧二羟丁酮修饰的鸟苷残基完全不许碱基配对。酶学数据的使用并非那么简单。T1核酸酶切割鸟苷残基的3'末端的磷酸二酯键。胰核糖核酸酶相似地作用在胞苷(cytidylic)和尿苷(uridylic)残基上。无论哪种情况，在部分水解的条件下，仅一部分化学键被切断，表示单链区域。提供给计算机程序的指令是识别的核苷酸和下一个核苷酸（向3'末端方向）不能同时碱基配对。实际上，对包含这两个核苷酸的堆积区域赋予非常大的正能量。因此，我们在碱基配对区域的末端允许一个可及的键。信息被直接嵌入到能量函数中，不浪费时间计算不需要的结构的能量。

该程序的几个特殊版本已被编写。一个不允许分叉环，并计算一个最佳结构，该结构是一连串的发夹结构。这种结构称为开放式结构（open），它们在识别局部结构中很有用。这个程序非常快。另一种版本允许从序列中切除指定的子序列，并折叠剩余的核苷酸。当我们已经对某些子序列的结构充满信心并正在寻找长距离的相互作用时，可以使用该方法。

这里有两类输出结果。一种是结构的实际计算机绘图。稍后在图5中给出一个示例。这个程序也产生一个详细的输出结果，列出了所有连续堆积区和所有不稳定环的集合，以及相应的负能量和正能量。

结果与讨论（Results and Discussion）

A. 折叠免疫球蛋白γ1重链信使RNA片段的459个核苷酸（Folding 459 nucleotides from an immunoglobulin γ1 heavy chain messenger RNA fragment）

为了说明我们的算法及其相关计算机程序的强大功能，我们使用Studnicka等人（8）的程序的稍加修改的版本，折叠了与Rogers等人（11）所用的相同的免疫球蛋白重链γ1信使RNA片段。这459个核苷酸片段的折叠如图三中所示。计算时考虑了分叉环的失稳效应（参见方法）。我们的计算结构的自由能为-181.4 kcal/mole，比Rogers等人使用相同能量计算发现的最小自由能提高了15％（-158.5 kcal/mole）。此外，我们的程序运行快且容易使用。仅需134秒即可计算出这459个残基片段的最稳定二级结构。只有一个计算周期，仅需要核苷酸序列作为输入。因此，Studnicka等人所需要的许多术语，例如“组数”、“分支迁移”和“主区域”，不会出现在我们更简单但功能更强大的算法中。没有为了节省时间而做出折衷方案；它们是不需要的。故此，将最小自由能提高15％也就不足为奇了，因为我们正是按照上述方法解决了这个问题。最值得注意的是，我们的计算结构和Rogers等人提出的最佳折叠几乎完全没有相似性（11）。在图3中，在两个提出的折叠中出现的碱基配对区域用较粗的线表示。在我们的结构中，只有29%的碱基配对可以在Rogers等人的相应的碱基配对中找到。我们并不声称我们生成的结构比他们的“好”。给出该示例主要是为了说明我们方法的强大功能。它提出的更深层次的问题是在已经生成的折叠中，如何选择一个而不是另一个。

B. 次优级和额外信息（Suboptimality and additional information）

“最优结构”（optimal structure）这个词指的是最小自由能的折叠。这种结构不一定是独一无二的。如果分子的大小等于或大于5S核糖体RNA，并且使用了相当复杂的能量函数，计算出的能量只有kcal/mole的十分之一，因此不可能有两个或多个具有完全相同能量的最佳结构，但是它可以发生。可以设计一个回溯（traceback）程序来寻找最小能量的所有解决方案，但我们还没有这样做。那是徒劳的。对于一个200个核苷酸的序列，即使消除了琐碎的亚结构，在最小自由能的5%或10%范围内可能有许多结构。从对寡核苷酸（oligoribonucleotides）和多聚核苷酸（polyribonucleotides）的研究中推断出的能量计算，不能被认为精确到足以区分接近最小能量的各种结构。展示使用特定算法任意生成的几个次优结构的列表是没有意义的。为什么显示这些结构而不是其他结构？在计算有意义的替代结构之前，必须提供额外信息。在我们对大肠杆菌16S核糖体RNA的研究中，这些额外信息来自于与特定化学试剂具有反应性或易受核糖核酸酶攻击的核苷酸数据（参见方法）。因此，这些核苷酸不参与碱基配对或位于碱基配对区的末端。我们的策略是通过适当使用能量函数将这些信息直接合并到算法中（参见方法）。通过考虑到辅助信息从而折叠分子产生了改进的结构。然后仔细检查预测的结构，需要通过生物学的评估来选择将包含在最终模型中的特征。

在上述免疫球蛋白重链信使RNA片段的情况下，除了核苷酸序列之外，我们没有其他信息，并且将我们的结构作为对Rogers等人提出的折叠的一种暂定的改进。然而，很明显，这个序列有更多不同的结构，能量接近-181.4kcal/mole，这可能在生物学上更有意义。

C. 应用于大肠杆菌16S核糖体RNA（Application to E.aoli 16S ribosomal RNA）

本文所描述的计算机程序的主要应用是促进对大肠杆菌16S核糖体RNA二级结构的研究。关于核糖体小亚基内RNA链折叠及其相关结构组织的实验数据已经积累（12-18和P.Stiegler，P.Carbon，C.Ehresmann，未发表的结果）。这些拓扑结果提供了正确建模所需的生物学的信息。

一次折叠16S核糖体RNA中的所有1542个核苷酸超出了我们的计算机程序的能力，该程序可以折叠最多800个核苷酸的分子。因此，我们选择了对16S核糖体RNA序列进行细分，不是任意细分，而是根据描绘分子中不同结构域的拓扑学研究分类。例如，我们在这里给出了前574个核苷酸的计算机折叠结果。该部分形成了核糖体蛋白S4的RNA相互作用位点，该部分能够单独结合到16S核糖体RNA上（14，15）。以酶为探针的拓扑学研究也表明，即使没有蛋白质，这个结构域也存在于RNA分子中（14，16）。

使用我们最复杂的程序折叠了574个核苷酸的RNA片段，该程序考虑了该序列的所有化学和酶信息（12、13、16）。结果如图4所示。箭头指向通过化学修饰或酶数据指出的可及性。这是可以生成的与辅助信息一致的最稳定结构。然而，计算结构缺乏“生物现实性”（biological reality）。它不能完全满足对这个RNA结构域的详细结构研究得到的实验观察，而该实验结果无法充分量化以直接用于该程序。即使如此，也有许多结构的模序（motif）可以从这个相当复杂的结构中挑选出来。实际上，图4中所示的所有碱基配对中约有一半被保存在16S核糖体RNA二级结构的暂定模型中（P.Stiegler，P.Carbon，M.Zuker，J.P.Ebel，C.Ehresmann，正在准备的手稿）。如方法部分所述，使用从基本程序派生的专用程序执行优化。一项任务是搜索局部结构。如果有人认为在体内核糖核酸分子是从5'末端到3'末端依次合成的，那么这样的研究是合理的。因此，正如最近对转运RNA（20）的研究所建议的那样，顺序合成可以确保分子的正确折叠。在前574个核苷酸或子序列上寻找局部结构是通过寻找一个最优的开放（open）结构来完成的。每个发夹的大小被任意限制，以迫使计算机显示一系列相对较短的发夹结构。该程序的另一个版本被用来评估“切除”已经高置信度选择的局部结构后的长程相互作用的可能性。在图4中，由较粗的线绘制的区域是那些为16S核糖体RNA二级结构的暂定模型选择的局部结构和长程相互作用。其中一些结构也可以在最近提出的大肠杆菌16S核糖体RNA（21,22）的两个二级结构中发现。计算机模型显示的其他碱基配对区域被拒绝，由于来自蛋白质S4-RNA结合位点的拓扑研究的特定证据。例如，核苷酸567到574被认为与位于16S核糖体RNA中心的序列（核苷酸818到897）相互作用，因此不能与核苷酸22到27碱基配对。

判断一个指定碱基配对区域是否有效的一个重要标准是检查该区域的二级结构是否保留在其他相关或无关物种或生物体的16S核糖体RNA分子中（21,23）。通常，补偿性碱基变化成对发生，因此碱基配对得以保留，从而允许保留相同或相似的二级结构模序（motif）。例如，在普通变形杆菌中，核苷酸122到141和核苷酸220到237（见图4）的长程碱基配对区域通过补偿性碱基变化而被保留（P.Carbon和C.Ehresmann，个人通信）。因此，16S核糖体RNA的最终模型中包含了这种自由能非常低的结构（-30.5 kcal/mole）。尝试着将一些系统发育信息直接纳入到程序中。我们重新折叠了大肠杆菌16S核糖体RNA的前574个核苷酸，对通过补偿性变化保存在玉米（Zea mays）叶绿体（23）的亲缘关系近的16S核糖体RNA序列中的每个碱基对给予-3.0 kcal/mole能量加成（bonus）。计算机生成的模型显示了一些额外的拟序结构（coherent structures），但是图4中的其他区域丢失了。显然，这里需要更复杂的方法。但是，该示例强调了该算法如何轻松地适应专门的搜索。

图5显示了大肠杆菌16S核糖体RNA的核苷酸571到765的预测结构。大约80%的所示碱基配对通过了严格的审查，并被最终模型保留下来。该图显示了计算机输出的形式，并包含通过酶促研究指定为“可及”的核苷酸碱基对的示例。

D. 结论（Conclusions）

纯技术方面，没有太大的改进空间。我们算法的计算时间与核苷酸数的立方成正比，并且除了比例常数的减少外，这种性能无法进一步提高。可用的计算机存储限制了我们的算法最多可折叠800个核苷酸。忽略分叉环失稳效应的程序可以更有效地重写以折叠多达1000个核苷酸。然而，正如我们刚刚看到的，仅仅折叠600个核苷酸就有问题。大肠杆菌16S核糖体RNA的前574个核苷酸的折叠在技术上没有问题。在我们的电脑上只需244秒，是可以负担得起的。尽管如此，所得到的结构并不完全被我们所接受，尽管它具有最小的自由能，并且符合所有关于化学反应性和酶敏感性的证据。我们不同意那些放弃所有复杂折叠算法，只在调查中使用最基本的计算机搜索碱基配对的人。我们的方法已被证明是研究大肠杆菌16S核糖体RNA二级结构的一个有价值并且灵活的工具。由于我们可以快速轻松地测试相当长的片段或相关子序列的众多折叠可能性，因此消除了许多繁琐的工作。

首先需要改进的是热力学计算。最近，Ninio（24）朝着更有意义的能量计算迈出了一步。他的热力学规则被设计成使三叶草模型（cloverleaf model）在能量上最有利于转运RNA的结构计算。之所以这样做，是因为Salser（6）编写的规则不足以预测所有已知tRNA中约一半的三叶草结构。Ninio的规则在应用于较大的RNA序列时是否更有意义是一个开放性问题。无论如何，我们并不认为我们使用的热力学规则是错误的，它可能仅仅是不够充分的。它们构成了折叠分子的必要基础，但不能指望它们本身就足够预测结构了。一个仅仅基于构象规则和热力学的程序本身不会产生具有生物学意义的分子折叠。有太多具有相似能量的不同结构。算法中还必须包含更多不同种类的附加信息。对于16S核糖体RNA，寻找一种有效的方法将其与相关RNA的系统发育比较纳入优化过程将会是很有用的。

原文致谢和参考文献如下：

ACKNOWLEDGEMENTS
We wish to thank Dr. C. Ehresmann for stimulating discussion on RNA
folding and for communicating unpublished results. P.S. thanks NATO for a
long term fellowship.
NRCC number 18755

Technical details on the computer program will be made available from the
first author to whom reprints should also be requested.
Present address :
Laboratoire de Biochimie, Institut de Biologie Mol§culaire et Cellulaire
du CNRS, 15, rue Rene Descartes, 67084 STRASBOURG Cedex - FRANCE.

REFERENCES
1. TINOCO, I., UHLENBECK, O.C. and LEVINE, M.D. (1971) Nature 230, 362-367
2. UHLENBECK, O.C , BORER, P.N., DENGLER, B. and TINOCO, I. (1973) J. Mol.
Biol. 2i» 483-496
3. GRALLA, J. and CROTHERS, D.M. (1973) J. Mol. Biol. 73_, 497-511
4. GRALLA, J. and CROTHERS, D.M. (1973) J. Mol. Biol. 78^ 301-319
5. TINOCO, I. BORER, P.N., DENGLER, B., LEVINE, M.D., UHLENBECK, O.C ,
CROTHERS, D.M. and GRALLA, J. (1973) Nature New Biol. 246, 40-41
6. SALSER, W. (1977) Cold Spring Harbor Symp. Quant. Biol. 42, 985-1002
7. PIPAS, J.M. and McMAHON, J.E. (1975) Proc. Nat. Acad. Sci. USA 7^2, 2017-
2021
8. STUDNICKA, CM., RAHN, CM., CUMMINGS, I.W. and SALSER, W.A. (1978)
Nucleic Acids Res. 5_, 3365-3387
9. NUSSINOV, R., PIECZENIK, C , GRIGGS, J.R. and KLEITMAN, D.J. (1978) SIAM
J. Appl. Math. 3_5, 68-82
10. WATERMAN, M.S. and SMITH, T.F. (1978) Mathematical Biosciences ^2_, 257-266
11. ROGERS, J., CLARKE, P. and SALSER, W. (1979) Nucleic Acids Res. 6_, 3305-
3321
12. CHAPMAN, N.M. and NOLLER, H.F. (1977) J. Mol. Biol. 2°9_. 131-149
13. HERR, W., CHAPMAN, N.M. and NOLLER, H. (1979) J. Mol. Biol. 130, 433-449
14. EHRESMANN, C , STIEGLER, P., CARBON, P., UNGEWICKELL, E. and GARRET, R.A.
(1977) FEBS Lett. 8^, 188-192
15. UNGEWICKELL, E., GARRETT, R.A., EHRESMANN, C , STIEGLER, P. and CARBON, P.
(1977) FEBS Lett. £U_, 193-198
16. EHRESMANN, C , STIEGLER, P., CARBON, P., UNGEWICKELL, E. and GARRETT, R.A.
(1980) Eur. J. Biochem. _[03, 439-446
17. UNGEWICKELL, E., EHRESMANN, C , STIEGLER, P. and GARRETT, R.A. (1975)
Nucleic Acids Res. 2_, 1867-1888
18. RINKE, J., YUKI, A. and BRIMACOMBE, R. (1976) Eur. J. Biochem. jv4, 77-89
19. CARBON, P., EHRESMANN, C , EHRESMANN, B. and EBEL, J.P. (1979) Eur. J.
Biochem. KX), 399-4)0
20. BOYLE, J., ROBILLARD, G.T. and KIM, S-H. (1980) J. Mol. Biol. 139, 601-625
21. WOESE, C.R., MAGRUM, L.J., GUPTA, R., SIEGEL, R.B., STAHL, D.A., KOP, J.,
CRAWFORD, N., BROSIUS, J., GUTTEL, R., HOGAN, J.J. and NOLLER, H.F. (1980)
Nucleic Acids Res. £, 2275-2293
22. GLOTZ, C. and BRIMACOMBE, R. (1980) Nucleic Acids Res. 8^, 2377-2395
23. SCHWARZ, Zs and KOSSEL, H. (1980) Nature ^83, 739-742
24. NINIO, J. (1979) Biochimie 61, 1133-1150.