第3回 2019/6/24

• 講師：寺井 悟郎

Table of contents

1. RNA secondary structure
2. Nussinov algorithm
3. Zuker algorithm
4. McCaskill algorithm

1. RNA secondary structure

遺伝子の実態はタンパク質であり、RNAはDNAとタンパク質の橋渡しをしているだけだと考えられていたが、ncRNAなど、RNA単独でも様々、かつ重要な機能を果たすことがわかってきた。

例えば、ncRNAのうち長さが21~23塩基程度のmiRNAは、まずヘアピン構造を持つ長いRNA(pri-miRNA)として転写される。

その後、ヘアピン構造の基部を切断するRNA分解酵素Droshaにより、中間産物のpre-miRNAが切り出される。pri-miRNAはさらにDicerによって切断されて二本鎖のmiRNAとなった後、その片方、もしくは両方がmature miRNAとして機能する。

このように、miRNAが正しいプロセシングを受けるためには、特徴的な構造が不可欠である。

また、mRNAの結合部位にヘアピン構造などの二次構造があるとRISCが近づけず、RNAiが起こらないことも知られている。

これらの機能を理解するためにも、RNAの構造予測や構造比較が重要になっている。

しかし、RNAの3次元構造決定はタンパク質よりも困難な場合が多く、既知のデータも少ないので、現在の研究の中心は、どの塩基とどの塩基が立体構造中で結合しているかを表現した二次構造と呼ばれる一種のグラフ構造の予測や比較である。

擬似ノット(pseidoknot)

これ以降3つの2次構造予測アルゴリズムを紹介するが、それらはどれも擬似ノット(pseidoknot)なしの2次構造のみを対象とすることに注意が必要である。

擬似ノットを含む2次構造の予測についてはこの論文が詳しそう。

定義

以下のいずれかを満たす2個の塩基対 \((i,j),(i',j')\in A\) は擬似ノットと呼ばれる。 - \(i<i'<j<j'\) - \(i'<i<j'<j\)

擬似ノットがある場合、少なくとも1つの対は交差することになる。

2. Nussinov algorithm

RNA二次構造予測問題は「考えられる様々な2次構造から生物学的にもっとも正しい構造を選ぶ問題」であり、一般にエネルギー最小化問題として定義される。（細胞内では、RNAはエネルギー的に安定な構造をとると考えられるから。）

Nussinov algorithmもその一つであり、「結合する塩基対はどれも同じ程度の負のエネルギーを持つ」という最も簡単な仮定をおく。（つまり、塩基対が多ければ多いほど生物学的に正しいであろうという仮定）

The basic idea

より短い部分配列の最適構造を伸ばしたり組み合わせることで目的の配列の最適構造を求める。以下の変数を定義する。

ここでDP matrixを用意するが、このDP matrixの斜めの線は、長さ \(n\) の部分配列の最大塩基対数を表している。

Algorithm

Implementation: link

アルゴリズムは、以下のようにかける。

• Initialization:
  $$ \begin{array}{l}{\gamma(i, i)=0 \text { for } i=1 \text { to } L} \\ {\gamma(i, i-1)=0 \text { for } i=2 \text { to } L}\end{array} $$
  
