本文提出一種結合多序列比較,基于堿基對思想將RNA序列進行預處理,并利用動態規劃法尋找特征值進行子序列劃分,運用最小自由能原理進行結構預測的新方法。本文從公共數據庫選取數據,對算法的準確性進行驗證,用堿基對距離、爬山距離及形態學距離對新方法運算結果實行評估,同時將評估結果與多種預測工具的預測效果相比較,分析新方法的優勢及改進方向。
引用本文： 孫迎, 陸宏偉, 桂堅斌, 宋雪坤. RNA二級結構預測方法研究. 生物醫學工程學雜志, 2014, 31(5): 1065-1069. doi: 10.7507/1001-5515.20140200 復制
引言
據研究表明，RNA分子功能與其結構有著密切的聯系，RNA結構的研究已成為熱點[1]。但利用實驗方法獲得RNA二級結構的代價過高，因此，研究者開始借助計算統計方法從理論上預測RNA二級結構[2]。目前，基于不同的預測理論，主要可分為最小自由能算法和比較序列分析方法，相應的預測軟件分別有RNAforester與RNAstructure等。
本文提出一種將序列比較與最小自由能原理相結合的新算法,命名為SFORESTER。我們從公共數據庫選取數據作為測試樣本，對新算法編程并實行測試。實驗結果表明，SFORESTER算法能夠有效預測RNA二級結構。
1 RNA二級結構預測方法
1.1 比較序列分析方法
比較序列分析法原理是基于序列的上下文語義分析，構建同系列的二級結構模型，并經過多次訓練使之更加優化。由于比較序列分析方法的預測結果僅次于實驗測定，故被認為是值得信賴的計算方法。但方法需要一定數量的相關序列樣本[3]。
序列分析法可以按照序列比對與結構預測的先后順序分為三種[4]。第一種為先比對后預測的方法，此方法的代表軟件有RNAforester。另一種為先預測后比對，這時需得到很多次優結構，但不能確定這些結構中包含了多少真實結構。最后一種為結構預測與序列比對同時進行，主要有Sankoff算法[5]等，基于該算法的軟件有Foldalign[6]和Dynalign[7]，兩者都要求限制子結構的大小和形狀。
1.2 最小自由能方法
最小自由能法的核心是：假定RNA結構的自由能具備可加和性且相對獨立性，即單元結構的自由能是所有組成單元的能量之和，同時各個單元間的自由能又為一確定值，這樣可由式(1)計算總能量的最小值[8],即
| $\begin{align} & {{E}_{i,j}}=min[{{E}_{i+1,j-1}}+{{\alpha }_{ij}},min({{E}_{i+k,j}}+{{\beta }_{k}}), \\ & min({{E}_{i,j-k}}+{{\beta }_{k}}),min({{E}_{i+k,j-l}}+{{\gamma }_{k+l}}), \\ & min({{E}_{i+k,j\prime }}+{{E}_{i\prime ,j-l}}+{{\varepsilon }_{k+l+i\prime -j\prime }}),{{\delta }_{j-i}}], \\ \end{align}$ |
式中αij表示i、 j配對時的堆積能，βk、γk、εk、δk則分別表示突環、內環、多分支環和發卡環的能量，計算結果E1,n為該序列的最小自由能，可以通過回溯得到相應的二級結構，基于最小自由能算法的軟件有RNAstructure、Mfold等。
1.3 RNA二級結構預測精度的評價
有很多方法可以度量預測的準確率，這些度量方法包括敏感性(sensitivity) 、特異性(specificity) 、馬休茲相互作用系數(Matthews correlation coefficient)等參數；另外，還包括假陰性(false negative) 和假陽性(false positive) 等概念[9]。為了更直觀地反映精度，本文采用堿基配對距離(base pair distance)、爬山距離(mountain distance)、形態學距離(morphological distance) 三個參數來評估精度。
堿基對距離是用來衡量RNA二級結構構象之間相似程度的一種方法。其原理為，設有一條長為n堿基的RNA序列，序列的二級結構通常可以表示為一個對稱矩陣Iij，1≤j＜i≤n，其中，當且僅當第i個堿基和第j個堿基配對時，Iij=1；不配對則Iij=0。對于兩個RNA二級結構構象I1=Iij1和I2=Iij2，可以定義堿基對距離如下：
| $D({{I}_{1}},{{I}_{2}})=\sum\limits_{1\le i ＜j\le n}{\left| {{I}_{ij}}^{1}-{{I}_{ij}}^{2} \right|}$ |
堿基對距離可以衡量兩個RNA二級結構構象I1和I2中不同堿基對的個數。堿基對距離計算較為簡單，適用于大規模的RNA二級結構比較分析。
爬山距離同樣可以用作二級結構構象測距。假設集合Sn是設定長度為n堿基的RNA序列的二級結構構象，每一個構象用“點-括號”對方式表示。則對于每一個S∈Sn，定義轉換fS=[fS(1),fS(2),…,fS(n)]，式中fS(i)為從“點-括號”對中第1個位置到第i個位置上所有左括號的個數與右括號的個數之差，則有如下定義[10]：
| $f{{\prime }_{S}}\left( i \right)=\sum\limits_{k=1}^{i}{({{W}_{S}}\left( k \right),{{d}_{M}}({{S}_{1}},{{S}_{2}})}=\left| f{{\prime }_{{{S}_{1}}-f{{\prime }_{{{S}_{2}}\left| 1 \right|}}}} \right|,$ |
式中dM為距離，WS(k)為：
| ${{W}_{S}}\left( k \right)=\left\{ \begin{align} & 1/\left( l-k \right),\left( k,l \right)\in B \\ & -1/\left( k-l \right),\left( l,k \right)\in {{B}_{S}},k\in \left[ 1,\ldots ,n \right] \\ & 0.otherwise \\ \end{align} \right.$ |
形態學距離的辨識度與爬山距離相近,其原理為：定義一個RNA二級結構構象為堿基對集合S，(i，j)∈S，表示堿基i和j配對，對于任意兩個構象S1，S2，存在形態學距離：
| ${{d}_{MD}}\left( {{S}_{1}},{{S}_{2}} \right)=max\left\{ \begin{align} & d{{\prime }_{MD}}\left( {{S}_{1}},{{S}_{2}} \right) \\ & d{{\prime }_{MD}}({{S}_{2}},{{S}_{1}}) \\ \end{align} \right.$ |
式中
| $d{{\prime }_{M\text{ }D}}({{S}_{1}},{{S}_{2}})=\sum\limits_{({{i}_{1}},{{j}_{1}})\in {{S}_{1}}}{\underset{({{i}_{2}},{{j}_{2}})\in {{S}_{2}}}{\mathop{min}}\,}max~\left\{ \begin{align} & \left| {{i}_{1}}-{{i}_{2}} \right| \\ & \left| {{j}_{2}}-{{j}_{2}} \right| \\ \end{align} \right.$ |
上述三個參數分別可以對堿基對、多分支環以及內環等結構做出精細評價。
2 比較序列與最小自由能結合的預測法
由于比較序列分析和最小自由能法各有長處，因此，能夠結合兩者的優缺點，取長補短，發展聯合算法為眾多研究者所關注,而如何和諧地融合兩種算法，充分利用相互的分析優勢，是這方面研究的一大難點[11]。
本文提出一種啟發式的序列對比預測并結合最小自由能法的預測算法，命名為SFORESTER算法。
2.1 算法的提出
為觀察RNA二級結構的預測，本文將實驗驗證的二級結構[12]全長序列進行堿基標記，并以700 nt左右的長度分區[13]，得到獨立二級結構片段。將二級結構子片段與基于最小自由能原理預測軟件Mfold的預測結果進行配對AUGC比率值計算。設：預測軟件的數據為Value(Mfold)，真實結構中的數據為Value(True)，按式(7)進行計算，得到各項指標差異占片段序列總長的比率：
| $\begin{align} & Value(Result)=[Value(Mfold)-Value(True)]/ \\ & Value(Total\text{ }Length) \\ \end{align}$ |
將Mfold工具預測值與實驗得到的真實值對比，由式(7)得出配對堿基G、A、C、U，堿基對GC、GU占片段序列總長的比例，如表 1所示。
表1
各配對堿基占序列總長的比例(%)
Table1.
The proportion of paired bases accounted in sequence length(%)
各項數值越小，表明預測的二級結構中配對的堿基越接近實測結構的配對數。實驗結果同時表明，預測的準確性隨著預測序列長度的增加逐漸下降。實驗數據中有個別數據差異較大，與我們劃分獨立結構時的分割原則有關。
我們由上述實驗經驗提出，將比較序列分析方法與最小自由能原理相結合，合理分割序列，達到充分利用自由能最小進行結構預測的目的。
2.2 算法的思想
RNA序列中，特定的不同堿基之間可以通過互補配對而形成所謂的堿基對，本文將這種通過自身折疊而形成的互補堿基稱為STACK。因此，STACK是組成RNA二級結構的最基本單元。而整個RNA序列就分為兩類區域：① 非配對區或非結構區，即由不和其它任何堿基配對的堿基所組成的序列片段；② 配對區或結構區，即STACK。
本文算法研究的基本指導思想為分析RNA確定序列中的配對區STACK特性。分析方法如下：
假設序列中有兩個堿基配對(A[i]、A[i]′與B[j]、B[j]′)，如果這兩個配對沒有相交的部分，那么由堿基配對的起始位置P組成以下3種關系[14]:
1) 平行
| $P(A[i])＜P(A[i]′)＜P(B[j])＜P(B[j]′)或P(B[j])＜P(B[j]′)＜P(A[i])＜P(A[i]′) $ |
2) 嵌套
| $P(A[i])＜P(B[j])＜P(B[j]′)＜P(A[i]′)或P(B[j])＜P(A[i])＜P(A[i]′＜P(B[j]′) $ |
3) 偽結
| $P(A[i])＜B[j])＜P(A[i]′)＜P(B[j]′)或P(B[j])＜P(A[i])＜P(B[j]′)＜P(A[i]′) $ |
式(8)~(10)中，A[i]′與B[j]′分別是A[i]與B[i]的堿基配對序列。
如果兩個堿基配對之間有相交的部分，那么它們也存在以下3種關系:
1) 相接
A[i]與B[j]相交，A[i]′與B[j]′相交，且滿足條件:
P(B[j])- P(A[i])= P(A[i]′)- P(B[j]′)
2) 假相接
A[i]與B[j]相交，A[i]′與B[j]′相交，且滿足條件:
P(B[j])- P(A[i])≠ P(A[i]′)- P(B[j]′)
3) 互斥
A[i]與B[j]相交，A[i]′與B[j]′不相交；或者A[i]與B[j]不相交，A[i]′與B[j]′相交。
由上可見，我們只要找到兩條序列上的最佳匹配STACK對(STACK1，STACK2)，其中STACK1是一條序列 SEQ1上的STACK，STACK2是另一條序列SEQ2上的STACK。
如果將SEQ1和SEQ2分別依據STACK1、STACK2部分拆分成若干條子序列，然后對這些子序列進行結構比對，就能保證SEQ1子序列上的STACK不會與STACK1相會交錯，也就不會產生偽結。同樣，SEQ2子序列上的STACK也不會與STACK2成為偽結。這樣，序列上匹配的結構信息就不會出現偽結現象。
在尋找特征值的過程中，序列拆分是一個非常關鍵的步驟。其主要思想就是建立在完整序列的組成特征和序列中最佳匹配STACK對的基礎上，將序列拆分成若干個部分，這樣可以將長序列的結構比對問題轉化為若干個短序列的結構預測。算法不僅能避免由于存在偽結而帶來的負面影響，解決長序列的結構比對問題，還可利用若干個短序列進行有效的最小自由能計算。
2.3 SFORESTER算法
本文提出的新算法首先基于STACK的比較序列方法。算法的第一部分采用動態規劃思想對序列進行全局比對，使用迭代法計算特征值得分，剔除得分低項，保留分值高項，并存入一個得分數組中，根據分值位置回溯尋找比對片段，得到片段上匹配堿基對，這樣保證了子序列不會相互交錯。第二部分根據序列的特征值，對子序列進行折疊，根據最小自由能原理獲取各類結構，最終得到整個條序列的二級結構。
算法的主要步驟是：
(1)檢測RNA多序列比對中或待比對的同源RNA序列集合中所包含的保守STACK；
(2)選取高可信度的保守STACK作為初始對象；
(3)以這些對象的起點設為特征值，把輸入的RNA序列切割成多個子片段，然后利用最小自由能原理對這些子片段分別進行結構預測；
(4)對預測得到的子結構進行連接，從而得到完整的RNA二級結構。
算法流程如圖 1所示。
圖1
SFORESTER算法流程圖
Figure1.
SFORESTER Algorithm flow chart
圖 1中，F子函數為利用多序列比對思想進行特征值計算，具體步驟描述如下：
(1)檢測序列集中所有可能的保守STACK。
(2)利用信噪比度量方法評估檢測到的保守STACK。
(3)設置信噪比閾值W，選取信噪比值大于W的保守STACK作為結構預測特征值依據。
程序流程如圖 2所示。
圖2
F子函數程序流程
Figure2.
F Sub-function flow chart
S子函數的步驟如下：
(1)依據所得特征值，將序列集分割若干子片段。
(2)由最小自由能原理折疊序列，根據F子函數中的特征數來定位，為所求序列劃定一個初步折疊范圍，在算法的循環中將對該范圍做出適當微調；根據堿基配對的原則以及氫鍵的強弱來計算每一條子序列的二級結構。
(3)重復上述步驟，直到所有子片段全部被處理，即得到所有子片段各自包含的二級結構信息。
程序流程如圖 3所示。
圖3
S子函數算法流程圖
Figure3.
S Sub-function flow chart
3 實驗結果
本文選用5s rRNA Archaea內的三組序列作為被比對序列，它們來自三個不同的古細菌：泉古菌門、廣古菌門以及納古菌門。預測軟件有SFORESTER、RNAforester與RNAstructure。
3.1 程序的準確性分析
使用堿基對距離、爬山距離以及形態學距離計算上述算法的準確性，結果與RNAforester、RNAstructure的預測精確度作比較，分析在使用SFOR-ESTER算法后精度是否提高。表 2顯示了3種算法的結果。
由表中計算結果可見，由于使用的基礎理論不同，RNAforester與RNAstructure預測軟件的差異很大。在堿基對方面，RNAforester的準確性較高，并且穩定，而形態學距離上相對較低。SFORESTER算法不論從堿基對方面還是形態學距離的準確性均較高，說明結合最小自由能來折疊序列可以更為精確，在莖環結構支鏈的細節部分有更好的預測。
表2
SFORESTER算法精度與比對
Table2.
SFORESTER algorithm and Accuracy assessment
3.2 結果分析
本文提出將比較序列分析和最小自由能分析方法相結合的SFORESTER算法，利用信噪比動態評估STACK取值，將信噪比設為過程中堿基配對與非配對數之比，通過動態調整比值，可得到各個STACK區域在全局上的特征值，分析這些特征值，可達到在全序列調整各子序列STACK的位置，從全局上優化了比對結果，提高了比對精度。
SFORESTER算法利用二級結構保守性、最小自由能原理配合動態規劃法，多方法結合從整體到細節預測RNA二級結構顯得更為全面。
然而，從SFORESTER算法的執行過程與評估中，也能看到算法的完善方向。首先，SFORESTER算法進行多序列比對與最小自由能下的動態規劃需要額外的時間步長。因此，如何簡化該算法，是日后的完善目標之一。其次，進行更多種類的RNA實驗，由于序列比較算法中需要進行大量已驗知RNA二級結構的統計計算，因此需要一個完善的RNA二級結構數據庫。不僅數量上需要完善，已驗證的RNA二級結構的分類規劃也需要完善[15]。
綜上所述，SFORESTER算法作為新穎的RNA二級結構預測方法，在RNA分子結構的預測上有較好的結果。考慮到不同類型的RNA結果特性不同，所以很難將算法統一化，但可以根據不同的RNA來細分，分析其特性并發展相關算法和軟件，這也是將來預測算法發展的另一個方向。
引言
據研究表明，RNA分子功能與其結構有著密切的聯系，RNA結構的研究已成為熱點[1]。但利用實驗方法獲得RNA二級結構的代價過高，因此，研究者開始借助計算統計方法從理論上預測RNA二級結構[2]。目前，基于不同的預測理論，主要可分為最小自由能算法和比較序列分析方法，相應的預測軟件分別有RNAforester與RNAstructure等。
本文提出一種將序列比較與最小自由能原理相結合的新算法,命名為SFORESTER。我們從公共數據庫選取數據作為測試樣本，對新算法編程并實行測試。實驗結果表明，SFORESTER算法能夠有效預測RNA二級結構。
1 RNA二級結構預測方法
1.1 比較序列分析方法
比較序列分析法原理是基于序列的上下文語義分析，構建同系列的二級結構模型，并經過多次訓練使之更加優化。由于比較序列分析方法的預測結果僅次于實驗測定，故被認為是值得信賴的計算方法。但方法需要一定數量的相關序列樣本[3]。
序列分析法可以按照序列比對與結構預測的先后順序分為三種[4]。第一種為先比對后預測的方法，此方法的代表軟件有RNAforester。另一種為先預測后比對，這時需得到很多次優結構，但不能確定這些結構中包含了多少真實結構。最后一種為結構預測與序列比對同時進行，主要有Sankoff算法[5]等，基于該算法的軟件有Foldalign[6]和Dynalign[7]，兩者都要求限制子結構的大小和形狀。
1.2 最小自由能方法
最小自由能法的核心是：假定RNA結構的自由能具備可加和性且相對獨立性，即單元結構的自由能是所有組成單元的能量之和，同時各個單元間的自由能又為一確定值，這樣可由式(1)計算總能量的最小值[8],即
| $\begin{align} & {{E}_{i,j}}=min[{{E}_{i+1,j-1}}+{{\alpha }_{ij}},min({{E}_{i+k,j}}+{{\beta }_{k}}), \\ & min({{E}_{i,j-k}}+{{\beta }_{k}}),min({{E}_{i+k,j-l}}+{{\gamma }_{k+l}}), \\ & min({{E}_{i+k,j\prime }}+{{E}_{i\prime ,j-l}}+{{\varepsilon }_{k+l+i\prime -j\prime }}),{{\delta }_{j-i}}], \\ \end{align}$ |
式中αij表示i、 j配對時的堆積能，βk、γk、εk、δk則分別表示突環、內環、多分支環和發卡環的能量，計算結果E1,n為該序列的最小自由能，可以通過回溯得到相應的二級結構，基于最小自由能算法的軟件有RNAstructure、Mfold等。
1.3 RNA二級結構預測精度的評價
有很多方法可以度量預測的準確率，這些度量方法包括敏感性(sensitivity) 、特異性(specificity) 、馬休茲相互作用系數(Matthews correlation coefficient)等參數；另外，還包括假陰性(false negative) 和假陽性(false positive) 等概念[9]。為了更直觀地反映精度，本文采用堿基配對距離(base pair distance)、爬山距離(mountain distance)、形態學距離(morphological distance) 三個參數來評估精度。
堿基對距離是用來衡量RNA二級結構構象之間相似程度的一種方法。其原理為，設有一條長為n堿基的RNA序列，序列的二級結構通常可以表示為一個對稱矩陣Iij，1≤j＜i≤n，其中，當且僅當第i個堿基和第j個堿基配對時，Iij=1；不配對則Iij=0。對于兩個RNA二級結構構象I1=Iij1和I2=Iij2，可以定義堿基對距離如下：
| $D({{I}_{1}},{{I}_{2}})=\sum\limits_{1\le i ＜j\le n}{\left| {{I}_{ij}}^{1}-{{I}_{ij}}^{2} \right|}$ |
堿基對距離可以衡量兩個RNA二級結構構象I1和I2中不同堿基對的個數。堿基對距離計算較為簡單，適用于大規模的RNA二級結構比較分析。
爬山距離同樣可以用作二級結構構象測距。假設集合Sn是設定長度為n堿基的RNA序列的二級結構構象，每一個構象用“點-括號”對方式表示。則對于每一個S∈Sn，定義轉換fS=[fS(1),fS(2),…,fS(n)]，式中fS(i)為從“點-括號”對中第1個位置到第i個位置上所有左括號的個數與右括號的個數之差，則有如下定義[10]：
| $f{{\prime }_{S}}\left( i \right)=\sum\limits_{k=1}^{i}{({{W}_{S}}\left( k \right),{{d}_{M}}({{S}_{1}},{{S}_{2}})}=\left| f{{\prime }_{{{S}_{1}}-f{{\prime }_{{{S}_{2}}\left| 1 \right|}}}} \right|,$ |
式中dM為距離，WS(k)為：
| ${{W}_{S}}\left( k \right)=\left\{ \begin{align} & 1/\left( l-k \right),\left( k,l \right)\in B \\ & -1/\left( k-l \right),\left( l,k \right)\in {{B}_{S}},k\in \left[ 1,\ldots ,n \right] \\ & 0.otherwise \\ \end{align} \right.$ |
形態學距離的辨識度與爬山距離相近,其原理為：定義一個RNA二級結構構象為堿基對集合S，(i，j)∈S，表示堿基i和j配對，對于任意兩個構象S1，S2，存在形態學距離：
| ${{d}_{MD}}\left( {{S}_{1}},{{S}_{2}} \right)=max\left\{ \begin{align} & d{{\prime }_{MD}}\left( {{S}_{1}},{{S}_{2}} \right) \\ & d{{\prime }_{MD}}({{S}_{2}},{{S}_{1}}) \\ \end{align} \right.$ |
式中
| $d{{\prime }_{M\text{ }D}}({{S}_{1}},{{S}_{2}})=\sum\limits_{({{i}_{1}},{{j}_{1}})\in {{S}_{1}}}{\underset{({{i}_{2}},{{j}_{2}})\in {{S}_{2}}}{\mathop{min}}\,}max~\left\{ \begin{align} & \left| {{i}_{1}}-{{i}_{2}} \right| \\ & \left| {{j}_{2}}-{{j}_{2}} \right| \\ \end{align} \right.$ |
上述三個參數分別可以對堿基對、多分支環以及內環等結構做出精細評價。
2 比較序列與最小自由能結合的預測法
由于比較序列分析和最小自由能法各有長處，因此，能夠結合兩者的優缺點，取長補短，發展聯合算法為眾多研究者所關注,而如何和諧地融合兩種算法，充分利用相互的分析優勢，是這方面研究的一大難點[11]。
本文提出一種啟發式的序列對比預測并結合最小自由能法的預測算法，命名為SFORESTER算法。
2.1 算法的提出
為觀察RNA二級結構的預測，本文將實驗驗證的二級結構[12]全長序列進行堿基標記，并以700 nt左右的長度分區[13]，得到獨立二級結構片段。將二級結構子片段與基于最小自由能原理預測軟件Mfold的預測結果進行配對AUGC比率值計算。設：預測軟件的數據為Value(Mfold)，真實結構中的數據為Value(True)，按式(7)進行計算，得到各項指標差異占片段序列總長的比率：
| $\begin{align} & Value(Result)=[Value(Mfold)-Value(True)]/ \\ & Value(Total\text{ }Length) \\ \end{align}$ |
將Mfold工具預測值與實驗得到的真實值對比，由式(7)得出配對堿基G、A、C、U，堿基對GC、GU占片段序列總長的比例，如表 1所示。
表1
各配對堿基占序列總長的比例(%)
Table1.
The proportion of paired bases accounted in sequence length(%)
各項數值越小，表明預測的二級結構中配對的堿基越接近實測結構的配對數。實驗結果同時表明，預測的準確性隨著預測序列長度的增加逐漸下降。實驗數據中有個別數據差異較大，與我們劃分獨立結構時的分割原則有關。
我們由上述實驗經驗提出，將比較序列分析方法與最小自由能原理相結合，合理分割序列，達到充分利用自由能最小進行結構預測的目的。
2.2 算法的思想
RNA序列中，特定的不同堿基之間可以通過互補配對而形成所謂的堿基對，本文將這種通過自身折疊而形成的互補堿基稱為STACK。因此，STACK是組成RNA二級結構的最基本單元。而整個RNA序列就分為兩類區域：① 非配對區或非結構區，即由不和其它任何堿基配對的堿基所組成的序列片段；② 配對區或結構區，即STACK。
本文算法研究的基本指導思想為分析RNA確定序列中的配對區STACK特性。分析方法如下：
假設序列中有兩個堿基配對(A[i]、A[i]′與B[j]、B[j]′)，如果這兩個配對沒有相交的部分，那么由堿基配對的起始位置P組成以下3種關系[14]:
1) 平行
| $P(A[i])＜P(A[i]′)＜P(B[j])＜P(B[j]′)或P(B[j])＜P(B[j]′)＜P(A[i])＜P(A[i]′) $ |
2) 嵌套
| $P(A[i])＜P(B[j])＜P(B[j]′)＜P(A[i]′)或P(B[j])＜P(A[i])＜P(A[i]′＜P(B[j]′) $ |
3) 偽結
| $P(A[i])＜B[j])＜P(A[i]′)＜P(B[j]′)或P(B[j])＜P(A[i])＜P(B[j]′)＜P(A[i]′) $ |
式(8)~(10)中，A[i]′與B[j]′分別是A[i]與B[i]的堿基配對序列。
如果兩個堿基配對之間有相交的部分，那么它們也存在以下3種關系:
1) 相接
A[i]與B[j]相交，A[i]′與B[j]′相交，且滿足條件:
P(B[j])- P(A[i])= P(A[i]′)- P(B[j]′)
2) 假相接
A[i]與B[j]相交，A[i]′與B[j]′相交，且滿足條件:
P(B[j])- P(A[i])≠ P(A[i]′)- P(B[j]′)
3) 互斥
A[i]與B[j]相交，A[i]′與B[j]′不相交；或者A[i]與B[j]不相交，A[i]′與B[j]′相交。
由上可見，我們只要找到兩條序列上的最佳匹配STACK對(STACK1，STACK2)，其中STACK1是一條序列 SEQ1上的STACK，STACK2是另一條序列SEQ2上的STACK。
如果將SEQ1和SEQ2分別依據STACK1、STACK2部分拆分成若干條子序列，然后對這些子序列進行結構比對，就能保證SEQ1子序列上的STACK不會與STACK1相會交錯，也就不會產生偽結。同樣，SEQ2子序列上的STACK也不會與STACK2成為偽結。這樣，序列上匹配的結構信息就不會出現偽結現象。
在尋找特征值的過程中，序列拆分是一個非常關鍵的步驟。其主要思想就是建立在完整序列的組成特征和序列中最佳匹配STACK對的基礎上，將序列拆分成若干個部分，這樣可以將長序列的結構比對問題轉化為若干個短序列的結構預測。算法不僅能避免由于存在偽結而帶來的負面影響，解決長序列的結構比對問題，還可利用若干個短序列進行有效的最小自由能計算。
2.3 SFORESTER算法
本文提出的新算法首先基于STACK的比較序列方法。算法的第一部分采用動態規劃思想對序列進行全局比對，使用迭代法計算特征值得分，剔除得分低項，保留分值高項，并存入一個得分數組中，根據分值位置回溯尋找比對片段，得到片段上匹配堿基對，這樣保證了子序列不會相互交錯。第二部分根據序列的特征值，對子序列進行折疊，根據最小自由能原理獲取各類結構，最終得到整個條序列的二級結構。
算法的主要步驟是：
(1)檢測RNA多序列比對中或待比對的同源RNA序列集合中所包含的保守STACK；
(2)選取高可信度的保守STACK作為初始對象；
(3)以這些對象的起點設為特征值，把輸入的RNA序列切割成多個子片段，然后利用最小自由能原理對這些子片段分別進行結構預測；
(4)對預測得到的子結構進行連接，從而得到完整的RNA二級結構。
算法流程如圖 1所示。
圖1
SFORESTER算法流程圖
Figure1.
SFORESTER Algorithm flow chart
圖 1中，F子函數為利用多序列比對思想進行特征值計算，具體步驟描述如下：
(1)檢測序列集中所有可能的保守STACK。
(2)利用信噪比度量方法評估檢測到的保守STACK。
(3)設置信噪比閾值W，選取信噪比值大于W的保守STACK作為結構預測特征值依據。
程序流程如圖 2所示。
圖2
F子函數程序流程
Figure2.
F Sub-function flow chart
S子函數的步驟如下：
(1)依據所得特征值，將序列集分割若干子片段。
(2)由最小自由能原理折疊序列，根據F子函數中的特征數來定位，為所求序列劃定一個初步折疊范圍，在算法的循環中將對該范圍做出適當微調；根據堿基配對的原則以及氫鍵的強弱來計算每一條子序列的二級結構。
(3)重復上述步驟，直到所有子片段全部被處理，即得到所有子片段各自包含的二級結構信息。
程序流程如圖 3所示。
圖3
S子函數算法流程圖
Figure3.
S Sub-function flow chart
3 實驗結果
本文選用5s rRNA Archaea內的三組序列作為被比對序列，它們來自三個不同的古細菌：泉古菌門、廣古菌門以及納古菌門。預測軟件有SFORESTER、RNAforester與RNAstructure。
3.1 程序的準確性分析
使用堿基對距離、爬山距離以及形態學距離計算上述算法的準確性，結果與RNAforester、RNAstructure的預測精確度作比較，分析在使用SFOR-ESTER算法后精度是否提高。表 2顯示了3種算法的結果。
由表中計算結果可見，由于使用的基礎理論不同，RNAforester與RNAstructure預測軟件的差異很大。在堿基對方面，RNAforester的準確性較高，并且穩定，而形態學距離上相對較低。SFORESTER算法不論從堿基對方面還是形態學距離的準確性均較高，說明結合最小自由能來折疊序列可以更為精確，在莖環結構支鏈的細節部分有更好的預測。
表2
SFORESTER算法精度與比對
Table2.
SFORESTER algorithm and Accuracy assessment
3.2 結果分析
本文提出將比較序列分析和最小自由能分析方法相結合的SFORESTER算法，利用信噪比動態評估STACK取值，將信噪比設為過程中堿基配對與非配對數之比，通過動態調整比值，可得到各個STACK區域在全局上的特征值，分析這些特征值，可達到在全序列調整各子序列STACK的位置，從全局上優化了比對結果，提高了比對精度。
SFORESTER算法利用二級結構保守性、最小自由能原理配合動態規劃法，多方法結合從整體到細節預測RNA二級結構顯得更為全面。
然而，從SFORESTER算法的執行過程與評估中，也能看到算法的完善方向。首先，SFORESTER算法進行多序列比對與最小自由能下的動態規劃需要額外的時間步長。因此，如何簡化該算法，是日后的完善目標之一。其次，進行更多種類的RNA實驗，由于序列比較算法中需要進行大量已驗知RNA二級結構的統計計算，因此需要一個完善的RNA二級結構數據庫。不僅數量上需要完善，已驗證的RNA二級結構的分類規劃也需要完善[15]。
綜上所述，SFORESTER算法作為新穎的RNA二級結構預測方法，在RNA分子結構的預測上有較好的結果。考慮到不同類型的RNA結果特性不同，所以很難將算法統一化，但可以根據不同的RNA來細分，分析其特性并發展相關算法和軟件，這也是將來預測算法發展的另一個方向。
