第1回 2019/6/10

• 講師：浅井 潔
• 資料：バイオインフォマティクス 確率モデルによる遺伝子配列解析 第９章

Lecture Schedule

• 講義HP

講義概要

今回の講義は、上記の資料（浅井先生が翻訳に関わられている本）の文章をそのまま追うだけであった。

変形文法

ゲノム配列解析論Ⅰでは、「配列のアラインメント」や「隠れマルコフモデル」など、生物配列を一次元の記号列として扱っていた。この仮定は、計算するには扱いやすいが、生物配列の性質を十分に表現できていない。

そこで、タンパク質や塩基の三次元的な折れ畳み（一次元の配列上では隣り合っていない残基の広範囲な物理的相互作用）を考慮に入れた確率モデルを作ることを考える。なお、その例としてRNAの二次構造解析を見る。

生物配列の性質を取り入れると言っても、結局は配列を「ある文法にのっとった文字列」として見る（RNAの二次構造解析では、確率自由文法(SCFG)を利用する）ため、形式言語の基本を理解する必要はある。

そのため、今日は

1. 確率をもたない形の変形文法を概観
2. 遠距離の相関と制約のある配列の完全な確率的モデル化のための形式システムとしての確率文法の導入
3. 確率文脈自由文法に対する一般化したアラインメントアルゴリズムの提示

といった流れで確率文脈自由文法に基づく確率モデルをRNA二次構造問題に適用するための土台を準備する。

1. 変形文法

変形文法は、有限個の記号と、\(\alpha\rightarrow\beta\) の形の有限個の書き換え規則(rewriting rule)から成る。\(\alpha\) と \(\beta\) はともに記号の列であり、書き換え規則は生成規則(production rule)とも呼ばれる。

記号には、以下の2種類があり、左辺の \(\alpha\) は少なくとも1個の非終端記号を含み、一般に右辺の非終端記号や終端記号の新しい列に変換される。

• 抽象的な非終端記号(nonterminal symbol)
• 観測文字列に実際に現れる終端記号(terminal symbol)

例えば、2文字の終端記号のアルファベット \(\{a, b\}\) とひとつの非終端記号 \(S\)、処理を終わらせるための特殊な空終端記号 \(\epsilon\) を用いて \(a\) と \(b\) からなる任意の文字列を生成する変形文法は以下のように書かれる。

部分文字列を選び、許された書き換え規則のひとつに従ってそれを書き換えるプロセス（非終端記号がなくなるまで続く）の結果生じる文字列の系列を文法からの導出(derivation)と呼ぶ。

配列解析の問題に変形文法を用いる時は、特定の配列がその文法に適合する（その文法から生成された）かどうかが問題となるため、文字列の導出が存在するかを逆向きにたどることになる。

与えられた配列に対して有効な導出を見つけることを構文解析(parsing)と呼んでおり、その意味で「導出は配列の構文解析」とも呼ばれる。

チョムスキー階層

• \(W\) を任意の非終端記号
• \(a\) を任意の終端記号
• \(\alpha\) と \(\gamma\) をヌル記号を含む非終端記号や終端記号の任意の系列
• \(\beta\) をヌル記号を含まない非終端記号や終端記号の任意の系列

と定義すると、それぞれの階層は以下のように説明できる。

オートマトン

上記の図形のように、それぞれの文法はオートマトンと呼ばれる抽象的な計算機構に対応している。オートマトンは配列を受理あるいは棄却する構文解析器として記述される。

2. 正規文法

正規文法(regular grammar)の全ての生成規則は \(W\rightarrow aW\) または \(W\rightarrow a\) の形をしている。しばしば \(W\rightarrow\epsilon\) という生成規則を追加することがあるが、正規文法は常に \(\epsilon\) を吸収するように拡張できることが知られている（例：正規文法に近い \(S\rightarrow aS|bS|\epsilon\) は \(S\rightarrow aS|bS|a|b\) と同値）ので、アルゴリズムに深刻な問題とはならない。