• Recursion:
  $$ \gamma(i, j)=\max \left\{\begin{array}{c} {\gamma(i+1, j)} \\ {\gamma(i, j-1)} \\ {\gamma(i+1, j-1)+\delta(i, j)} \\ {\max _{i \leq k<j}[\gamma(i, k)+\gamma(k+1, j)]} \end{array}\right. $$

Sample

※ G-U can form a base-pair.

この計算では、一つの要素を決めるのに \(O(L)\) の計算量がかかり、それが \(O(L^2)\) あるので、計算量は \(O(L^3)\) である。一方メモリは要素数に比例するので、\(O(L^2)\) である。

Bifurcation

上の再帰式の最後の行は Bifurcationと呼ばれ、「ある部分配列とある部分配列の最大塩基対数の和が全体配列の最大塩基対数になる時」に非常に重要な役割を果たす。

Traceback

トレースバックが複数ある場合は、最適構造が複数あることを意味している。以下のトレースバックアルゴリズムでは、そのうち一つのみを求めることができる。

#=== Initialization ===
stack = []
record = []
stack.append((1,L))
#=== Recursion ===
while(len(stack)>0):
  i,j = stack.pop()
  if (i >= j): continue
  elif γ(i+1,j) == γ(i,j): stack.append(i+1,j)
  elif γ(i,j-1) == γ(i,j): stack.append(i,j-1)
  elif γ(i+1,j-1) + δ(i,j) == γ(i,j):
    record.append((i,j))
    stack.append(i+1,j-1)
  else:
    for k in range(i,j):
      if γ(i,k) + γ(k+1,j) == γ(i,j):
        stack.append(k+1,j)
        stack.append(i,k)
        break

3. Zuker algorithm

Zuker algorithmでは、Nussinov algorithmの時より複雑にエネルギーを計算することになる。

二次構造のエネルギーは、ループ(loop)と呼ばれる部分構造（塩基対で閉じられた構造）の自由エネルギーの和として近似される。

例えば上の図形だと、以下のような様々なループが存在し、それぞれのエネルギーの総和で全体のエネルギーが表されるとする。

• 青色のhairpin構造は一つの塩基対で閉じられたループ
• 黄色のbulge構造は二つの塩基対によって閉じられたループ
• シアンのmulti-branched構造は三つの塩基対によって閉じられたループ

ここで、それぞれのループの自由エネルギーを以下のように定義する。

• \(k\) は塩基対を形成している塩基対の数
• \(u\) は塩基対を形成していない塩基の数

\(a,b,c\) は定数であり、またそれ以外の \(F_1,F_2\) という関数も予め実験的に求まっているものとする。（ex. \(F_1(C-G) = -3.4 \mathrm{kcal/mol},F_1(A-U) = -0.9 \mathrm{kcal/mol}\)）

Algorithm

Implementation: link

Zuker algorithm は基本的にはNussinov algorithmと似ているが、以下の4つのDP matrixを持つ。

Recursion

\(W(i,j)\)

\(V(i,j)\)

続いて \(V(i,j)\) の計算についてだが、\((i,j)\) 塩基対がどの構造を作るかによって場合分けが必要であり、以下のようになる。

\(M(i,j)\)

\(M(i,j)\) はマルチループを考えるため、\(i,j\) 間に2つの「塩基対に閉じられた部分構造」があることを保証する必要がある。したがって、\(M1(i,j)\) という変数（\(i,j\) 間に2つの「塩基対に閉じられた部分構造」がある）を用意している。

\(M1(i,j)\)

以下の再帰式の意味を解釈する。ここで、マルチループのエネルギー関数が

• \(k\): 塩基対を形成している塩基対の数
• \(u\): 塩基対を形成していない塩基の数

で計算できたことを思い出す。

• \(V(i,j)+b\) は、\(i,j\) が塩基対を組んでいる場合を仮定しているので、「塩基対を形成している塩基対の数」に比例したエネルギー係数 \(b\) を足している。
• \(M 1(i+1, j)+c\) は、\(i\) を1つ外に拡大した（\(i+1\rightarrow i\)）場合を考えている。この時 \(i\) 番目の塩基は「塩基対を形成していない」ので、\(c\) を足している。

最後に、\(V(i,j)\) の式で \(M(i,j)\) は「\(M(i+1,j-1) + a + b\)」の形で現れるが、これは \(i,j\) が塩基対を組んでいるので \(b\) を、マルチループ特有のエネルギーとして \(a\) を足していると考えれば理解できる。

Initialization

Traceback

※なお、講義ではマルチループを考えない（ \(M,M1\) を取り除く）場合を考えたので、ここでの説明は \(W,V\) のみのDPmatrixがある場合となっている。

Nussinov algorithmのトレースバックと似ているが、どちらの行列をたどっているかを記憶する変数が必要であり、さらに行列 \(V\) のトレースバックを新たに考える必要がある。

#=== Initialization ===
stack = []
record = []
stack.append((1,L,"w"))
#=== Recursion ===
while(len(stack)>0):
  i,j,m = stack.pop()
  elif(m=="w"):
    if W(i+1,j) == W(i,j): stack.append(i+1,j,"w")
    elif W(i,j-1) == W(i,j): stack.append(i,j-1,"w")
    elif W(i,j) == V(i,j): stack.append(i,j,"v")
    else:
      for k in range(i,j):
        if W(i,k) + W(k+1,j) == W(i,j):
          stack.append(k+1,j,"w")
          stack.append(i,k,"w")
          break
  elif(m=="v"):
    # hairpin loop
    if V(i,j) == F_1(i,j): record.append((i,j))
    else:
      flag = 1
      for h in range(i+1,j):
        for l in range(h+1,j):
          if V(i,j) == F_2(i,j,h,l) + V(h,l):
            record.append((h,l))
            stack.append(h,l,"w") # w??
            flag = 0
            break
      if flag==1:      
        for k in range(i+2,j-2):
          if V(i,j) == W(i+1,k) + W(k+1,j-1):
            satck.append(i+1,k,"w")
            stack.append(k+1,j-1,"w")
            break

ここまで、与えられた配列の二次構造を予測してきた。

4. McCaskill Algorithm

ここからは、与えられた配列の塩基対確率行列を計算する方法を学ぶ。

ex.) "GGGAAAUCC" という配列を考える。すると、以下の \(20\) 通りの配列を取り得る。（ワトソンクリック対とWobble対のみが対合できるという条件下で）

