生物情報ソフトウェア特論 （８） RNA二次構造予測

Presentation on theme: "生物情報ソフトウェア特論 （８） RNA二次構造予測"— Presentation transcript:

1 生物情報ソフトウェア特論 （８） RNA二次構造予測
阿久津　達也 京都大学　化学研究所 バイオインフォマティクスセンター

2 講義予定 第１回: 文字列マッチング・データ構造 第２回： たたみ込みとハッシュに基づくマッチング 第３回： 近似文字列マッチング
第１回:　文字列マッチング・データ構造 第２回：　たたみ込みとハッシュに基づくマッチング 第３回：　近似文字列マッチング 第４回：　配列解析 第５回：　木構造の比較：順序木 第６回：　木構造の比較：無順序木 第７回：　文法圧縮 第８回：　RNA二次構造予測 第９回：　タンパク質立体構造の予測と比較 第１０回：　固定パラメータアルゴリズムと部分k木 第１１回：　グラフの比較と列挙 第１２回：　ニューラルネットワークの離散モデル

3 RNA二次構造予測

4 RNA二次構造予測 RNA二次構造予測（基本版） 塩基対 入力： RNA配列 a=a[1]…a[n] 出力： 以下を満たし、スコア 二次構造
U G C 塩基対 二次構造 二次構造でない RNA二次構造予測（基本版） 入力： RNA配列 a=a[1]…a[n] 出力： 以下を満たし、スコア 　Σ(i,j)∈M μ(a[i],a[j]) が最小となる塩基対 の集合 M={(i,j)|1≤i+1<j≤n,{a[i],a[j]}∈B} (a[i],a[j]) ,(a[h],a[k]) ∈M となる i ≤h ≤j ≤k が存在しない 塩基対： B={{a,u},{g,c}｝ スコア関数（最も単純なもの） μ(a[i],a[j])= -1 if {a[i],a[j]} ∈B μ(a[i],a[j])=　0 otherwise スコアが最小でないものも二次構造とよび、最小のものを最適二次構造とよぶこともある {g,u} も塩基対に含まれる場合がある

5 RNA二次構造の二種類の表現

6 RNA二次構造の例 RNA配列

7 Nussinovアルゴリズム

8 予測アルゴリズム（Nussinovアルゴリズム）
入力配列： a=a[1]…a[n] 動的計画法 初期化 メインループ 最適解（＝ －塩基対の個数）

9 アルゴリズムの説明 メインループ

10 Nussinovアルゴリズムの解析 定理： 上記アルゴリズムは O(n3) 時間で最適解を計算
略証： テーブル W(i,j) のサイズはO(n2)。1個のテーブル要素の計算にO(n)時間（最後の行）。

11 RNA二次構造予測と 確率文脈自由文法

12 RNA二次構造予測と確率文脈自由文法 （１）
確率文脈自由文法（SCFG）： 導出確率が最大となる構文解析木を計算 ⇒ 確率の代わりにスコアを用いる score(Y)+score(Z) 文法表現としては X→aYu, X→XY などではなく、S→aSu, S→SS などが正式 RNAの場合はスコアは１ではなく、－１に置き換えることが必要

13 RNA二次構造予測と確率文脈自由文法 （２）
スコア最大（≒確率最大）の構文解析木 ⇔ 最適二次構造 実際、NussinovアルゴリズムはCYKアルゴリズム（文脈自由文法の構文解析アルゴリズム）に類似

14 Valiantアルゴリズムの利用 [Akutsu: J. Comb. Opt. 1999] [Zakov et al.: Alg. Mol. Biol. 2011]

15 二次構造予測の計算量の改良 文脈自由文法の構文解析 高速行列乗算に基づくRNA二次構造予測
高速行列乗算に基づく Valiant アルゴリズムを用いれば O(nω) 時間 O(nω) は n×n の行列乗算にかかる計算時間 高速行列乗算に基づくRNA二次構造予測 基本的に Valiant アルゴリズムを適用 しかし、行列乗算の基本演算を (+,×) から (max,+) に変える必要 (max,+) の行列演算は、ほんの少し O(n3) より良くなるだけ O(n3((log log n)/(log n))1/2)時間 [Akutsu: J. Comb. Opt. 1999] O(n3(log3(log n))/log n)時間 [Zakov et al.: Alg. Mol. Biol. 2011] SCFGの内側アルゴリズム[Akutsu: J.Comb.Opt. 1999]、外側アルゴリズム[Zakov et al.:Alg. Mol. Boil. 2011] は(+,×)演算で済むので O(nω) 時間で可能 Four-Russian アルゴリズムに基づくRNA二次構造予測 O(n3/log n)時間 [Frid, Gusfield: Proc. WABI 2009] ω は二十数年ぶりに から へ、さらに、2.327 へ改善された　[Wiliianms: Proc. STOC 2012]

16 Valiantアルゴリズムの概略（１） 基本的に分割統治 W(i,j)=Σk W(i,k)×W(k+1,j) の計算（青のベクトルと赤の
　　ベクトルの乗算）を 　　高速化 行列乗算を適用する 　　ため、複数の W(i,j) 　　の計算をまとめて 　　実行 左下三角は計算不要

17 Valiantアルゴリズムの概略（２） 基本戦略 白と黄が計算済みとして、青を計算（青は一部計算済み）
ピンクの部分の行列積を計算後、青と加算（⇒結果は緑） 緑と黄色からなる行列を作り、再帰計算により青を計算 　　　　　　　　　　　　（左下の白は不要なのですべて０としてOK） 高速 乗算 再帰計算

18 Valiantアルゴリズムの概略（３） 時間計算量 T2(n)= T2(n/2)+ 2T3(3n/4)+ T4(n)+ O(n2)
T3(n)= M(n/3)+ T2(2n/3)+ O(n2) T4(n)= 2M(n/4)+ T2(n/2)+ O(n2) ⇒ T2(n)= 4T2(n/2)+ 4M(n/4)+ O(n2) M(n)はn×n行列の乗算の計算量

19 Valiantアルゴリズムの概略（４） メインルーチン： 下図のとおり 時間計算量：

20 まとめ RNA二次構造予測 補足 動的計画法により O(n3)時間 Valiant アルゴリズムなどの利用により少しだけ改善
Valiantアルゴリズムは、RNAアラインメント・構造同時予測問題や結合RNA二次構造予測問題にも適用可 [Zakov et al.: Alg. Mol. Biol. 2011] 擬似ノットなしRNA二次構造予測の（クリークによる）下限は Ω(nω-ε) 時間[Abboud et al.: FOCS 2015]。現状では上限 とのギャップ。 ⇒　肯定的に解決された(O(n2.8…) time) [Bringmann et al.: FOCS 2016]

21 RNA二次構造予測の O(n2.8…polylog(n))時間アルゴリズム
原論文：　[Bringmann et al.: FOCS 2016]

22 背景 Valiant アルゴリズムにより二次構造は計算可能
ただし、行列乗算において、(×,+)を(+,max)に置き換える必要がある ⇒ 高速乗算が利用できない 一般の行列に対して、(+,max)（ (+,min)も等価）の行列乗算がO(n3-ε)時間で計算できるかは未解決 一方、行列の要素の絶対値が M以下であれば、 　　　　　時間で計算できることは既知(#1) [Alon et al.:JCSS, 1997] RNA二次構造予測では行列に bounded difference (BD) という制約（隣接要素の差異が小さい（W 以下））が存在(#2) #1と#2の性質を利用し、さらに乱拓アルゴリズムを巧みに用いることにより、BD行列(+,min)乗算の 　　　　　　時間を達成 さらに、脱乱拓可が可能（ 　　　　　　時間） は を表す ω は通常((×,+))の高速行列乗算の計算量のべき　(ω≦2.3729)　

23 アルゴリズムの概略 Phase 1 Phase 2: 以下をρ回繰り返す Phase 3 (C=A・B を計算)
（i,k,j）がΔの整数倍の要素についてのみ計算を行い、C=A・B の近似 D を計算 ⇒ すべての要素の誤差が 4ΔW 以内 Phase 2: 以下をρ回繰り返す ランダムに(ir,jr)を選び、絶対値が48ΔW以下の要素からなるAr, Br を作成し、 Ar・Br の結果をもとに D を更新 　　（　　　　　　時間行列乗算を利用） ⇒ 正しく計算されない要素の個数は Phase 3 正しく計算されなかった(i,k,j)（ただし、Δの倍数）に対し、 それを含む近傍（Δ×Δ×Δ）について再計算

24 Phase 1 I(i’) := {i ∊ [1..n] | i’-Δ < i ≦ i’}
補題１：　|Ci,j – Di,j| ≦ 4ΔW n=8, Δ=2

25 Phase 2: (+,max)行列乗算に関する性質
A’: 行列 Ai,k の各列 k の要素に同じ数 Fk を加算 B’: 行列 Bk,j の各行 k の要素から同じ数 Fk を減算 すると、A’・B’=A・B が成立 Fk をうまく定義することにより、（重要な部分について）要素の絶対値の小さな行列を得る　⇒ 　　　時間行列乗算を適用

26 Phase ２：アルゴリズム 基本的に が成立するので、 が、min{Ai,k+Bk,j} を与える k について成立すれば
より、Ci,j　を正しく計算 （少し余裕をみて）どの程度の個数の (i,k,j) についてρ回のうち１回も が成立しないかを見積もる

27 Phase ２：解析 補題２ (i,k,jr), (ir,k,jr), (ir,k,j) が weakly relevant であれば、
strongly relevant: k は mink を与える weakly relevant: mink となる k の候補の数の上限を与える 補題２ (i,k,jr), (ir,k,jr), (ir,k,j) が weakly relevant であれば、 (i,k,j) は weakly uncovered でない 補題３ weakly relevant で、かつ、weakly covered でない (i,k,j) の個数は、（非常に高い確率で） 証明の方針：補題２をもとに、（kを固定した場合の）weakly relevant な辺 (i,j)からなる二部 グラフにおける長さ４の閉路の個数に関する結果を利用

28 Phase 3：アルゴリズム 補題３から、 　　　　　　　　　　　個の (i,k,j) について、再計算すれば十分 　　　　　　⇒ しかし、そのような (i,k,j) の探索にO(n3)時間かかる そこで、まずΔの倍数の (i,k,j) についてチェックし、該当しそうなものについて、その近傍を再計算 ⇒ 少し条件を変える必要 近傍の再計算 (i,k,j)=(4,6,8) Δ=2

29 Phase ３：正当性 補題４: 各（i,k,j)について以下が成立 定義２: （i,k,j)∊
I(i’)×I(k’)×I(j’)について 補題４:　各（i,k,j)について以下が成立 strongly relevant ⇒ approximately relevant approximately relevant ⇒ weakly relevant strongly uncovered ⇒ approximately uncovered approximately uncovered ⇒ weakly uncovered よって、Phase 3により、 ・strongly relevant はすべて対応（ただし、既にstrongly uncoveredでないのは対応済み） ・approximately relevant の個数は weakly relevant の個数以下 ・strongly relevant で strongly uncovered であるものはすべて対応 ・（approximately relevant のうち） approximately uncovered の個数は weakly uncovered の個数以下 strongly relevant なものがもれなく計算されていればOK ⇒　Phase 3までで、すべて計算される