正規文法は配列を左（右）から右（左）に発生するため、終端記号の間の遠距離の相関を効果的に記述一時配列モデルであると言うことができる。

有限状態オートマトン

正規文法に対応する構文解析オートマトンは有限状態オートマトン(finite state automaton)である。

有限状態オートマトンはいくつかの状態(state)からなるモデルで、状態は状態遷移(state transition)によって互いに結合している。状態と状態遷移は同値な正規文法の非終端記号と生成規則に対応している。

ムーア機械とミーリー機械

有限状態オートマトンには以下の二つのタイプの機械があるが、相互変換可能である。

• ムーア機械(Moore machine)：状態で記号を受理するもの
• ミーリー機械(Mealy machine)：遷移時に記号を受理するもの

決定性オートマトンと非決定性オートマトン

有限状態オートマトンにはさらに以下の二つの性質があり、任意の非決定性有限状態オートマトンは決定性有限オートマトンに変換可能であることが知られている。

• 決定性：任意の状態で任意の入力記号を受理する1個以上の遷移はない。
• 非決定性：ある状態で、ある記号列を受理する複数の遷移がある。

決定性有限状態オートマトンによる構文解析は極めて効率的であり、高速なBLASTデータベース検索プログラムに生かされている。

一方の非決定性有限状態オートマトンは配列を棄却する前に全ての代替経路を調べないといけないが、効率的に行うことができ、grepやsed,awk,vi などのUNIXのユーティリティにおけるテキストのパターン称号プログラムは高度に効率的な非決定性有限状態オートマトンを用いている。（※UNIXの「正規表現」は正規文法と同値である。）

正規文法でできないこと

正規文法が記述できない言語 \(L\) の古典的な2個の例題：

1. \(L\) が \(aa,bb,abba,abaaba\) といった、逆から読んでも同一な文字列全てを含む（回文言語）
2. \(L\) が \(aa,bb,abab,aabaab\) といった、同一の2個の部分からなる文字列全てを含む（複写言語）

※ なお、RNAの二次構造予測の場合、回文構造を読むことは極めて大事である。

3. 文脈自由文法

回文構造は、チョムスキー階層の次の段階、文脈自由文法(CFG)によって扱われる。RNAの二次構造は一種の回文構造になっており、「入れ子になった位置（ステムループの架橋構造）の塩基対の強い相互作用が維持される限り、配列自体はどうなっても良い」と言う問題を提示している。

文脈自由文法の生成規則の左辺は正規文法と同様に1個の非終端記号でなければならないが、右辺は終端記号と非終端記号の任意の組み合わせであって構わない。ゆえに、右辺において相関のある塩基対を生成することができる。

回文を生成するCFGの例は以下のようなものである。

正規文法は文字を左から右に生成するが、文脈自由文法は文字列を外側から内側に生成するため、入れ子になった相関だけが捕らえられる。

※ 複写言語は入れ子の制限を破るため、文脈自由言語ではない。

プッシュダウンオートマトンは構文解析木を一次元でやっている（プッシュとポップを繰り返せばできる。）

構文解析木

文脈自由文法の配列に対するアラインメント（＝構文解析）には構文解析木(parse tree)と呼ばれるエレガントな表現がある。任意の部分木(subtree)は配列の連続した断片を導出するという重要な性質を持っている。

プッシュダウン・オートマトン

CFGのための構文解析オートマトンはプッシュダウン・オートマトン(push down automaton)と呼ばれるが、これを用いることで、構文解析木が行なっていることを一次元で行うことができる。

プッシュダウン・オートマトンはスタックをメモリとして保持しており、プッシュとポップを使うことで構文解析する。

アルゴリズムは以下の通り。なお、開始非終端記号をプッシュすることで初期化され、入力記号がなくなるまで以下のステップを逐次的に繰り返す。入力記号がなくなった時にスタックが空なら、その配列の構文解析は成功したことになる

非決定性文脈自由文法に対しては、オートマトンの全ての有効な動作が網羅的に試されるため、プッシュダウン・オートマトンは効率の良い認識器ではない。