まず、簡単のためにそれぞれの配列になる確率が等確率であるとする。

ここで、2番目の"G"と8番目の"C"が対合する配列は上の赤で示した4通りである。したがって、「2番目の"G"と8番目の"C"が対合する確率は \(4/20\) 」 と結論付けることができる。

もちろんこの仮定はあまりにも単純すぎであり、本来は異なるエネルギーを持つ構造のうち、安定なものを取りやすいはずである。

したがって、それぞれの構造ごとにエネルギー \(E(\sigma_i)\) を割り当て、それに従って \(P(\sigma_i)\) が計算される、とするべきである。

Boltzmann distribution

\(E(\sigma_i)\) から \(E(\sigma_i)\) を計算する式が必要になるが、RNAの二次構造の確率は、ボルツマン分布(Boltzmann distribution)に従うと仮定する。

この時、上の式の変数はそれぞれ以下の意味を持つ。

• \(e^{-E(\sigma) / k T}\) をボルツマン因子(Boltzmann factor)
• \(k\) を定数
• \(T\) を温度(Temperature)

※エネルギーがボルツマン分布に従うという仮定は物理の世界ではかなり一般的な概念であり、粒子の集合が合計で \(Z\) のエネルギーを持つ時に、個々の粒子にエネルギーが等確率で分配されるというモデルを考えるとわかる。

従って、エネルギーから確率への変換は、以下のように表すことができる。

DP matrix

短い配列の場合は先ほどの表のように全ての取り得る配列を列挙することができるが、配列が長くなるとかなり厳しい。

そこで、これをMcCaskill algorithmで今まで同様に動的計画法によって求める。

なお、McCaskill Algorithm の計算量は \(O(L^3)\) で、メモリは \(O(L^2)\) 必要で、Nacinohu Algorithm や Zuker Algorithm と等しい。

注意点

値の丸め込み

McCaskill Algorithm では、「ある2次構造のエネルギーはその部分構造のエネルギーの和として求めることができる。」とするが、この時以下のように掛け算が出てくる。

すると、DP行列に入る値が非常に大きく or 小さくなってしまう。そこで、ログを取るもしくはスケーリングする必要が出てくる。

2次構造のダブルカウント

分配関数 \(Z\) は、「全ての取り得る2次構造のボルツマン因子の和」として定義されていた。従って、同じ2次構造のボルツマン因子を二度足し合わせることはできない。

部分構造のボルツマン因子から全体構造のボルツマン因子を計算する際には、「異なる部分構造の分け方が同じ二次構造のボルツマン因子になり得る」ということに注意が必要である。

Nacinohu Algorithm や Zuker Algorithm では \(\mathrm{Max}\) や \(\mathrm{Min}\) を取っていたためこのようなダブルカウントは問題にならなかったが、今回のアルゴリズムでは注意しなくてはならない。

Algorithm

McCaskill algorithm では、5つのDP matrixを考えることになる。

各変数の役割

再帰式

具体的な解説

ここから、一つずつそれぞれの役割や、再帰式の意味を詳しく見ていく。

\(Z_{ij}\)

• 第1項: \(i\) から \(j\) までが一つも塩基対を組まないような2構造のボルツマン因子に対応。塩基対が一つもなければエネルギーは \(0\) なので、\(e^0=1\) となる。
• 第2項: \(i\) から \(k\) までの部分配列の全ての2次構造のボルツマン因子の和に、\(k+1\) がその下流の塩基と塩基対を組んで、そのさらに下流には塩基対が存在しない構造のボルツマン因子を掛け合わせたもの。

\(Z_{i j}^{1}\)

\(i\) がその下流の塩基と塩基対を組んで、さらにその下流には塩基対が存在しない構造のボルツマン因子の和。

