10.3 共分散モデル：SCFGに基づくRNAのプロファイル
進化的に関係したRNAのファミリー，例えば転移RNAかグループI触媒イントロンのようなもの，があるとする．これらは共通の二次構造といくつかの一次配列モチーフを持つから，それを用いてデータベースから相同性のあるRNAを検索したい．第5章では，タンパク質とDNA配列ファミリーの共通部分をモデル化するのに，HMMに基づくプロファイルを用いた．しかし第9章で示したように，HMMは一次構造のモデルであってRNA二次構造の制約を効果的に扱えない．この節では共分散モデル(CM:covariance model)と呼ばれるSCFGに基づくRNA構造プロファイルについて述べる．これはHMMに対するプロファイルHMMに相当する，SCFGに於ける同様なモデルである．プロファイルHMMは，配列のマルチプルアラインメントのモデル化に適した線形のHMM構造を規定しているが，CMはRNAの共通二次構造のモデル化に適した木の形のSCFG構造を規定している．
ここでは[Eddy, S. R.& Durbin, R.1994]による共分散モデルに従うが，同じ一般的なアイデアとアルゴリズムは[Sakakibara, Y. et al.1994]によって同時期に独立に開発されたSCFGに基づくRNAモデルでも同様である．
CMは詳細で複雑な確率モデルである．最初に，小さなRNAのアラインメントのより単純なモデルで直感的に説明する．
ギャップなしRNAアラインメントのSCFGモデル
図10.11は，RNA共通構造と，共通構造に一致するRNAファミリーのギャップなしアラインメントの例である．このSCFGに基づくモデルによるマルチプルアラインメントを記述するには，様々なタイプの二次構造と配列の要素を生成するそれぞれ異なった非終端記号が必要である．
図10.11 挿入や欠失のないRNAファミリー共通構造を上に示す．同一の構造を持つ5個の異なった生物種の配列例のマルチプルアラインメントを構造の下に示す．塩基対の位置を箱で囲み，塩基対のパターンを線で結んだ．アラインメントの最後の行は，RNA構造アラインメントの注釈付けに用いるフォーマットで共通構造表現である．
塩基対を構成する列は，塩基対の両方の塩基を発生する非終端記号を出力することで，ペアワイズにモデル化されている．
一本鎖の列は，可能な限り左から順に非終端記号を出力するようにモデル化されている．ステムの\(3’\)側のバルジや内側ループについては，右から順に非終端記号を出力することが必要になる．枝分かれ非終端記号は，複数のステムとマルチループに分かれるときに使われる．初期非終端記号と枝分かれの直後に生まれる子供として，特殊な開始非終端記号を定義する1．
また，確率1で\(\epsilon\)を発生して生成を終了させる，特殊な終了非終端記号を定義する．これらの状態の生成規則を以下にまとめる．$W$は6個の状態の任意の1個を表す一般化非終端記号である．
\[
\begin{eqnarray}
\begin{array}{llll}
P & \rightarrow & aWb & {\rm ペアワイズ(16ペアの出力確率)}\\
L & \rightarrow & aW & {\rm 左側出力(4個の単独出力確率)}\\
R & \rightarrow & Wa & {\rm 右側出力(4個の単独出力確率)}\\
B & \rightarrow & SS & {\rm 枝分かれ(確率1)}\\
S & \rightarrow & W & {\rm 開始(確率1)}\\
E & \rightarrow & \epsilon & {\rm 終了(確率1)}
\end{array}
\end{eqnarray}
\]
下に詳細を示すように，図10.11の構造は1個のSCFGに集約される．分かりやすいように，それぞれの非終端記号に対して可能な生成を1個だけ(図10.11の構造とヒトの配列に対応する生成)を示した．ペアワイズの生成は16種類の生成とその確率が可能な16種類のペアに対して存在する．右側と左側の生成は4個である．
\[
\begin{array}{lllllllll}
& & & {\bf Stem1} & & {\bf Stem2} & \\
S_1 & \rightarrow & L_2 \ldots & S_5 \rightarrow & P_6 & S_{15} \rightarrow & L_{16}\\
L_2 & \rightarrow & aL_3 \ldots & P_6 \rightarrow & gP_7c \ldots & L_{16} \rightarrow & uP_{17}\ldots\\
L_3 & \rightarrow & aB_4 \ldots & P_7 \rightarrow & aR_8u \ldots & P_{17} \rightarrow & gP_{18}c\ldots\\
B_4 & \rightarrow & S_{5}S_{15} & R_8 \rightarrow & P_9a \dots & P_{18} \rightarrow & gL_{19}c\ldots\\
& & & P_9 \rightarrow & cL_{10}g \ldots & L_{19} \rightarrow & cP_{20}\ldots\\
& & & L_{10} \rightarrow & uL_{11} \ldots & P_{20} \rightarrow & gL_{21}c\ldots\\
& & & L_{11} \rightarrow & uL_{12} \ldots & L_{21} \rightarrow & aL_{22}\ldots\\
& & & L_{12} \rightarrow & cL_{13} \ldots & L_{22} \rightarrow & cL_{23}\ldots\\
& & & L_{13} \rightarrow & gL_{14} \ldots & L_{23} \rightarrow & aL_{24}\ldots\\
& & & L_{14} \rightarrow & \epsilon & L_{24} \rightarrow & \epsilon
\end{array}
\]
このモデルには重要な性質がある．それは非終端記号が木形状につながっていることである．SCFGの木はRNAの構造と構文解析木の構造を正確に反映している(これはSCFGの必要な性質ではない．例えば，先に見たRNA折れ畳みSCFGではそうなっていない)．これにより，モデル化しようとしているRNAファミリーの構造を直感的かつ簡潔に反映したSCFGの便利なグラフィック表現が可能である．上の文法からも判るように，単純なRNAのSCFGでさえ生成規則の形に書くことはとても煩わしい．同じSCFGのグラフィック表現を図10.12に示す．
モデルには全部で24個の非終端記号があり，24塩基のRNAアラインメントをモデル化している．この2個の数は必ずしも正確に同じである必要はないが，非終端記号の数は大まかに言って，アラインメントの長さとともに線形に増加する．それぞれのペアに1個の非終端記号が，それぞれの1本鎖塩基に1個の非終端記号が必要であり，非終端記号$B$,$S$,$E$を組み合わせてモデルを完成させる．
生成規則のリストと図10.12のグラフィックモデルとの間には重要な違いがある．SCFGは典型的には状態遷移で記号を出力する(すなわち，\(W_1 \rightarrow aW_2b\)というように記号出力と新たな非終端記号への移動を同時に行う)．遷移と出力を分離し，先立つ状態遷移と独立に，状態から記号を出力するようにすることもできる．この状態出力(ミーリー機械)と記号出力(ムーア機械)の区別は第9章で議論した．これまで述べたほとんどのHMMはプロファイルHMMを含めてムーア機械を用いてきた．第5章のプロファイルHMMのSCFGに基づく拡張として，CMでも状態出力の形式を用いることにする．同様に，CMの非終端記号を「状態」とも呼ぶ．そうすると，図10.12のギャップなしSCFGには，ペアワイズ状態毎に16個の出力確率と左側・右側状態毎に4個の出力確率があり，すべての遷移確率は1である(ギャップなしモデルでは代替経路はない．挿入・欠失を許すモデルを開発するときには，遷移確率はより興味深いものとなる．
RNA構造アラインメントをSCFGに置き換えるには曖昧性がある．例えば，へアピンループは右側状態から右から左へと生成することができる．左側状態はプロファイルHMMの扱いと同様なので，可能な限り左側状態を使うことにする．
挿入や欠失がない限り，図10.12のモデルは図10.11のRNAファミリーの合理的なモデルである．このモデルは挿入／欠失状態がなく一致状態のみからなるHMMと比較可能である．このモデルは，完全な確率モデルとして，RNAのデータベース検索に用いることができる2．
しかし，多くの相同的な構造を見逃すだろう．ここで，一致状態だけからなるHMM(すなわちギャップなし重み行列)を挿入・欠失状態でプロファイルHMMに拡張したように，挿入や欠失に耐えられるようにギャップなしSCFGを拡張したCMモデルの話題に移ることにしよう．
図10.12 ギャップなしのRNA SCFGのグラフィック表現の例を左図に示す．\(P\)とラベルされた箱は16種のペアワイズの生成確率を表す．\(L\)および\(R\)とラベルされた箱は4種の左側生成と右側生成を表す．\(S\),\(B\),\(E\)とラベルされた箱は，開始，二股分かれ，終了の非終端記号を表す．SCFGの木構造により近い表現でRNA共通構造を書き換えて示す(中央)．この構造の構文解析木を，SCFGの状態にRNA塩基を割り当てて示す．
[問題10.10]
ギャップなしRNAモデルの生成規則のリストを，HMMと同様に，前の状態に独立に記号が出力されるように書き換えよ．
これは図10.12のSCFGのグラフィック表現に対応する形式的な確率変形文法である．
共分散モデルの設計
CM設計のゴールははっきりしているが，どのようにゴールに到達するかは自明ではない．最初に，ギャップなしRNAモデルで議論した構文解析木と同様に，RNA共通構造のCMを構築する．次に，CMは任意の位置に任意の長さの挿入や欠失を許すようにする．プロファイルHMMが共通部分に対して挿入と欠失を扱ったときと同じ戦略を，
CMは用いる．プロファイルHMMは，それぞれの位置のマルチプルアラインメントのモデル化に，3種類の状態(一致，挿入，欠失)の組を繰り返し用いることを思い出せ．
プロファイルHMMは，ギャップなし共通部分モデルの定型的な展開と考えることができる．プロファイルHMMでは，ギャップなしモデルのすべての一致状態は，一致，欠失，挿入状態に展開される．同様に，CMはギャップなし共通部分モデルを個別状態の定型的なパターンに展開する．このモデル構造の繰り返し部分をノード(node)と呼ぶことにする．プロファイルHMMは，1種類しかない3状態ノードの，線形の列である．CMはノードの枝分かれした木で，ノードにはいくつかのタイプがあり，異なった数の状態がある．
一本鎖の共通部分の左側ノードを，HMMノードと同様に一致・挿入・欠失状態に展開し，右側単独塩基ノードも展開する．このようにして，左側ノードは\(ML,IL,D\)の3状態となり，右側ノードは\(MR,IR,D\)の3状態となる．
ペアワイズノードを展開するときは，考慮すべき挿入・欠失の可能性がいくつかある．欠失では，塩基対の両塩基が取り除かれるか，\(5’\)か\(3’\)側の一方のみ削除して相手の塩基をバルジとして残すかである．塩基対ステムでの挿入は\((5’\)側，\(3’\)側，両方のどれかで起る．CMはペアワイズノードを6個の状態，\(MP\)状態，\(D\)状態(塩基対を完全に欠失させるもの)，\(ML\)状態と\(MR\)状態\(5’\)と\(3’\)塩基を取り除く単塩基の欠失)，塩基対の\(5’\)側か\(3’\)側で挿入を許す\(IL\)状態・\(IR\)状態，に展開する．
根の開始ノードは，開始状態と\(5’\)・\(3’\)側の挿入状態\(IL\)・\(IR\)に展開される．二股枝分かれの左側の子供の開始ノードは，
1個の\(S\)状態だけに展開される．二股枝分かれの右側の子供の開始ノードは，1個の\(S\)状態と左挿入\(IL\)状態に展開される．この挿入状態の配置は，曖昧に割り当てられた任意の位置の挿入を保証している3．
共通木の二股枝分かれノードと終了ノードはCMでは単に\(B\)状態と\(E\)状態となる．
ここで状態は状態遷移で結ばれる．プロファイルHMMと同様に，状態は，そのノードのすべての挿入状態と，次のノードのすべての非挿入状態と結ばれる．1塩基を越える挿入を許すために，挿入状態は自分自身への状態遷移を持っている．ペアワイズノードでは，\(IL\)は\(IR\)に結ばれているが，逆向きは結ばれず，挿入がモデルに沿って単一の経路を曖昧性なく割り当てられるようになっている．図10.13に状態遷移の結合関係をまとめて図示する．
完全なCMは状態の有向グラフで，背後の共通構造木に従って描かれている．CMの「主軸」は共通構造木と完全に一致するが，CMは欠失と挿入を扱う代替状態を通る経路を許す．
これは可能なRNA構造プロファイルのうちたったひとつの設計である．ステムにおける挿入は一本鎖ではなく塩基対に生じる(すなわち，ステムの長さが変化する)．
従って，ペアワイズ挿入(\(IP\))状態を，ペアワイズノードに含めることもできたのである．\(ML\)・\(MR\)状態をペアワイズノードから取り除き，その代わりに，ペア中の単独塩基の欠失(バルジを生じる)を，ペアの完全な欠失の後で\(IL\)・\(IR\)状態によってバルジが挿入されるという風にモデル化することもできたのである．洗練された設計では，RNAの長い挿入はしばしば構造を持っている事実をモデル化することも，試みられるかも知れない．無限の計算機資源が与えられれば，それぞれの挿入状態を，最初の節で述べた一般化したSCFGのRNA折れ畳みモデルで置き換えることさえも考えることができる．
図10.12 CMの状態(小さな四角の箱)はRNA共通構造の木(大きな箱)のノードに従ってグループ化されている．
ここの例は厳密に仮想的なもので，8種類すべてのCMノードを同じ絵に置くためだけに設計した．状態遷移は矢印で表した．共通経路の「主軸」は太い矢印で描かれている．この主軸は，共通の木自身である．波線は，その下にモデルの示されていないノードがあることを示す．
RNAアラインメントからのCMの構築
RNA配列のアラインメント，二次構造の注釈，どの列が挿入でどの列が共通部分なのかの注釈が与えられると，構造があらかじめ備わった状態でCMを詳細に定義できる．アラインメントの非挿入の列の構造注釈を使い，最初に共通構造の木が構築される．そしてこの木のノードはCM状態で埋められ，状態は上で議論したように状態遷移で結ばれる．アラインメントにおける列の木のノードに対するこの割り当てに基づき，個々の記号が個々のCM状態に割り当てられ，アラインメントのそれぞれの配列に一意なCM構文解析木が割り当てられる．これらの構文解析木における出力と遷移の回数が数えられる．これらの観測された度数は，CMの記号出力確率と遷移確率の推定に使われる．通常は事前分布としてDirichlet分布もしくは混合Dirichlet分布と事後平均推定も同時に用いられる．CMの構造はRNA共通構造を直接反映するから，手続きに曖昧性はなく，非常に高速である．
どの列が共通構造か注釈付けされていなければ，50%以上ギャップ記号がある列を挿入状態と見なすといった，簡単なヒューリスティックを使うことができる．もし共通構造自身が分からなければ，これはより難しい問題であり，より詳しくは章の後半で議論する．
このようにして，構造が注釈付けされたRNA配列のマルチプルアラインメントから，CMを構築することができる．ここで，アラインメントアルゴリズムの記述に話題を移そう．
アラインメントアルゴリズムは，データベース検索，構造に基づく新しいアラインメント，当初アラインメントされず構造化もされていない配列からのCMの学習，
に必要である．
CMアラインメントアルゴリズム
第9章ではチョムスキー標準形のSCFGに適用する一般的なSCFGアルゴリズムを与えた．しかし前に述べたように，RNA解析のためにチョムスキー標準形のSCFGを用いることは煩わしい．チョムスキー標準形は，\(W \rightarrow cWg\)のようなペアワイズの生成を許さない．従って，CMに特有な同様のアラインメントアルゴリズムの集合を与える．
記法
CMは\(W_1,\ldots,W_M\)と書かれる\(M\)個の異なった状態(非終端記号)から成っている．
\(v,y,z\)を状態\(W_v,W_y,W_z\)の添字とする．
\(P,L,R,D,S,B,E\)は，ペアワイズの出力，左側(\(5’\))出力，右側(\(3’\))出力，欠失，開始，二股枝分かれ，終了の7種類の状態を表すラベルとする．
\(W_1\)はCM全体の開始状態である．
CMは多岐RNA二次構造を反映した多岐モデルであるから，通常1個以上の終了状態がある．
7種の状態のタイプには，表10.01に示すように，記号出力と状態遷移の確率が割り当てられている．
状態\(v\)で左側及び右側に記号を出力する回数を\(\Delta_v^L\)と\(\Delta_v^R\)と定義する．
これらを導入することにより，アルゴリズムの記述が単純になる．第6章でのSankoff-Cedergrenアルゴリズム(N次元DPのアルゴリズム)記述でも，同様にアルゴリズムを単純化する記法を用いた．
記法と実装の便利のために，CMのそれぞれの状態は付加的な情報を持っている．7個の形の生成規則を表す
\(P,L,R,D,S,B,E\)の中から値をとる状態のタイプ(state type)を，\(s_v\)とする．
\({\cal C}_v\)を，\(W_v\)から状態遷移できる状態\(W_y\)の添字のリストで表される状態の子供の集合とする．
\({\cal P}_v\)を，\(W_v\)に状態遷移できる状態\(W_y\)の添字のリストで表される状態の親の集合とする．
CMでは二股枝分かれ(\(B\))状態は特別に扱われる．二股枝分かれ状態\(W_v\)は，選ばれた2個の\(S\)状態\(W_y\)と\(W_z\)に
常に確率1で遷移する．\(B\)状態の子供のリスト\({\cal C}_v\)は2個の\(S\)状態のペア\((y,z)\)である．
\(S\)状態である2個の子供の親のリスト\({\cal P}_y\)と\({\cal P}_z\)は\({v}\)である．なぜなら，それぞれの状態に遷移するのは\(W_v\)だけだからである．
二股枝分かれで遷移がひとつしか選べないことはチョムスキー標準形とは違うが，その様子は二股枝分かれ規則\(W_v \rightarrow W_yW_z\)の\(y\),\(z\)の確率的な選択によって記述できる．RNAモデルでは，二股枝分かれ状態は，構造上のマルチループあるいは多重ステムを記述する場合にだけ必要である．モデルの大部分は\(P,L,R\)状態から成る．二股枝分かれのこの制限は，CMアルゴリズムの計算複雑性を大幅に減少させている．
開始(\(S\))と欠失(\(D\))の状態は，アラインメントアルゴリズムと同様に扱われる．唯一の違いは，構造を考慮することである．開始状態は二股枝分かれの直接の子供としてと，根状態\(W_1\)にだけ現れる．欠失状態はCMの\(P,L,R\)ノードの中に現れる．
CMには3個の付加的な制限がある．第一に，それぞれの状態は1個の生成規則のタイプを用いる(\(s_v\)は状態と生成規則の両方のタイプを示す)．
第二に，HMMプロファイルと同様，状態は全結合されていない．\({\cal C}_v\)の結合状態の数は定数であって，状態の総数\(M\)には依存しない．この制限が，チョムスキー標準形のSCFGと比べて，アラインメントアルゴリズムの複雑性をさらに減少させている．
最後に，重要な制限を課す．すべての\(y \in {\cal C}_v\)に対して\(y>v\)となるように，挿入状態を除く状態はナンバーを振られている．
挿入状態については，\(y \in {\cal C}_v\)に対して\(y\geq v\)である．この条件は出力しない状態(\(S,D,B\))について重要であり，出力のない循環が起らないことを保証する．HMMに対する欠失状態に関する同様の制限については第3章で述べた．
ここで，RNAのCMを扱うために重要なアルゴリズムを検討しよう．
\[
\begin{array}{llllll}
状態 & 生成 & \Delta_v^L & \Delta_v^R & 出力 & 遷移\\
P & W_v \rightarrow x_iW_yx_j & 1 & 1 & e_v(x_i,x_j) & t_v(y)\\
L & W_v \rightarrow x_iW_y & 1 & 0 & e_v(x_i) & t_v(y)\\
R & W_v \rightarrow W_yx_j & 0 & 1 &e_v(x_j)b & t_v(y)\\
D & W_v \rightarrow W_y & 0 & 0 &1 & t_v(y)\\
S & W_v \rightarrow W_y & 0 & 0 &1 & t_v(y)\\
B & W_v \rightarrow W_yW_z & 0 & 0 &1 & 1\\
E & W_v \rightarrow \epsilon & 0 & 0 &1 & 1\\
\end{array}
\]
スコア付け：内側アルゴリズム
観測された\(L\)個の記号からなるRNA配列\(x=x_1,\ldots,x_i,\ldots,x_j,\ldots,x_L\)があるとする．与えられたCM，$¥theta$に対する配列の尤度を，\(x\)のすべての可能な構造について合計したもの，\(P(x|\theta)\)を計算する，
スコア付け問題を最初に考えよう．この確率は内側アルゴリズムを用いて計算される．
内側アルゴリズムは，3次元のDP行列を\(\alpha_v(i,j)\)で再帰的に書き埋めていく．\(\alpha_v(i,j)\)は，状態\(v\)を根にもつ部分配列\(x_i,\ldots,x_j\)のすべての構文解析部分木について合計された確率である．\(\alpha_v(i+1,i)\)は，長さ0のヌル部分配列の確率である．
出力のない状態\(D,S,B\)が存在するので，これは境界条件に含まれていなければならない．記法を簡単にするために，\(e_v(x_i,x_j)\)をすべての出力確率を表すのに用いる．
\(L\)状態については\(e_v(x_i,x_j)=e_v(x_i)\)，\(R\)状態については\(e_v(x_i,x_j)=e_v(x_j)\)，出力にない状態については\(e_v(x_i,x_j)=1\)である．
CMの内側アルゴリズム
初期化：for \(j=0\) to \(L, v=M\) to \(1\);
\[
\alpha_v(j+1,j) =
\left\{
\begin{array}{ll}
s_v = E: & 1;\\
s_v \in S,D: & \sum_{y \in {\cal C}_v}t_v(y)\alpha_y(j+1,j);\\
s_v =B: & \alpha_y(j+1,j)\alpha_z(j+1,j);\\
s_v \in P,L,R & 0.\\
\end{array}\right.\]
再帰： for \(j=1\) to \(L\), \(i=j\) to \(1\), \(v=M\) to \(1\);
\[
\begin{eqnarray}
\alpha_v(i,j) =
\left\{
\begin{array}{ll}
s_v = E: & 0;\\
s_v = P,j=i: & 0;\\
s_v = B: & \sum_{k=i-1}^{j}\alpha_y(i,k)\alpha_z(k+1,j);\\
{\rm otherwize:}\\
& e_v(x_i,x_j)\sum_{y \in {\cal C}_v}t_v(y)
\alpha_y(i+\Delta_v^L,j-\Delta_v^R) \\
\end{array}
\right.
\end{eqnarray}
\]
上のアルゴリズムが終了したとき，\(\alpha_1(1,L)\)が\(P(x|\theta)\)になる．もし\(b\)個の二股枝分かれ状態と\(a\)個のその他の状態があれば(\(M=a+b\))，このアルゴリズムの複雑性のオーダーはメモリが\(O(L^2M)\)，時間が\(O(aML^2+bML^3)\)である．
CMの外側アルゴリズム
次の節で説明する内側・外側パラメータ推定のためには，外側アルゴリズムが必要である．外側アルゴリズムは，部分配列\(x_i,\ldots,x_j\)を除く配列全体を生成する状態$v$を根にもつすべての構文解析木の確率の合計\(\beta_v(i,j)\)を計算する．
外側アルゴリズムでは内側の\(\alpha\)項を最初に計算しておく必要がある4．
3次元DP行列のすべての要素をゼロに初期化した後，外側アルゴリズムは，以下のようになる．
CMの外側アルゴリズム
初期化：
\[
\beta_1(1,L) = 1
\]
再帰： for \(i=1\) to \(L+1\), \(j=L\) to \(i-1\), \(v=2\) to \(M\);
\[
\begin{eqnarray}
\beta_v(i,j) =
\left\{
\begin{array}{l}
\mbox{ for } s_v= S,P_v=y,C_y={v,z}:\\
{\rm for }s_v= S, P_v=y, C_y={v,z}:\\
{\rm }\sum_{k=j}^L\beta_y(i,k)\alpha_z(j+1,k);\\
{\rm for }s_v= S, P_v=y, C_y={z,v}:\\
{\rm }\sum_{k=1}^i\beta_y(k,j)\alpha_z(k,i-1);\\
{\rm for }s_v \in P,L,R,D,B,E:\\
{\rm }\sum_{y \in P_v}
e_y(x_{i-\Delta_y^L},x_{j+\Delta_y^R})
t_y(v)\beta_y(i-\Delta_y^L,j+\Delta_y^R).
\end{array}
\right.
\end{eqnarray}
\]
\(\lhd\)
外側アルゴリズムのメモリと時間の複雑性は，内側アルゴリズムと同一である．
CMパラメータの内側・外側期待値最大化
もしモデルの構造が分かっていて(つまりファミリーのRNA共通構造が分かっていて)，確率パラメータが不明なら，内側・外側アルゴリズムと呼ばれるEMアルゴリズムを用いて，確率パラメータをアラインメントされていない学習用RNAから推定できる．実際，このアルゴリズムはめったに使われることがない．ほとんどすべてのRNA共通構造は，比較配列解析を用いてマルチプルアラインメントから導かれる．構造の注釈が付いたマルチプルアラインメントから，直接パラメータの値付きのCMを導く方法に付いては既に述べた．
しかし，内側・外側アルゴリズムのCM版を完全性のために与える．他の方法で共通構造が得られる場合を想定できる．マルチプルアラインメントが完全に正しいとは仮定したくない場合もあり，そのような場合には，アラインメントの度数のデータをそのまま使いたくないから，内側・外側アルゴリズムがより適切となる．
配列\(x\)の導出で，状態\(v\)を部分列\(x_i,\ldots,x_j\)に対して使う確率は，\(\frac{1}{P(x|\theta)}\alpha_v(i,j)\beta_v(i,j)\)である．これらの項を様々な方向に足し合わせて，数多くの有用な確率や期待度数を得ることができる．その中には，チョムスキー標準形のSCFGでやったと同様に，出力確率及び遷移確率のパラメータのEM再推定に必要な期待度数も含まれる．
状態\(v\)が単一配列\(x\)に使われる回数の期待値は，
\[
c(v {\rm used}) = \frac{1}{P(x|\theta)}
\sum_{i=1}^{L+1}\sum_{j=i-1}^L\alpha_v(i,j)\beta(i,j)
\]
状態遷移\(t_v(y)\)が単一配列\(x\)に使われる回数の期待値は，
\[
c(v \rightarrow {\rm y used}) =
\frac{1}{P(x|\theta)}
\sum_{i=1}^{L+1}\sum_{j=i-1}^L\beta_v(i,j)e_v(x_i,x_j)
t_v(y)\alpha_y(i+\Delta_v^L,j-\Delta_v^R)
\]
観測された\(N\)個の独立な配列\(x^1,\ldots,x^h,\ldots,x^N\)を単一配列\(x\)の代わりに学習に使うときは通常RNA配列のファミリーからモデルを学習するときと同様である．それぞれの配列に対して計算された内側，外側変数\(\alpha^h，\beta^h\)を用いてそれぞれの配列の期待度数を合計する．
\[
\begin{eqnarray}
c(v & \rightarrow & {\rm y used}) \\
& = & \sum_{h=1}^N\frac{1}{P(x^h|\theta)}
\sum_{i=1}^{L+1}\sum_{j=i-1}^L
\beta_v^h(i,j)e_v(x_i^h,x_j^h)
t_v(y)\alpha_y^h(i+\Delta_v^L,j-\Delta_v^R)
\end{eqnarray}
\]
このようにして，\(N\)個の学習配列\(x^1,\ldots,x^N\)が与えられたときの，状態\(v\)から状態\(y\)へのCM遷移確率の内側・外側EM再推定の式は以下のようになる5．
\[
\hat{t}v(y) = \frac{\sum{h=1}^N\frac{1}{P(x^h|\theta)}\sum_{i=1}^{L+1}\sum_{j=i-1}^L\beta_v^h(i,j)e_v(x_i^h,x_j^h)t_v(y)\alpha_y^h(i+\Delta_v^L,j-\Delta_v^R)}{\sum_{h=1}^N\frac{1}{P(x^h|\theta)}\sum_{i=1}^{L+1}\sum_{j=i-1}^\alpha_v^h(i,j)\beta_v^h(i,j)}
\]
同様の議論で，状態$v$のCM出力確率の最推定の式は，
以下のようになる．
\(\delta()\)は括弧の中が真の時1で，偽の時0である．
\[
\begin{eqnarray}
{\rm for }s_v=P:\\
\hat{e}v(a,b) & = \frac
{\sum{h=1}^N\frac{1}{P(x^h|\theta)}
\sum_{i=1}^{L}\sum_{j=i}^L
\delta(x_i^h,x_j^h=a,b)
\beta_v^h(i,j)\alpha_v^h(i,j)}
{\sum_{h=1}^N\frac{1}{P(x^h|\theta)}
\sum_{i=1}^{L}\sum_{j=i}^L
\beta_v^h(i,j)\alpha_v^h(i,j)}\\
{\rm for }s_v=L:\\
\hat{e}v(a) & = \frac
{\sum{h=1}^N\frac{1}{P(x^h|\theta)}
\sum_{i=1}^{L}\sum_{j=i}^L
\delta(x_i^h=a)
\beta_v^h(i,j)\alpha_v^h(i,j)}
{\sum_{h=1}^N\frac{1}{P(x^h|\theta)}
\sum_{i=1}^{L}\sum_{j=i}^L
\beta_v^h(i,j)\alpha_v^h(i,j)}\\
{\rm for }s_v=R:\\
\hat{e}v(a) & = \frac
{\sum{h=1}^N\frac{1}{P(x^h|\theta)}
\sum_{i=1}^{L}\sum_{j=i}^L
\delta(x_j^h=a)
\beta_v^h(i,j)\alpha_v^h(i,j)}
{\sum_{h=1}^N\frac{1}{P(x^h|\theta)}
\sum_{i=1}^{L}\sum_{j=i}^L
\beta_v^h(i,j)\alpha_v^h(i,j)}\\
\end{eqnarray}
\]
内側・外側アルゴリズムで計算される量は，
その他の量を推定するのにも使うことができる．
例えば，
単一配列$x$の$x_i$と$x_j$が任意のペアワイズの生成によって
塩基対となる確率は，以下のようになる．
\[
P(x_i,x_j {rm paired})
= \frac{1}{P(x|\theta)}\sum_{v|s_v=P}\alpha_v(i,j)\beta_v(i,j)
\]
データベース検索：CYKアルゴリズム
非常に長い配列(例えば完全ゲノム)を与えられ，RNAモデルに一致する部分配列を見つけることが課題であると仮定しよう．これまでに与えたアルゴリズムは大域アラインメントには適していたが，部分配列のCMに対する局所アラインメントには向いていない．時間及びメモリを，データベースの配列長\(L\)の2乗や3乗のスケールで必要とするようなアルゴリズムは，明らかに不適切である．モデルにアラインメントされる最も長い部分配列の長さを一定値\(D\)に制限し，DP行列の座標変換を採用することによって，データベース検索のための効率的なCYKアルゴリズム(または内側アルゴリズム，あるいは外側アルゴリズム)を実装できる．DP行列は添字\(v,i,j\)の代わりに\(v,j,d\)(\(d\)は部分配列\(i,\ldots,j\)の長さ\(d=j-i+1,d\leq D\))を使うように座標変換される．
位置\(j\)で終わる長さ\(0,\ldots D\)の部分配列に対して，最適アラインメントの行スコアを再帰的に計算することが，この変換された座標では簡単であることを図10.14に示す．
SCFGアラインメントの標準CYKアルゴリズムはモデル$¥theta$に対する配列$S$の最適構文解析$¥hat{¥pi}$の確率の対数$P(S,¥hat{¥pi}|¥theta)$を返す．このスコアはアラインメントされる配列の長さに依存する関数である．従って，このスコアを用いて，データベース検索において異なった長さの重なりあった部分配列の中から，最適な部分配列を選ぶことは困難である．HMMの場合に述べたように，この問題に対する良い解決法は，ランダムな配列の「ヌル」モデルに対する対数オッズ比を計算することである．もしランダムモデルが独立同一分布(iid)で，帰無仮説での配列の尤度が個々の塩基の頻度$f_a$の積となるなら，HMMアラインメントのスコア付けと同様に，対数確率の出力の項を塩基対や単一塩基の対数オッズ比で置き換えて，CYKアルゴリズムが直接対数オッズ比を計算するようにできる6．
下のアルゴリズムでは，対数確率\(\log{e}\)ではなく，対数オッズ比の出力スコア$¥log{¥hat{e}}$を用いる．
\[
\begin{eqnarray}
\begin{array}{llll}
{\rm for }s_v = P: & \log\hat{e}_v(a,b) & = & \log(e_v(a,b)/f_af_b);\\
{\rm for }s_v = L: & \log\hat{e}_v(a,b) & = & \log(e_v(a)/f_a);\\
{\rm for }s_v = R: & \log\hat{e}_v(a,b) & = & \log(e_v(b)/f_b);
\end{array}
\end{eqnarray}
CYKデータベース検索アルゴリズムは以下の通りである．
¥bigskip
¥noindent
{¥bf アルゴリズム： データベース検索のためのCYK}
¥begin{eqnarray}
初期化： & {¥rm for¥ }j=0{¥rm ¥ to¥ }L, v=M {¥rm ¥ to¥ }1:¥¥
¥gamma_v(j,0) = &
¥left¥{
¥begin{array}{l}
{¥rm for¥ }s_v = E: ¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }0;¥¥
{¥rm for¥ }s_v ¥in D,S:¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }¥max_{y ¥in {¥cal C}v}[¥gamma_y(j,0)+¥log{t_v(y)}];¥¥
{¥rm for¥ }s_v = B, {¥cal C}_v = (y,z):¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }¥gamma_y(j,0) + ¥gamma_z(j,0);¥¥
{¥rm otherwise:}¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }-¥infty.
¥end{array}
¥right.¥¥
再帰： & {¥rm for¥ }j=1{¥rm ¥ to¥ }L,
d=1{¥rm ¥ to¥ }D (and d¥leq j), v=M {¥rm ¥ to¥ }1:¥¥
¥gamma_v(j,d) = &
¥left¥{
¥begin{array}{l}
{¥rm for¥ }s_v=E: ¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }-¥infty;¥¥
{¥rm for¥ }s_v=P and d<2:¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }-¥infty;¥¥
{¥rm for¥ }s_v=B, {¥cal C}_v = (y,z):¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }¥max{0¥leq k¥leq d}
[¥gamma_y(j-k,d-k) + ¥gamma_z(j,k)];¥¥
{¥rm otherwise:}¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ }¥max_{y ¥in {¥cal C}_v}
[¥gamma_y(j-¥Gamma_v^R,d-¥Gamma_v^L-¥Gamma_v^R)+¥log{t_v(y)}]¥¥
{¥rm ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ }+¥log¥hat{e}_v(x_i,x_j).
¥end{array}
¥right.¥¥
¥end{eqnarray}
図10.15 (a)開始位置を\(i\)，終了位置を\(j\)で添字付けした，\(L=10\)のデータベース配列に対する標準CYK DP行列の1レベル(モデルの1状態につき1個，\(M\)個の異なったレベルが3次元行列にある)．最大一致部分配列の長さを\(D=5\)と制限したときに，計算する必要がある行列の部分を白くして示した．データベース配列を(\(j\)を増加させながら)なめる検索アルゴリズムでの，行列要素の計算順序を矢印で示した．(b)同じCYK計算のための座標変換で，終了位置$j$と部分配列長\(d\) (\(d=j-i+1\))で添字付けされている．(b)の行列は\(L\times L\)ではなく\(D\times L\)であるから，この座標系を使った方が，必要なメモリ量が$L$に無関係ななめらかなCYKデータベース検索アルゴリズムを実装することが容易である．
図10.15 4個の異なった状態タイプに対するCYKデータベース検索アルゴリズムの再帰の様子．3次元DP行列のうち1個のレベルだけを示す．例えば左上では，
\(v\)とラベルされた要素の値\(\gamma_v(j,d)\)は\(v\)が結合しているそれぞれの状態\(y\)に相当する，\(y\)とラベルされた1個以上の配列要素\(\gamma_y(j,d-1)\)に依存している．これはDP行列の右に異なった形で示されている．もし\(v\)が1塩基を左側から生成すれば，\(j\)で終結する長さ\(d\)の部分文字列に対する状態\(v\)を根にもつ構文解析部分木が，\(y,j,d-1\)に対する部分木を付け加えることによって構築できる．\(R\)(右上)\(P\)(左下)に対する計算も同様である．\(v\)が二股枝分かれ(\(B\))状態であるときの計算は，最適な枝分かれ点を選ぶことに依存している．DP行列で\(y\)と\(z\)とラベルされた要素がそのような点の1個を示している．同様に結合されたその他の要素の集合をを灰色で示した．
図10.15に，このアルゴリズムの再帰のうちキーとなる段階を図示した．
いくつかの実装上の詳細は重要である．\(L\0行すべてを初期化する代わりに，初期化と繰り返しのステップを混ぜ合わせ，ある行\(j\)で最初\(\gamma_v(j,0)\)を初期化し，それから\(d=1,\ldots,D\)に対して\(\gamma_v(j,d)\)を計算した方が良い(長さ0の部分配列に対する初期化の計算は配列に独立であるから，\(\gamma_v(0,0)\)はすべての$v$に対して一度だけ計算すれば良く，その値はその後の\(\gamma_v(j,0)\)の初期化の時コピーすれば良い)．上のアルゴリズムと図10.15では，二股枝分かれ(\(B\))状態を除けば，行\(j\)のすべてのスコアは，行\(j\)と行\(j-1\)のスコアだけに依存することが明らかである．\(B\)状態はCMの構造の定義にしたがって2個の\(S\)状態に枝分かれするから，二股枝分かれ状態は，以前の\(D\)行の開始状態のスコアのみに依存する．このようにして，\(S\)状態のスコアの\(D+1\)行とその他の状態の2行のスコアだけがメモリ上に記憶する必要がある．
図10.16 CYKデータベース検索行列に重なり合わない2個の一致をを示す小さな例．2個の一致(位置5から10と13から19)は行列の左の箱として図示されている．
スコアの高い行列の要素\(\gamma_0(j,d)\)は黒く示されている．完全なアラインメントを再構築するためにメモリ上になければならない配列要素を灰色で示した．灰色の三角形は，要素\((j,d)\)の座標が，完全な配列をメモリ上に持たなくともそれぞれの一致の開始点$i$をどのように決定するかもしめしている．簡単のために，3次元行列の1レベルのみが示されている．
行\(j\)のスコア\(\gamma_0(j,d)\)は，位置$j$で終わるモデルに対する完全なアラインメント(すなわち，根状態\(v=0\)から始まる構文解析木)の対数オッズ比のスコアである．図10.16に示すように，一致の開始位置\(i\)は\(d\)から\(i=j-d+1\)と計算できる．これは，アラインメントの開始位置を得るのに，DP行列のトレースバックを必要としないということである．高スコアの\(\gamma_0(j,d)\)を見つけた後，CYK検索アルゴリズムが，スコアだけでなくこの高スコアを与える部分配列の開始位置$i$と終了位置\(j\)を，直接報告することができる．
実装では，一定の閾値を越えるすべての一致を報告するようにできる．しかし，高スコアアラインメントは，良いスコアで開始点や終了地点がわずかだけ違う，「影」を持っている．無制限の数の配列を一定のメモリでスキャンする単純なスコアの後処理アルゴリズムは以下のようである．行\(j\)を計算した後，行\(j\)における\(d\)に対する最高スコア\(\gamma_0(j,d)\)が決定される．もしこれが報告する閾値を越えていれば，これはリストに保存される．もしこれがリストにある以前の一致と重なりあっていれば，低い方スコアの一致が捨てられる．終了点\(j\)が現在の最小の開始点\(j-D\)より小さい全てのヒットは、重なり合わない一致として報告される．CYK検索アルゴリズムの時間複雑性は，\(M_a\)個の非枝分かれ状態と\(M_b\)個の枝分かれ状態のあるモデル，長さ\(L\)塩基のデータベース，最大一致長\(D\)塩基，に対して\(O(M_aLD+M_bLD^2)\)である．メモリの複雑性は\(O(M_aD+M_bD^2)\)である．計算時間はデータベースの大きさに対して線形に増え，必要なメモリ量はデータベースの大きさによらない．
[問題10.11]
内側アルゴリズムや外側アルゴリズムでは，行列での同様な座標変換を用いることができる．CYKと比べて，内側アルゴリズムは，部分配列に対するすべての可能な構造とアラインメントの確率を合計するという利点があるが，計算はCYKよりも複雑である．部分配列に対する長さが\(D\)以下の一致を探す内側アルゴリズムを与えよ．
構造アラインメント：トレースバック付CYKアルゴリズム
メモリの効率性の面から，行列要素の大部分は計算の途中で捨てられるから，上で述べた検索アルゴリズムではトレースバックができない．従って，アラインメントを回復することができない．スコアと，開始地点，終了地点のみが回復できる．メモリを消費して良ければ，トレースバックを行い，部分一致配列に対するSCFG最適構文解析木を回復するように，アルゴリズムを改良することは容易である．CMの構造はRNAの共通二次構造を反映するから，CM構文解析木は，モデルに対する最適アラインメントと最適二次構造予測の両方を表している．構文解析木の\(P\)状態(ペアワイズ)に対する割り付けは，予測された塩基対を示している．
第2章で述べた動的計画法のように，トレースバックポインタの行列を使うかスコア計算を再構築するかすれば，CYKのトレースバックを実装できる．一致をトレースバックすることを保証するには，トレースバックポインタかスコアの\(D+1\)個のすべての行を，メモリ上に置かねばならない．もし上のようなオーバーラップ処理アルゴリズムを使えば，ある一致が後で調べる一致とオーバーラップしていないことを確かめるためにさらに\(D\)行を必要とするから，完全な\(2D+1\)行がメモリに保持しなければならない．もしスコアとアラインメントの両方が欲しければ，CYK局所検索アルゴリズムとトレースバック付きCYK大域アラインメントアルゴリズムの両方を実装することが効率的であり，単に以下の2段階を行えば良い．最初に一致部分配列を探す局所検索を行う．次にトレースバック付きの大域アラインメント用のモデルにその部分配列を最適にアラインメントする．
トレースバックは位置\(j\)で終わる長さ\(d\)の高スコア部分配列\(\gamma_0(j,d)\0から始まり，後向きにたどっていく．配列に基づく局所ではなく大域のアラインメントに関しては，トレースバックは\(\gamma_0(L,L)\)から始まる．以前の章で述べたHMMの動的計画法のトレースバックによる回復やより単純なRNAモデルからのSCFGのトレースバックと，ここに述べた処理は基本的に同一であり，詳細はここでは述べない．
[問題10.12]
CYKアルゴリズムを改良し，最適構文解析木を回復するのを助けるためにそれぞれのセルのトレースバックの情報を保持するようにせよ．二股枝分かれ状態からトレースバックするために，保持する必要のある最低限の情報は何か？その他の任意の状態からトレースバックするために，保持する必要のある最低限の情報は何か？
CMを用いた「自動化」比較配列解析
アラインメントされていない
基本的な考え方は次の2段階を繰り返すことである．
(a)与えられた現段階のアラインメントから，最適な(あるいはほとんど最適な)CM構造を構築する．
(b)現段階のCMに対して最適なマルチプルアラインメントを構築する．
(a)に関しては，与えられたRNAマルチプルアラインメントに対して，共通構造を見つけるにはいくつかの方法が可能である．比較配列解析法から考え出されたヒューリスティック法が[Eddy, S. R. & Durbin, R.1994]によって使われた．アラインメントされた列のすべてのペアに対して，相互情報量\(M_{ij}\)が最初に計算される．次に，DP折れ畳みアルゴリズム(基本的にはNussinovアルゴリズム)を使って，\(M_{ij}\)の合計が最大化される構造木が見つけられる．このアルゴリズムが最大の合計$S_{ij}$を見つけるためのDP行列書き込み段階は以下のように書ける．
\[
\begin{array}{ll}
S_{i,j} = \max
\left\{ \begin{array}{ll}
S_{i+1,j} & 列 i が非ペア\\
S_{i,j-1} & 列 j が非ペア\\
S_{i+1,j-1}+M_{i,j} & 列 i, j がペア\\
\max_{i<k<j}S_{i,k}+S_{k+1,j} & 二股枝分かれ\\
\end{array}
\right.
\end{array}
\]
この行列のトレースバックは共通構造木を産み出し，CMに拡張される．この手法の利点は，相互情報量を用いることで，比較解析で最も支持される塩基対だけがペアにされることである．この方法のひとつの欠点はペアを作りすぎることである．なぜなら\(M_{ij} \geq 0\)であり，見せかけのペアを加えることによっても\(S_{ij}\)はわずかに増加するからである．2番目の欠点は共変関係だけを探すために，配列上の変異がほとんどない高度に保存された構造やドメインを推定し損なう可能性があることである．
次に，アラインメントされた配列に共通構造を当てはめ，それぞれの配列の構文解析木を見つけ，構文解析木から遷移度数，出力度数を数え，通常の方法で度数を確率変数に変換する．
注釈なしのRNAマルチプルアラインメントを使って，最適な(最大事後確率)CM構造とそのパラメータを同時に抽出する，厳密なCM構築アルゴリズムも存在する[Eddy, S. R., unpublished]．このアルゴリズムは，第5章で述べたプロファイルHMMに対するMAP構築アルゴリズムの拡張である．
現段階のモデルに対して最適なマルチプルアラインメントを構築するという(b)に関しては，CYKアラインメントアルゴリズムをそれぞれの配列に適用する．HMMの場合と同様に，配列のマルチプルアラインメントは，共通モデルに対する個別のアラインメントの集合から導かれる．
Eステップで内側・外側アルゴリズムの代わりにCYKアルゴリズムを用い，Mステップで(パラメータだけでなく，モデルの構造とパラメータを同時に再推定する)モデル構築アルゴリズムを用いることを除けば，これはEMアルゴリズムと似ている．
アルゴリズムは，人間の比較解析者が行うのとほとんど同じやり方で学習を行う．いくつかの場合には，同じ答えを与える．
当初アラインメントされていない配列のアラインメントを適当に決めることから(ランダムもしくは配列のみに基づくアラインメント)始めて，構造を推定し，推定された構造からアラインメントをやり直し，構造を最推定するという手続きを，一貫性のある解が得られるまで繰り返す．これは，[Eddy, S. R. & Durbin, R.1994]で使われた100本の転写RNA配列のような，理想的なデータセットに対してはうまく行く．
しかし，そのように理想的なデータセットはあまりない．アルゴリズムがうまく行くためには，数多くの小さなRNAが必要である．塩基対の大部分を明らかにするような十分な一次配列上のばらつきも必要である．共通二次構造は高度に保存されていなければならない．配列は局所的ではなく，全体としてアラインメントできなければならない．アルゴリズムは局所最適に陥りやすく，しかもプロファイルHMMに用いたようなシミュレーテッドアニーリング(第5章)を使うには計算時間がかかりすぎる．これらは単に克服すべき技術的な問題と言うこともできるが，自動比較解析という点では，SCFGに基づく解析は未だに実用と言うより理論的な興味の対象である．CM法の最も実用的な応用は，次の節で述べるように，相同的なRNA構造のデータベース検索である．
CMアルゴリズムの実用的な応用の例
実用的な配列解析問題にこの理論がどのように適用されるかをもっと良く理解するために，CM法の実際への応用を見てみることは意味がある．
大部分のゲノムでは，最も大きな遺伝子ファミリーはタンパク質遺伝子ファミリーではなく，構造RNA遺伝子のファミリー，転移RNA(tRNA)である．例えば，酵母Saccharomyces cerevisiae(出芽酵母)では，274個のtRNA遺伝子があり，ヒトゲノムにはおおよそ1500個のtRNA遺伝子がある．ゲノム配列中のtRNA遺伝子を見つける数多くのプログラムが開発されてきた．これらのtRNA発見プログラムは通常，大規模なゲノム注釈付けプロジェクトの一部である．このような手動で注意深く最適化されたプログラムの偽の正解の割合は，100万塩基のDNAに対して\(0.2\)-\(0.4\)個の間違った予測を与える程度である[Fichant, G. A.& Burks, C. 1991;Pavesi, A. et al. 1994]．この偽の正解の割合は酵母のような小さなゲノム(14 Mb)に対しては受け入れ可能だが，3000 Mbのヒトゲノムでは，このようなプログラムは1000個前後の偽の正解を産み出すから，予測されたtRNAの多くの部分が偽の正解ということになる．
転移RNAはCM法をテストするのに理想的な候補である．これらは短く(通常$75$-$95$塩基)，一次配列の同一性は低く，「クローバ葉」状の良く保存された二次構造のを持ち，何千ものtRNA配列を確率モデルの学習のために用いることができる．オルガネラやウイルスのtRNAを含む，様々な生物種の1415個のtRNAからなるマルチプルアラインメントが，tRNAのCMを構築するための出発点であった[Steinberg, S.,Misch, A.& Sprinzl, M. 1993]．
多くの真核生物のtRNA遺伝子は短いイントロンを持つことに「気が付く」ようにモデルを学習するために，38本の配列は短いイントロンをアンチコドンループにもっている．様々な位置にあるより長い多くのイントロン(通常はcatalytic グループIとグループIIのイントロン)はアラインメントから除かれている．なぜなら，これらはこのアルゴリズムで探すには長すぎるからである．モデルは，共通二次構造の注釈を用いて，このアラインメントから直接に作られる．モデル構築のための計算時間はほとんど無視できる．得られたモデルは285状態を持つ．状態は共通モデルで72ノード，3個の二股枝分かれ状態，28個のペアワイズ状態，33個の左側生成状態，1個の右側生成状態，7個の開始状態にグループ分けされる．28個のペアワイズノードはtRNAのクローバ葉の4個のヘリックスと可変アームヘリックスの28個の共通塩基対に相当する．
CYKデータベース検索アルゴリズムの1個の実装(COVEプログラムパッケージのCOVELS)は，最大一致の長さを\(D=150\)の制限したとき，SGI Indigo2 R4400ワークステーションを用いて1秒当り20塩基のDNAを検索した[Eddy, S. R. & Durbin, R. 1994]．探索行列で必要なメモリ空間は500Kである．従って，tRNA検索では，メモリ所要量よりもCPU時間が問題である．MasPar並列計算機上に実装された並列アルゴリズムは，毎秒2000塩基まで検索を高速化したが，特別なハードウェアを用いることは実用的ではない．
そこで，2個の異なった既存のtRNA発見プログラムと高速前処理フィルターを用いた，TRNASCAN-SEと呼ばれるハイブリッドプログラムが書かれた[Lowe, T. M. & Eddy, S. R. 1997]．どちらか一方か両方のプログラムで候補とされたものが第二段階でtRNAのCMに対して一致させられ，統計的に有意なもの(対数確率スコアが\(\geq 20\)ビット)が報告される．表10.02にこのハイブリッドプログラムと，他のtRNA発見プログラム，単なるCM，と比べた性能をまとめた．
¥begin{table}
¥centering
{¥small
¥begin{tabular}{lrrr}
¥hline
&速度 &真の正解 &偽の正解¥¥
プログラム &(bp/sec) &(¥%) &(per Mb)¥¥
¥hline
TRNASCAN 1.3 [¥AINDEX{Fichant, G. A.}{Fichant}
¥& ¥AINDEX{Burks, C.}{Burks} 1991]
&400 &95.1 &$0.37$¥¥
POL3SCAN [¥AINDEX{Pavesi, A.}{Pavesi} {¥em et al.} 1994]
&373000 &88.8 &$0.23$¥¥
CM alone [¥AINDEX{Eddy, S. R.}{Eddy} ¥& ¥AINDEX{Durbin, R.}{Durbin} 1994]
&20 &99.8 &$<0.002$¥¥
TRASCAN-SE [¥AINDEX{Lowe, T. M.}{Lowe} ¥& ¥AINDEX{Eddy, S. R.}{Eddy} 1997]
&30000 &99.5 &$<0.00007$¥¥
¥hline
¥end{tabular}
}
¥caption{tRNA遺伝子検出手法の比較}¥label{tab:10.02}
¥end{table}
(ほとんど)自動的なCMに基づくアプローチであるTRNASCAN-SEは，手動で改良された他の2個の方法よりも，精度は大きく感度は少し改良されている．3000Mbのヒトゲノム中の偽の正解tRNAの期待値は1000以上から1以下に減少した．
精度と感度に加えて，CMアプローチのもう一つの利点はその一般性である．RNAのどんなアラインメントからでもモデルを作ることができる．例えば，TRNASCAN-SEでは，セレノシステインtRNA(大部分のtRNAとは様々な点で異なっている)に対する特別なモデルを簡単に作ることができて，プログラムがこれらのtRNA遺伝子を検出するようにすることができる．
RNA解析でのCM法の主要な欠点は，大きなメモリとCPU時間を必要とすることである．TRNASCAN-SEでは，他の既存のプログラムを前処理のフィルターとして用いることにより，時間の問題は克服できる．しかし，このアプローチは一般的ではない．なぜなら，他の多くのRNAファミリーではそのようなプログラムが利用できないからである．より大きなRNAにCM法を用いるには，より良いアルゴリズムとより良い計算機の片方か両方が必要である
- 開始非終端記号に枝分かれする必要はないが，その方がこれ以降のアルゴリズムを単純にできる．こうするひとつの理由は，枝分かれ開始の結合を切り放し，構造のそれぞれの枝を独立したRNAドメインの独立したSCFGとして扱うことができることである． ↩
- もしこのSCFGをRNA解析に用いるのなら，内側・外側およびCYKアルゴリズムをアラインメントに用いるのは無駄である．ギャップがなければ，1個のアラインメントだけが可能である．ギャップなしのRNAモチーフは$O(L)$時間で確率モデルのスコアが計算できる． ↩
- もう一つの可能な配置は，\(S\)と\(IR\)を左の子供に割り当て，\(S\)を右の子供に割り当てることである．選択の余地があるときには，プロファイルHMMと同様な左側生成を行うことを原則とするので，この配置を選んだ． ↩
- 実際には，内側アルゴリズムの\(s_v=S\)である\(\alpha_v\)だけが必要である．このことはメモリが限られているときは有用である． ↩
- \(B\)状態の遷移確率は定義により1であるから，最推定する必要はない． ↩
- 完全な確率モデルであるためには，ランダムモデルにも継続時間長分布を指定しなければならないが，この項は通常無視される． ↩