文脈自由文法に対するより洗練された多項式時間(polynomial time)のCocke-Younger-Kasami(CYK)構文解析アルゴリズムが存在する。

4. 文脈依存文法

複写言語は文脈自由言語ではなく、文脈依存文法を必要とする。

文脈依存文法は、複写言語の交差して相互作用している記号の組を直接には生成しない。その代わり、非終端記号 \(W\) が交差しない相互作用の組としてそれらを生成し、それから文法が局所的な文脈を調べることによって非終端記号を並び替える。

並び替えの規則は 「 \(\hat{}\) のついた非終端記号を \(\hat{}\) のつかない非終端記号と入れ替えること」である。なお、任意の生成規則はその左辺が現れた時いつでも適用可能であるので、非終端記号が正しく並び替えられるまでは終端記号を生成し始めないように注意深く文法を作らないといけない。

例

\(cc,acca,abaccaba, bbabccbbab\) のような、\(a\) と \(b\) からなる文字列の2個のコピーの間に2個の \(c\) が挟まった全ての文字列からなる複写言語を考える。

• 初期化
  • \(S\rightarrow CW\)
• 非終端記号の生成
  • \(W\rightarrow A\hat{A}W|B\hat{B}W|C\)
• 非終端記号の並べ替え
  • \(\hat{A}B\rightarrow B\hat{A}\)
  • \(\hat{A}A\rightarrow A\hat{A}\)
  • \(\hat{B}A\rightarrow A\hat{B}\)
  • \(\hat{B}B\rightarrow B\hat{B}\)
• 終端記号の生成
  • \(CA\rightarrow aC\)
  • \(CB\rightarrow bC\)
  • \(\hat{A}C\rightarrow Ca\)
  • \(\hat{B}C\rightarrow Cb\)
• 終了
  • \(CC\rightarrow cc\)

上記の文法から \(aabccaab\) を生成する導出は以下のようになる。

文脈依存文法のための構文解析オートマトンは線形有界オートマトン(linear bounded automaton)である。

線形有界オートマトンは、「導出が開始非終端記号にたどり着く」か、「たどり着かずに全ての可能な導出を調べ尽くす」かするまで、観測文字列に対する全ての可能な導出を規則正しく逆向きにたどる仕組みである。

任意の文脈依存文法は単調文法(monotonic grammar)であり、「\(u\rightarrow v\) という生成規則を持つ時 \(|u|\leq |v|\) かつ \(u\) は少なくとも1つの非終端記号を含まなければならない」という制約を満たしている。

したがって、線形有界オートマトンがトレースバックするときに調べるべき可能な導出は必ず有限になる。しかしながら、可能な導出の数は指数的に大きい。

文脈依存文法に対する一般的な多項式時間のアルゴリズムは知られていないため、文脈依存文法の実用的な応用を考える場合にはこれが深刻な問題となる。

NP問題と「手に負えなさ」

「解を見つける多項式時間のアルゴリズムは知られていないが、解が正しいかどうかは多項式時間で検査できる」問題は、NP問題(nondeterministic polynomial problem)と呼ばれている。

文脈依存文法の構文解析を含む多くのアルゴリズムがNP問題である。NP問題は「手に負えない(intractable)」とも呼ばれるが、多くのモデル化の問題では「一般的にはNPだが、扱いやすい特殊な場合が存在する」ことがあり、それを理解することが重要である。

句構造文法とチューリング機械

• 句構造文法(phase structure grammar)は変形文法で、生成規則の左辺と右辺が記号の任意の組み合わせであって構わないものである。
• 同値な構文解析オートマトンはチューリング機械(Turing machine)である。
• ある文字列に句構造文法からの有効な導出があるかどうかを有限時間で決定する保証された一般的なアルゴリズムは存在しない。

5. 確率文法

