第5回 2019/7/8
- 講師:浅井 潔
- 資料:バイオインフォマティクス 確率モデルによる遺伝子配列解析 第5章
講義概要
今回の講義は、上記の資料(浅井先生が翻訳に関わられている本)の文章をそのまま追うだけであった。
ギャップなしのスコア行列
ここでは、タンパク質ファミリーのマルチプルアラインメントを考える。
5p21- MTEYKLVVV'G'AGGV"GK"'S'ALTIQLIQNH'F'VDEYD'P'TIEDSY
1ctqA MTEYKLVVV'G'AGGV"GK"'S'ALTIQLIQNH'F'VDEYD'P'TIEDSY
1c1yA MREYKLVVL'G'SGGV"GK"'S'ALTVQFVQGI'F'VEKYD'P'TIEDSY
1kao- MREYKVVVL'G'SGGV"GK"'S'ALTVQFVTGT'F'IEKYD'P'TIEDFY
1huqA --QFKLVLL'G'ESAV"GK"'S'SLVLRFVKGQ'F'HEYQE'S'TIGAAF
1g16A ----KILLI'G'DSGV"GK"'S'CLLVRFVE--'-'-DKFN'P'I--DFK
1ek0A VTSIKLVLL'G'EAAV"GK"'S'SIVLRFVSND'F'AENKE'P'TIGAAF
3rabA ---FKILII'G'NSSV"GK"'T'SFLFRYADDS'F'TPAFV'S'TVGIDF
1mh1- ----KCVVV'G'DGAV"GK"'T'CLLISYTTNA'F'PGEYI'P'TVFDNY
2ngrA MQTIKCVVV'G'DGAV"GK"'T'CLLISYTTNK'F'PSEYV'P'TVFDNY
1tx4B ----KLVIV'G'DGAC"GK"'T'CLLIVNSKDQ'F'---YV'P'TVFENY
上の配列から見てわかるように、タンパク質ファミリーの配列同士には特徴がある。
タンパク質ファミリーのマルチプルアラインメントの特徴
- ギャップは互いに一列に並ぶ傾向があり配列によって所々抜けていることが少ない。
- サイトごとに保存度に差がある。
- サイトごとにアミノ酸の出現傾向に差がある。
- サイトごとにギャップの出現傾向に差がある。
そこで、まずは上記の性質1を利用して、分断化される非ギャップ領域のアラインメントを考える。
まず、位置 \(i\) において観察されるアミノ酸 \(a\) が、各位置 \(i\) について独立の確率 \(e_i(a)\) に従って現れると仮定した、最も簡単なモデルを考える。
この時、このモデルに従った配列 \(x\) の確率は、
と表される。ただし \(L\) はブロック長で、上の配列の場合 \(40\) である。
ここで、次の位置特異的スコア行列(Position Specific Score Matrix; PSSM)
変数 | 意味 |
---|---|
\(e_i(a)\) | 位置 \(i\) においてアミノ酸 \(a\) が出力される確率 |
\(q_a\) | アミノ酸 \(a\) の背景確率 |
を考える。この値は、サイトごとのスコア行列を意味している。
スコア行列が \(s(a,b)=\log\frac{p(a,b)}{q(a)q(b)}\) で表されていたことを思い出す。これとの比較より、今考えているのはスコア行列の第2引数 \(b\) が位置 \(i\) である時のスコア行列と言うことができる。
挿入状態と欠失状態をプロファイル HMM に加える
ここまでで、位置 \(i\) に固有のスコア行列 \(S_i\) を定義することができた。
しかし、これだけではまだタンパク質ファミリーの特徴を全て表現しきれていない。次に、ギャップを考慮する方法を考える。
ペアワイズアラインメント同様に単一のギャップスコア \(\gamma(g)\) を導入し、各位置についてギャップを許すことも可能だが、サイトごとのギャップの出現傾向を捉えるために、位置に依存したギャップスコアを導入する。
ギャップが起こる原因に関しては、「挿入」と「欠失」が考えられるが、これらは分けて取り扱う方が望ましい。
そこで、各位置においてそれぞれの隠れ状態を持ったHMMであるPrpfile HMMを考える。
Profile HMMs
状態 | 意味 | 画像 |
---|---|---|
\(M_i\) | 一致状態 | |
\(I_i\) | マルチプルアラインメントの \(i\) 番目のカラムに一致した塩基の後の挿入を表現する。 | |
\(D_i\) | ヌル状態。塩基を出力しない。欠失を表現するが、欠失のコストは \(M\rightarrow D\) 遷移、それに続くいくつかの \(D\rightarrow D\) 遷移、\(D\rightarrow M\) 遷移の和になる。 |
補足
- \(I_i\) は出力分布 \(e_{I_i}(a)\) を持つが、それらは普通バックグラウンドの分布 \(q(a)\) が割り当てられる。
- \(D\rightarrow I\) や \(I\rightarrow D\) の遷移は実際の配列解析ではめったに使われることがなく、それらを取り除いても一致のスコアには大した影響は与えない。(しかし、モデルの構築上の問題が生じてしまう)
マルチプルアラインメントからのプロファイル HMM の導出
プロファイルHMMのパラメータ
モデルの構造が定義できたので、各種パラメータの調整を考える。
出力確率や遷移確率が非ゼロだと仮定すると、プロファイルHMMは与えられたアルファベットからの任意の塩基配列をモデル化することができる。
そこで、プロファイルHMMが表す分布をファミリーの一員付近でピークに達するように、パラメータを調整したい。
分布の形を調整するために利用できるパラメータは,「確率パラメータ」と「モデル長」である。これらを最適に設定する方法を数多くあるが、まずはKrogh, A. et al.1994による基本的な方法を紹介する。
モデル長
これは非常に簡単で、「ギャップが半数より多いカラムは挿入」としてモデリングする。
確率パラメータ
マルチプルアラインメントは与えられているので、直接パラメータを推定することができる。
つまり,各遷移や出力が使われた回数を数え上げ、次式にしたがって確率を割り当てる。
ここで、\(k,l\) は状態の添え字であり、\(a_{kl},e_k\) は遷移確率と出力確率、\(A_{kl},E_k\) は対応する頻度である。
問題点
学習データが多い場合、この確率パラメータの推定方法は正確で矛盾のない確率の推定方法となる。
しかし、配列の数が少ない場合、ある遷移や出力が学習アラインメントに現れないとその確率が0になり、そのパスが未知の配列に対しても決して現れないことを意味する。
これを防ぐよく使われる方法は、「擬似度数(pseudocount)を観察頻度に加える」方法である。これは、事前分布を導入することに対応する。
プロファイル HMM による検索
HMMに対する一致をスコア化するには、\(2\) 通りの方法がある。つまり、Viterbiアルゴリズムと前向きアルゴリズムを利用することができる。
- 前者では,配列 \(x\) の最も尤らしいアラインメント \(\pi^{+}\) が、その確率 \(P(x,\pi^{+}|M)\) と共に得られる。
- 後者では,全ての可能なパス \(P(x|M)\) に渡って足し合わされた \(x\) の確率が得られる。
なお、いずれのケースでも、実用的な目的から、一致の候補を評価する時に考えるべき値は、標準的なランダムモデルが与える \(x\) の確率
に対する対数オッズ比である。
Viterbiアルゴリズム
動的計画法を回すために、以下の変数を定義する。
変数 | 意味 |
---|---|
\(V_j^M(i)\) | 状態 \(M_j\) によって出力された \(x_i\) で終わる最良パスのスコア |
\(V_j^I(i)\) | 状態 \(I_j\) によって出力された \(x_i\) で終わる最良パスのスコア |
\(V_j^D(i)\) | 状態 \(D_j\) によって出力された \(x_i\) で終わる最良パスのスコア |
すると、次のような漸化式を導出することができる。
漸化式
初期条件
ここで気をつけるべきことは、「アラインメントが欠失状態や挿入状態から開始したり終了したりできるようにすること」である。
このことを機械的に保証する最も簡単な方法は、開始状態を \(M_0\) とし、\(V_0^M(0) = 0\) とすることである。(\(M_0\) から \(I_0\) と \(D_1\) への遷移は許す。)
前向きアルゴリズム
前向きアルゴリズムのための再帰式はViterbi式と似ているが、\(\max()\) 操作が足し算によって置き換えられている。したがって、以下の変数を用いる。
変数 | 意味 |
---|---|
\(F_j^M(i)\) | \(V_j^M(i)\) に対する全ての可能性が足し合わされた対数オッズ比 |
\(F_j^I(i)\) | \(V_j^I(i)\) に対する全ての可能性が足し合わされた対数オッズ比 |
\(F_j^D(i)\) | \(V_j^D(i)\) に対する全ての可能性が足し合わされた対数オッズ比 |
すると、次のような漸化式を導出することができる。
漸化式
初期条件
初期化に関しては、\(F_0^M(0) = 0\) と初期化し、Viterbiと同じように行えば良い。
対数オッズ以外のスコア
HMM に関する初期の論文では、対数オッズスコアよりも、モデルが与える配列の対数尤度が直接利用されていた。(LLスコアと呼ばれている。\(LL(x) = \log P(x|M)\))
LLスコアは当然ながら配列長に依存するため、直接比較するには適当ではない。(配列長で割った値を使う手法もあるが、これもLLと配列長の関係が線形でないことから厳密な方法ではない。)
これを回避する方法は、配列長の関数として平均スコアと標準偏差を推定し、各配列の標準偏差を利用することである。これはZ-スコアと呼ばれる。
Z-スコアを計算するために、LLスコアまたは対数オッズスコアに滑らかな曲線を適合させる。(詳細はKrogh, A. et al.1994を参照)
そして、各配列長 (あるいは配列長周辺の短い区間) について標準偏差を推定し、標準偏差を単位として滑らかな曲線からの距離を各スコアについて計算する。これがZ-スコアである。
一致の検出
今や各配列がタンパク質ファミリーに属するかの閾値を見つけることが可能になった。
また、判別には対数オッズに基づくスコアがより適当であることも明確である。(∵ランダムモデルの確率で割ることによって、配列の残基組成に合わせた調節が行なわれているから。これを行なわないと、単純に頻出する塩基で組成された配列のスコアが高くなり、ノイズの分散を増加させることになる。)
アラインメント
一致の検出から離れると、プロファイルHMMの他の主要な使い道は、ファミリーに対して配列をアラインメントすること、より正確には、その配列をファミリーのマルチプルアラインメントに加えることである。
非大域アラインメントの変形としてのプロファイル HMM
ここでは、プロファイルHMMへの配列のViterbiアラインメントと、アフィンギャップコストを用いた2本の配列の大域DPによる配列比較に存在する密接な関連性を見た。(省略)
確率推定についての再考
ここでは、パラメータ推定について深掘りし、網羅的に紹介する。
なお、ここでは出力確率についてしか触れないが、類似の手法を遷移確率に対しても用いることができる。
最尤推定
最も素直なアプローチは、パラメータの最尤推定だと思われる。
少し表記を変えるが、アラインメントの位置 \(j\) における残基 \(a\) の観察頻度 \(c_{ja}\) について、対応するモデルパラメータ \(e_{M_j(a)}\) の最尤推定は次のようになる。
この方法に関する明らかな問題点は、訓練データの中にある結果が観察されていなければその確率が \(0\) と推定されてしまうことである。
この問題に対する最も簡単な方法は先に述べた通り、観察度数 \(c_{ja}\) に擬似度数を加えることである。
単純擬似度数
非常に簡単でよく用いられる擬似度数は、全ての度数に一定値を足し合わせることである(一様分布を事前分布に用いる)。
その一定値が \(1\) の場合を特別にLaplaceルールと呼ぶ。
少しだけ洗練された方法は、バックグラウンドの分布に比例した量を加えることである。
ここで、\(c_{ja}\) は実際の度数、\(A\) は実際の度数と比較して擬似度数に加える重みである。\(20\) 辺りの \(A\) の値がタンパク質のアラインメントには適しているとの意見もある。
この正規化において
- 利用可能なデータが非常に少ない場合、すなわち\(A\) と比較して全ての実度数が非常に小さい場合に、\(e_{M_j(a)}\) が近似的に \(q_a\) となる。
- 利用可能なデータが大量にある場合、この正規化の効果は有意ではなくなり、\(e_{M_j(a)}\) は本質的に最解と同じになる。
この擬似度数法は、ベイズの枠組みでは、パラメータ \(\alpha_a = Aq_a\) のDirichlet事前分布を仮定していることに対応している。
Dirichlet混合分布
この手法は、Brown, M. et al.1995によって提案された。
基本的なアイデアは、アラインメントの環境の種類に対応して、いくつかの異なる組の擬似度数 \(\alpha^1,\ldots,\alpha^K\) が存在するべきだろうというものである。
例えば、ひとつのパラメータの組みは、表面のループ環境に関係するもの、ひとつの組みは埋没した小さい残基の環境に関するもの、…といった具合いである。
与えられた度数 \(c_{ja}\) について、まず各事前分布がどれくらい観察データに適合しているかに基づいて尤もらしさを推定し、それらの事後確率にしたがってそれらの結果を結合する。
但し、\(P(k|c_j)\) は事後混合係数(posterior mixture coefficient)であり、ベイズの法則によって計算できる。
ここで、\(p_k\) は各混合要素の事前確率、\(P(c_j|k)\) はDirichlet混合分布 \(k\) に従った確率である。
置換行列の混合
Dirichlet分布の混合を利用した別のアプローチは、観察頻度と置換行列を用いて単一のDirichlet定式化における類似度数を調節することである。
これは理論的に的確に確かめられたアプローチではないが、非確率プロファイル法とDirichlet擬似度数法の特徴を結合するという直感的な意味を持つ。
最初のステップは、置換行列の要素 \(s(a,b)\) を条件付き確率 \(P(b|a)\) に変換することである。
これは、置換行列の要素が対数オッズ比として導出されている場合は \(s(a,b) = \left(\log P(a,b)/q_aq_b\right) = \log\left(P(b|a)/P(b)\right)\) なので、\(P(b|a) = q_be^{s(a,b)}\) となる。(対数オッズ比として導出されてない場合は後で紹介)
すると、与えられた条件付き確率 \(P(b|a)\) について
によって \(\alpha_{ja}\) を求め、これを用いて
とすれば出力確率が求まる。
この手法では、「アミノ酸 \(i\) の擬似度数 \(j\) への寄与は、カラムでのそのアミノ酸の度数とアミノ酸 \(j\) へ変化する確率に比例する」としているが、これは到って直感的なアイデアである。
任意の行列からの \(P(b|a)\) の導出
スコア行列 \(s(a,b)\) が対数オッズ行列として導出されていなくても、条件
- 行列が負に偏っていること:\(\sum_{ab}q_aq_bs(a,b)<0\)
- 少なくともひとつの正の要素を含んでいること
が満たされている限り、対数オッズ行列として解釈すると \(\lambda s(a,b)\) が正しく振る舞うようなスケール因子 \(\lambda\) を見つけることができる。
ここで求めたいのは、
についての値 \(r_{ab}\) の組である。但し、\(r_{ab}\) は \(a,b\) ペアについての確率として解釈することができる。この式は簡単に逆をとることができ、\(r_{ab} = q_aq_be^{\lambda s(a,b)}\) となる。
確率として正しくするために \(r_{ab}\) の和は \(1\) としなければならないため、
を満たす \(\lambda\) を見つける必要がある。
そのような \(\lambda\) のひとつは \(\lambda = 0\) だが、これは明らかに求めたいものではない。この式において、先に与えた二つの条件が別の明確な解の存在を保証する。
- 負への偏りの条件を用いれば、\(f(\lambda)\) が十分小さい \(\lambda\) について負であることが示せる。(\(f'(0)\) を計算し、\(\lambda =0\) について \(f(\lambda)\) を導出する)
- 少なくともひとつは存在する正の \(s(a,b)\) を用いて、\(f(\lambda)\) が十分大きな \(\lambda\) について正になることを示す。
- \(f(\lambda)\) の2次の導関数が正であることと1,2より、上式を満たす \(\lambda\) の値が \(\lambda=0\) の他にただひとつ存在することを示す。
得られる \(\lambda\) の値は、置換行列の自然スケーリング因子と呼ばれる。置換行列を確率的に捉えることによって \(\sum_{ab}r_{ab}\log(r_{ab}/q_aq_b)\) という行列についてのエントロピースコアが導かれる。これは、置換行列を特徴付けたり比較するのに便利な量である。
進化的祖先に基づいた推定
これまでに説明したものより原理的で直接的な手法として、「観察配列全てが共通の祖先から独立に派生していると仮定」し、その共通祖先の与えられた位置に存在する残基を推定することで、そこからその祖先の子孫について、各残基が観察される確率を推定する、という方法がある。
アラインメントのカラム \(j\) における残基 \(x_j^k\) を持つサンプル配列 \(x^k\) があると仮定する。(ギャップがある場合 \(x_j^k\) は配列 \(x^k\) における \(j\) 番目の残基ではないことに注意)
ここでも、置換行列から導出された条件付き確率 \(P(a|b)\) が必要となる。共通祖先における残基を \(y_j\) とすると、\(y_j = a\) である事後確率を計算するためにベイズの法則を利用できる。
※単なる塩基配列しか情報がない場合、共通祖先における残基の事前分布はアミノ酸のバックグラウンド確率になるため、その分布が必要になることに注意!!
この値を用いれば、HMMの出力確率を新しい配列に対する予測確率として計算することができる。
このアプローチに関するひとつの問題は、やはり「単一の置換行列を利用するため、カラム毎の保存度合いの違いを表現できない」ことである。
この問題に対する解決策として、置換行列を複数用意し、保存度合いの違いに応じて異なる置換行列を選択することである。
行列の選択の仕方は観察データの尤度を最大化する方法が最も直感的である。
ベイズ流のアプローチを用いることもできる。行列ファミリーパラメータ \(t\) を用いて
とし、これと \(t\) についての事前分布を用いて \(t\) の事後分布を得、それを先の \(e_{M_j(a)}\) の式で足し合わせることも可能であるが、それだと非常に多くの計算量が必要になってしまう。
講義はここまで。実際の本には以下の章が続く。
- 5.7 最適モデルの構築
- 5.8 学習配列の重み付け