30 時間計算量の解析 Phase 1: Phase 2: Phase 3: + よって、全体で ここで、 とすると、 を得る
ここで、　　　　　　　　　　　　　　　　　　とすると、 　　　　　　　　　　　　　を得る W が O(1) であれば、この値は 最悪、O(n3)時間で済むので、平均計算時間も同様 定理 塩基対スコアの和に基づくRNA二次構造予測問題は （高確率および平均で）　　　　　　　時間で解ける

31 補題３：証明（１） 補題５：G=(U⋃V,E) を |U|=|V|=n, |E|=m である二部グラフとする。m≧2n1.5 であれば、長さ4の閉路数は m4/(32n4) 個以上 補題３の証明 z=c(n2/ρ)ln(n) とする。 Gk を E(Gk)={(i,j) | (i,k,j)は weakly relevant} により初期化。その後、weakly uncovered でなくなったら、その辺を削除。 (i,j) を含む長さ 4 の閉路 i→jr→ir→j→i が存在し、(ir,jr) が選択されれば、(i,j) は削除される z 個以上の閉路に含まれる (i,j) が削除されない確率は よって、c>4 とすれば高確率で、z 個以上の閉路に含まれる (i,j)は、すべて Gk から削除されている

32 補題３：証明（２） 証明（続き） mk を（最終的な） Gk の辺数とする。　　　　　　　　　　　　　　　（#）が成立するので、mk≧2n1.5 を満たす k のみを考える そのような各 k においては、１個の辺について z 以下の閉路しかない（それ以上であれば Gk から削除されているはず） よって、 Gk における閉路の総数は mkz 個以下。一方、補題５を考えると (mk/n)4/32 < mkz が成立。よって、以下が成立 よって、これをすべての k について足し合わせ、それに (#) を加えて、以下を得る 補題３ weakly relevant で、かつ、weakly covered でない (i,k,j) の個数は、（非常に高い確率で）

33 まとめ 行列の隣接値の差が定数以内であれば、(+,min)による行列乗算が（高確率および平均で） 時間で解ける
⇒ Valiant アルゴリズムとの組み合わせでRNA二次構造予測も同様の計算時間 さらに、これらの結果は、Phase 2, Phase 3の改良により 　　　　　　　時間まで改良されている また、その脱乱拓化もでき、　　　　　　時間の決定性アルゴリズムも開発されている クリークによる下限Ω(nω-ε) [Abboud et al.: FOCS 2015]とはギャップが存在。また、通常の行列乗算がO(n2)に近づいたとしても、計算量は　　　　　　　と議論されている