パターンの規則に対する例外はどのような位置にも起こり得る。ここで、例外のためにパターンを増やしていくと、パターンの特異性が下がっていき、結果として関係のないランダムな配列と一致してしまうほどパターンは情報のないものとなってしまう。

そこで、例外を許すときに「全ての可能性を同等に扱う」のではなく、「強い共通部分より例外に対して少ないスコアを与える」ことでこの問題を解決する。

チョムスキー階層のどの文法も、確率的な形をとることによって配列の確率的モデル化の基礎として使える。確率的文法モデル \(\theta\) は異なった文字列 \(x\) を確率 \(P(x|\theta)\) で生成する。

確率文脈自由文法か句構造文法か

文脈依存文法や句構造文法に対する一般的な多項式時間のアルゴリズムは存在せず、実用的な応用があるとは思えないので今後議論しないが、確率的にこれらを扱う際には注意深く定式化しなければならないことだけ留意しておく必要がある。

例えば、以下の文脈依存文法を考える。

このとき、この文法から生成される言語は \(\{aa,ab,ba,bb\}\) であり、その確率は \(\{p_1p_4,p_1p_5,p_2p_4,(p_2p_3 + p_2p_5)\}\) である。単に \(S\) と \(W\) に対する生成の確率の和が \(1\) 、つまり \(p_1+p_2=1\) かつ \(p_3+p_4+p_5=1\) であるだけでは、\(p_1=0\) または \(p_3=0\) である特殊な場合を除いて言語全体に対して確率分布を与えたことにはならない。

隠れマルコフモデルは確率正規文法である

隠れマルコフモデル(HMM)は確率正規文法と同様である。

唯一の違いは2種類のモデルは伝統的に異なった表現能力を持つことで、

• HMMは通常ムーア機械として、遷移に関わりなく状態から記号を出力するように記述される。
• 確率正規文法の生成はミーリー機械として、新たな非終端記号への遷移から記号を出力するように記述される。

6. 配列モデルのための確率文脈文法

文脈自由文法(CFG)は書き換え規則の右辺で無制限の記号文字列を持っている。一般的なCFG構文解析アルゴリズムを表現するのに、制限された「標準形」を書き換え規則に採用するのは有用である。

そういった標準形で最も一般的なものとしてチョムスキー標準形が挙げられる。チョムスキー標準形は、次のような生成規則のみからなる文法である。なお、任意のCFGがチョムスキー標準形に変換できることが知られている。

さて、これで配列のモデルとしての確率文脈自由文法を書き下した後に、次の3つの問題を解かなければならない。なお、HMM（確率正規文法）の場合と比較させた。

内側アルゴリズム

内側アルゴリズムは、SCFGが与えられた時の配列の確率（スコア）を計算する。

まず、いくつかの記法を定義する。

• \(M\) 個の非終端記号 \(W = W_1,\ldots,W_M\) をもつ
• 開始非終端記号は \(W_1\) である
• \(v,y,z\) を非終端記号 \(W_v,W_y,W_z\) に対する添字とする。
• 生成規則は \(W_v\rightarrow W_yW_z,\ W_v\rightarrow a\)（\(a\) は終端記号のアルファベットの中の可能な記号の1個）の形をしている
• 生成に対する遷移(transition)の確率パラメータを \(t_v(y,z)\) と呼ぶ
• 生成に対する出力(emission)の確率パラメータを \(e_v(a)\) と呼ぶ
• 配列 \(x\) は \(L\) 個の記号を持ち、\(x_1,\ldots,x_L\) と書かれる
• \(i,j,k\) を配列 \(x\) 中の記号 \(x_i,x_j,x_k\) に対する添字とする

内側アルゴリズムは、「部分文字列 \(x_i,\ldots,x_j\) に対応する非終端記号 \(W_v\) に根を持つ構文解析部分木の確率 \(\alpha(i,j,v)\) 」を、すべての \(i,j,v\) に対して計算する。

計算は \(L\times L\times M\) の3次元のDP行列を必要とし、長さ \(1(i=j)\) の部分配列から始めてより長い部分文字列へと、開始非終端記号に根を持つ完全な構文解析木の確率が求まるまで続けられる。

