2013-04-06

Generalized Neighbor Joining for Phylogenetic Tree and Algebraic Statistics

代数統計学近隣結合法 NJ法

こちらの本の第18章(全体構成などはこちら)

Algebraic Statistics for Computational Biology

作者: L. Pachter,B. Sturmfels
出版社/メーカー: Cambridge University Press
発売日: 2005/08/22
メディア: ハードカバー
購入: 1人クリック: 9回
この商品を含むブログを見る

Neighbor Joining Methodの復習
- $k$ 個の対象(種とか)がある
- $\begin{pmatrix} k\\2 \end{pmatrix}=\frac{k(k-1)}{2}$ ペアについて(DNA配列を基に距離(非負の実数)を対応づける
- $\begin{pmatrix}k\\2 \end{pmatrix}$ 個の情報から「木」を再構成する
- 再構成するときの手順(アルゴリズム)にNJ法という名前がつけられている
- ちなみに、木における２個の対象間の距離はパスの長さで定義している
NJ法を一般化する
- $k$ 個の対象(種とか)がある
- $\begin{pmatrix} k\\m \end{pmatrix}=\frac{k!}{m!(k-m)!}$ 組合せについて(DNA配列を基に距離(非負の実数)を対応づける
- $\begin{pmatrix}k\\m \end{pmatrix}$ 個の情報から「木」を再構成する
- 再構成するときの手順(アルゴリズム)にGeneralized NJ法という名前をつけることにする
- ちなみに、木における２個の対象間の距離はパスの長さで定義している
ここで代数統計学が出てくるのは次の理由から
- $\begin{pmatrix}k\\m \end{pmatrix}$ 個の情報から、最尤な木を得ようと思うと、 $\begin{pmatrix}k\\m \end{pmatrix}$ 個の情報がうまく木にフィットしているときは、(まだ)何とかなるが、ずれがあると、数値解法ではうまくいかなくなりがち、特に高次元になると
- そんなときに「全部のサブ木」について正確に計算してやって、尤度関数を最大にする特異点を見つけるアルゴリズムを回したりするとよいことがあるそうで、その「全部」を網羅して…というあたりが代数統計的(らしい)
GNJ法の背景的説明
- 木がある
- 木のエッジは、木の「内部にあるか、辺縁にあるか」について値を付与できる
  - そのエッジが木を分けてできる２つのサブ木の葉ノードの数の最小値を取ることで
- エッジにつけた「深さ」の尺度に従って、ある閾値より深い枝は省略する、というルールを作ると木がつぶれる
- つぶれた先には葉ノード数が一定数以下になったサブ木がぶら下がっている
- この葉ノード数が２と限らず、何個かであるという状態を扱う、というところがGNJ的
実際にGNJ法で解としての木に行きつくにはどうするか
- たとえば $\begin{pmatrix}k\\3\end{pmatrix}$ のときに代数統計アプリ(こちら)を使って、尤度関数のlocal maximaを求めるという手法によって解としての木が得られる(そうだ)
トリオに関して解いてみる、という例が挙げられている
- ３種のDNA配列比較をしている
- ３種には「共通祖先塩基」があると仮定し、種分化の過程で、共通祖先塩基を保存しているか、それから変化するかを考える。その変化のパターンは、塩基変化をすべて平等にすれば、４ｘ４の推移確率行列で表せて、その行列は対角成分と非対角成分とにそれぞれ異なる値を入れた単純なものになる
- 祖先から分化してからの「遠さ」は、時間で測ってもよいし、それと対応すると考えるのが通例だが、置換の多さで測ってもよい。(対角成分の値を決めると非対角成分の値も決まってしまう(推移確率は足して１の制約があるから)
- 「遠さ」によって、４ｘ４の値が違う
- また、３種のDNA配列について観測できることと言えば、「共通祖先塩基」は観察不明だから、３種で一緒、３種で別々、３種のうち２種で一緒(これに３通りある)の全部で５パターンに分けられる
- ３種で一緒の中身は、３種とも「祖先塩基」を守っている場合と、３種ともで「祖先塩基」から変異したけれど、たまたま同じ、という場合とで構成される
- 他のパターンの場合も同様にしてやると、４ｘ４の推移確率行列を構成している変数(１個の行列は１変数で確定するようにモデルを立ててあるから、計３変数)によって、観測データに関する尤度関数ができる
- 尤度を最大にする〜対数尤度を最大にする、という話に切り替わる
- 最大探しはlocal maxima探しなら、微分０探しとみなしてよいから、これを計算代数統計アプリでやりましょう、と
- MapleとSingular(こちら)を用いてやる、というそれのソースの一部が公開になっている(こちら)
- こちらはSolving the Likelihood Equationsと題された代数的に最尤推定する話→明日の記事へ続く
  - 中を覗くと、groebner()という関数やら、myring()という関数やら、idealと呼ばれるオブジェクト(？）やらが登場していて、いかにも計算機代数統計らしさがにじむ
  - ちなみに、処理は20,000観測ベクトル(データのこと)の処理に400MHz Pentium IIで50CPU hoursかかった、とのこと
こう呼んでくると、この論文は、実はそういう発想だったんだなー、と今さらに思う