• \(i\) が下流と形成する塩基対で閉じられた部分構造のボルツマン因子が \(Z_{ik}^b\) に対応する。その下流の部分のエネルギーは \(0\) なので、\(e^0=1\) がかけられているが、明示的には式に現れない。

\(Z_{ij}^b\)

\(i\) と \(j\) が塩基対を組む時の全ての可能な2次構造のボルツマン因子の和であり、Zuker Algorithm の \(V(i,j)\) に対応。

• 第1項: 塩基対 \((i,j)\) がヘアピンループを閉じる時のボルツマン因子
• 第2項: 塩基対 \((i,j)\) がstacking,bulge loop, internal loopを閉じる時のボルツマン因子
• 第3項: 塩基対 \((i,j)\) がmulti-loopを閉じる時のボルツマン因子。\(i+1\) から \(k-1\) までの \(Z^m\) と \(k\) から \(j-1\) までの \(Z^{m1}\) を掛け合わせたもの。\(e^{-[(a+b) / k T]}\) は、\((i,j)\) がmulti-loopを閉じる時のエネルギーの寄与分を表す。

\(Z_{i j}^{m}\)

\(i\) から \(j\) までがmulti-loopの中にあって、かつ塩基対で閉じられた部分構造を1つ以上持つ構造のボルツマン因子の和。

• 第1項: \(i\) から \(j\) までがmulti-loopの中にあって、かつ塩基対で閉じられた部分構造を1つだけ持つ構造のボルツマン因子の和。
• 第2項: \(i\) から \(j\) までがmulti-loopの中にあって、かつ塩基対で閉じられた部分構造を2個以上持つ構造のボルツマン因子の和。

\(Z_{i j}^{m 1}\)

\(i\) から \(j\) までの部分配列がmulti-loopの中にあり、\(i\) が下流の塩基と塩基対を形成し、さらにその下流の塩基が塩基対が存在しない構造のボルツマン因子の和。

• \(e^{-[b / k T]}\) は、multi-loopの中にあって、multi-loopを閉じる塩基対一つ一つにかかるエネルギーの寄与分を表す。
• \(Z_{i k}^{b}\) は、\(i\) がその下流の塩基と塩基対を組んだ時に、その塩基対に閉じられた部分構造のボルツマン因子の和。
• \(e^{-[c(j-k) / k T]}\) は、さらに下流にある部分のエネルギーの和。

ここまで、分配関数 \(Z\) の計算の説明を行った。

塩基対確率の計算

続いて、塩基対確率の計算方法の説明に移る。

ボルツマン分布の仮定のもとでは、\(i\) 番目と \(j\) 番目の塩基対確率は

によって表すことができた。なお、\(\mathrm{S}(i, j)\) は、\(i\) と \(j\) が塩基対を組むような全ての二次構造の集合を示す。

配列の長さが長くなると、この集合 \(\mathrm{S}(i, j)\) のサイズが大きくなり、数え上げることができなくなる。

そこで、数え上げなくても効率的に計算することのできる方法が必要になる。

これを実現するために、配列全体を \((i,j)\) 塩基対を基準にして内側と外側に分けて考える。そして、内側のボルツマン因子と外側のボルツマン因子を計算して掛け合わせることで全体のボルツマン因子を計算する。

内側

ここで、\((i,j)\) が塩基対を組む時の内側部分のボルツマン因子は、先ほど分配関数を計算する過程で現れた \(Z^b_{i,j}\) そのものになる。

外側

一方、外側部分を \(W^b_{i,j}\) とする。これを正確に言うと「\((i,j)\) が塩基対を組む時の、\(i\) から \(j\) までの部分配列を除いた外側配列の、全ての可能な2次構造のボルツマン因子の和」になる。

これは、以下の式によって求めることができる。

• 第1項: \((i,j)\) 塩基対の上流と \((i,j)\) 塩基対の下流で塩基対を形成しないような全ての可能な2次構造のボルツマン因子の和。
• 第2項: \((i,j)\) 塩基対の1つ外側にある \((h,l)\) 塩基対がstackingを閉じる場合のボルツマン因子の和。
• 第3項: \((i,j)\) 塩基対の1つ外側にある \((h,l)\) 塩基対がmult-loopを閉じる場合のボルツマン因子の和。なお、大括弧の中身は、以下を意味する。
  • 第1行: 左側だけに部分構造がある場合
  • 第2行: 右側だけに部分構造がある場合
  • 第3行: 両側に部分構造がある場合

結論

以上により、以下の式で塩基対確率を求めることができる。