形式的には内側アルゴリズムは以下のように書かれる。

• 初期化：
  $$\alpha(i,i,v) = e_v(x_i)\quad (\forall i, v)$$
• 再帰：
  $$\alpha(i,j,v) = \sum_{y=1}^M\sum_{z=1}^M\sum_{k=i}^{j-1}\alpha(i,k,y)\alpha(k+1,j,z)t_v(y,z)$$
• 終了条件：
  $$P(x|\theta) = \alpha(1,L,1)$$

なお、この時再帰式で以下のように \(i,j\) をともに左から走査するとDPがうまく回らない（求め終わっていない部分問題を利用しようとしてしまう）

for (int i=1; i<=L-1; i++){
  for (int j=i+1; j<=L; j++){
    for (int v=1; v<=M; v++){
      再帰式
    }
  }
}

したがって、\(j-i\) を次第に大きくするように計算していく必要がある。

内側アルゴリズムはこのようにして配列のSCFGによる確率（スコア）を計算する。

• \(\alpha\) に対する3個の添字から明らかなように、記憶領域複雑度は \(O(L^2M)\) である。
• 配列の3個の位置の添字 \(i,j,k\) と文法の非終端記号の3個の添字 \(v,y,z\) に関する再帰繰り返しの回数から明らかなように、時間複雑度は \(O(L^3M^3)\) である。

外側アルゴリズム

参考：The Inside-Outside Algorithm

外側アルゴリズムは、「開始非終端記号 \(W_1\)（図では\(S\)）を根に持ち、\(x_1,\ldots,x_{i-1},W_v,x_{j+1},\ldots,x_L\) という記号列に構文解析できる確率 \(\beta(i,j,v)\) 」を、すべての \(i,j,v\) に対して計算するもので、以下の特徴を持つ。

• 計算は \(L\times L\times M\) の3次元のDP行列を必要とする。
• 外側 \(\beta(i,j,v)\) 確率を計算するのには内側アルゴリズムの計算結果 \(\alpha(i,j,v)\) を必要とする。
• \(W_v\) に \(x_1,\ldots,x_L\) が含まれるものから始め、内側に向かって進む。

形式的には外側アルゴリズムは以下のように書かれる。

• 初期化：
  • \(\beta(1,L,1) = 1\)
  • \(\beta(1,L,v) = 0\ (1<v\leq M)\)
• 再帰：
  $$\beta(i,j,v) = \sum_{y,z}\sum_{k=1}^{i-1}\alpha(k,i-1,z)\beta(k,j,y)t_y(z,v) + \sum_{y,z}\sum_{k=j+1}^{L}\alpha(j+1,k,z)\beta(i,k,y)t_y(v,z)$$
• 終了条件：
  $$P(x|\theta) = \sum_{v=1}^{M}\beta(i,i,v)e_v(x_i) \text{ for any } i$$

この再帰式を理解するには、以下の2パターンを考えれば良い。

• \(N_y\rightarrow N_zN_v\) の時（上図） \(x_1,\ldots,x_{k-1},W_y,x_{j+1},\ldots,x_L\) という形に構文解析できる確率は \(\beta(k,j,y)\) で求められているから、そこに「 \(N_y\rightarrow N_zN_v\) と遷移する確率 \(t_y(z,v)\) 」と 「部分文字列 \(x_k,\ldots,x_{i-1}\) に対応する非終端記号 \(W_z\) に根を持つ構文解析部分木の確率 \(\alpha(k,i-1,z)\) 」をかければこの場合の \(\beta(i,j,v)\) は求めることができる。あとは、\(k\) の位置を \(1\) から \(i-1\) まで動かして総和を求めれば良い。
• \(N_y\rightarrow N_vN_z\) の時（下図） \(x_1,\ldots,x_{i-1},W_y,x_{k+1},\ldots,x_L\) という形に構文解析できる確率は \(\beta(i,k,y)\) で求められているから、そこに「 \(N_y\rightarrow N_vN_z\) と遷移する確率 \(t_y(v,z)\) 」と 「部分文字列 \(x_{j+1},\ldots,x_k\) に対応する非終端記号 \(W_z\) に根を持つ構文解析部分木の確率 \(\alpha(j+1,k,z)\) 」をかければこの場合の \(\beta(i,j,v)\) は求めることができる。あとは、\(k\) の位置を \(j+1\) から \(L\) まで動かして総和を求めれば良い。

期待値最大化によるパラメータ再推定

HMMのEMによる学習で前向き変数、後ろ向き変数が使われたのと同様に、SCFGの確率パラメータの期待値最大化における再推定に内側変数 \(\alpha\) と外側変数 \(\beta\) を使うことができる。

これらを使うと、導出の中で状態 \(v\) が使われる回数の期待値は以下の式で表される。

これは、さらに以下のように拡張して、\(W_v\) が占められ生成規則 \(W_v\rightarrow W_yW_z\) が使われる回数の期待値を見つけることができる。

生成規則 \(W_v\rightarrow W_yW_z\) の確率のためのEM再推定の式は、

同様の式は他の生成規則 \(W_v\rightarrow a\) についても成り立つから、

これらの1個の観測配列 \(x\) からの再推定は、独立した複数の観測配列の場合に容易に拡張できる。期待される度数は単に全ての配列について和を取れば良い。

CYKのトレースバック

残った問題は、「配列に対する最適な構文解析木を求めること」である。これは、Coke-Younger-Kasam(CYK)アルゴリズム、内側アルゴリズムで確率の合計の代わりに最大の確率を与える操作を使うものによって解かれる。

このアルゴリズムは、最終的に \(P(x,\hat{\pi}|\theta)\)（\(\hat{\pi}\) は最も確率の高い構文解析木）になる変数 \(\gamma(i,j,v)\) を計算する。なお、この時3次元のDP行列でトレースバックして最適のアラインメントを回復するために必要な \((y,z,k)\) という3個の数字の組の後戻り変数 \(\tau(i,j,v)\) も保持しなければならない。

式で書くと、行列を埋めるアルゴリズムは以下のようになる。

CYKアルゴリズム

• 初期化：for \(i=1\) to \(L\), \(v=1\) to \(M\):
  $$ \begin{aligned} \gamma(i,i,v) &= \log e_v(x_i)\\ \tau(i,j,v) &= (0,0,0) \end{aligned} $$
• 再帰：for \(i=1\) to \(L-1\), \(j=i+1\) to \(L\), \(v=1\) to \(M\):
  $$ \gamma(i,j,v) = \max_{y,z}\max_{k=1,\ldots,j-1}\{\gamma(i,k,y),\gamma(k+1,j,z)+\log t_v(y,z)\}\\ \tau(i,j,v) = \arg\max_{(y,z,l),k=i,\ldots,j-1}\{\gamma(i,k,y),\gamma(k+1,j,z)+\log t_v(y,z)\} $$
• 終了：
  $$\log P(x,\hat{\pi}|\theta) = \gamma(1,L,1)$$

このあと、最善のアラインメントを回復するために、プッシュダウン・スタックから三つ組 \((i,j,v)\) をプッシュとポップしながら逆戻りを行う。

CYK逆戻りアルゴリズム

(1,L,1)をスタックにプッシュする。
while()
(i,j,v)をポップ
if τ(i,j,v)==(0,0,0): # (i==j)
  x_iをvの子供として割り当てる。
else:
  y,zをvの子供として構文解析木に付け加える。
  (k+1,j,z)をプッシュ
  (i,k,y)をプッシュ

SCFGアルゴリズムのまとめ

ここまでチョムスキー標準形の場合で説明したが、内側・外側アルゴリズムが適用可能なのはチョムスキー標準形のSCFGだけではない。（一般性と記法の便利さゆえ）

実際、RNAをモデル化する別のSCFGが存在する。