ミトコンドリアゲノムはあまり長くないので系統樹を確定させるには十分な情報を持っていないことがある.とくに非常に古い分岐では推定があいまいになることがある.したがって最尤系統樹以外の系統樹が実は真実である可能性を無視できない.ここでは,プログラムCONSELを用いて系統樹の確率値(p-value)を計算し,真実である可能性が無視できない系統樹をすべて列挙する.
1.さきほど実行したcmd371はwork-conc-user.llsというファイルを出力している.これは系統樹の対数尤度を分解してアミノ酸配列のひとつひとつの列にたいして対数尤度を計算したものであり,これをすべての系統樹に計算したサイズ105 x 3582の行列である.これを用いると,対数尤度のブートストラップ複製というものが計算できる.とくにマルチスケールブートストラップ法といって配列長を10通りかえて計算する方法をここでは用いている.(これらの手法の詳細はここでは説明しない.指示に従って計算すればよい.) コマンドは
cmd411 [enter]
である.結果はwork-conc-user.rmtというバイナリファイルに格納される.計算のログはwork-conc-rmt.logに書き出される.このコマンドの実行は時間がかかるので,バックグラウンドジョブとして
cmd411 & [enter]
として実行したほうが良いかもしれない.(あやまって何度も同じコマンドを実行しないこと! バックグラウンドジョブのことが良く理解できていないならばやらないこと.) この場合終了したかどうかはログファイルをみて最後に "exit normally"と出ていれば良い.
2.系統樹の確率値を計算するコマンドは
cmd431 [enter]
である.結果はwork-conc-tree.logに書かれる.確率値は次のように示されている.
# reading work-conc-tree.pv # 0 1 2 9+ 10 | 3 4 5 6 7 8 # rank item obs au np | bp pp kh sh wkh wsh | # 1 7 -12.7 0.944 0.559 | 0.560 1.000 0.803 1.000 0.803 1.000 | # 2 14 12.7 0.412 0.125 | 0.127 3e-06 0.197 0.879 0.197 0.801 | # 3 12 15.0 0.390 0.077 | 0.071 3e-07 0.154 0.854 0.154 0.732 | # 4 16 20.3 0.303 0.066 | 0.067 2e-09 0.129 0.816 0.129 0.714 | # 5 24 36.7 0.188 0.032 | 0.030 1e-16 0.077 0.534 0.077 0.559 | # 6 17 36.1 0.158 0.038 | 0.038 2e-16 0.081 0.549 0.081 0.566 | # 7 6 36.4 0.142 0.012 | 0.011 2e-16 0.048 0.536 0.048 0.416 | # 8 63 50.9 0.141 0.025 | 0.026 8e-23 0.061 0.314 0.061 0.465 | # 9 1 29.9 0.128 0.014 | 0.013 1e-13 0.036 0.639 0.036 0.345 | # 10 20 40.4 0.118 0.012 | 0.012 3e-18 0.052 0.473 0.052 0.469 | # 11 5 38.8 0.117 0.008 | 0.009 1e-17 0.042 0.490 0.042 0.382 | # 12 62 58.7 0.116 0.011 | 0.010 3e-26 0.032 0.209 0.032 0.366 | # 13 10 37.9 0.085 0.007 | 0.008 3e-17 0.037 0.510 0.037 0.344 | # 14 29 46.5 0.074 0.005 | 0.006 6e-21 0.038 0.376 0.038 0.362 | # 15 9 41.6 0.072 0.004 | 0.003 8e-19 0.026 0.447 0.026 0.283 | # 16 35 61.7 0.057 0.004 | 0.004 2e-27 0.025 0.179 0.025 0.287 |
確率値には色々な定義があるがここでは"au"と書かれた列を見る.この値が0.05以上ならば有意水準5%で棄却できないことをあらわす.ここでは105個中の上位16個が棄却できない.ただしauの大きい順に系統樹が並べられている.work-conc-user.tpl の番号で言うと7, 14, 12, 16, 24 17, 6, 63, 1, 20, 5, 62, 10, 29, 9, 35番めの系統樹である.
3.枝にも確率値を計算することができる.コマンドは
cmd451 [enter]
である.結果はwork-conc-edge.logに書かれる.確率値は次のように示される.
# reading work-conc-edge.pv # 0 1 2 9+ 10 | 3 4 5 6 7 8 # rank item obs au np | bp pp kh sh wkh wsh | # 1 5 -29.9 0.959 0.845 | 0.845 1.000 0.964 0.999 0.919 1.000 | # 2 19 -20.3 0.929 0.805 | 0.803 1.000 0.871 0.995 0.871 0.998 | # 3 1 -12.7 0.856 0.710 | 0.710 1.000 0.803 0.990 0.803 0.996 | # 4 6 12.7 0.251 0.159 | 0.161 3e-06 0.197 0.847 0.197 0.734 | # 5 9 15.0 0.178 0.130 | 0.127 3e-07 0.154 0.823 0.154 0.661 | # 6 20 36.7 0.167 0.032 | 0.030 1e-16 0.077 0.529 0.077 0.525 | # 7 7 20.3 0.165 0.137 | 0.137 2e-09 0.129 0.790 0.129 0.651 | # 8 17 36.1 0.153 0.038 | 0.038 2e-16 0.081 0.546 0.081 0.535 | # 9 16 58.7 0.106 0.011 | 0.010 3e-26 0.032 0.209 0.032 0.345 | # 10 12 50.9 0.098 0.027 | 0.028 8e-23 0.061 0.309 0.061 0.417 | # 11 3 29.9 0.094 0.035 | 0.034 1e-13 0.036 0.626 0.048 0.373 | # 12 14 40.4 0.082 0.022 | 0.022 3e-18 0.052 0.461 0.052 0.421 | # 13 23 46.5 0.071 0.035 | 0.038 6e-21 0.038 0.372 0.061 0.420 |
のようにauが0.05以上になるのは13個の枝であり,番号5, 19, 1, ..., 23である.ただし,この番号がどの枝を意味するかの対応表はwork-conc-ass.logに書かれている.ファイルの中ごろでまず葉の番号を定義してから枝のリストが続く.
# leaves: 25 25 1 NC_001325 2 NC_001602 3 NC_004023 4 NC_004029 5 NC_003428 6 NC_002008 7 NC_001700 8 NC_001640 9 NC_001808 10 NC_000845 11 NC_001567 12 NC_000889 13 NC_001601 14 NC_002503 15 NC_001913 16 NC_001643 17 NC_001644 18 NC_001807 19 NC_001645 20 NC_002083 21 NC_004025 22 NC_002658 23 NC_001569 24 NC_001610 25 NC_000891 # base edges: 25 25 25 1 2 1234567890123456789012345 1 ++++++++++++++++++++++--- ; (0.1429) 2 ++++++++++++++--------+-- ; (0.1429) 3 +++++++++++++++---------- ; (0.1429) 4 ++++++++++++++-++++++++-- ; (0.1429) 5 ++++++++++++++-++++++---- ; (0.1429) 6 +++++++++++++++++++++-+-- ; (0.1429) 7 --------------+------+--- ; (0.1429) 8 --------------+++++++---- ; (0.1429) 9 ---------------------++-- ; (0.1429) 10 ---------------+++++++--- ; (0.1429) 11 ++++++++++++++-------+--- ; (0.1429) 12 +++++++++++++++------++-- ; (0.1429) 13 --------------+++++++++-- ; (0.1429) 14 --------------+-------+-- ; (0.1429) 15 ---------------++++++-+-- ; (0.1429) 16 +++++++++++++++-------+-- ; (0.0857) 17 --------------++++++++--- ; (0.0857) 18 ++++++++++++++-+++++++--- ; (0.0857) 19 +++++++++++++++++++++---- ; (0.0857) 20 +++++++++++++++------+--- ; (0.0857) 21 ++++++++++++++-------++-- ; (0.0857) 22 ---------------++++++++-- ; (0.0857) 23 --------------+------++-- ; (0.0857) 24 --------------+++++++-+-- ; (0.0857) 25 ++++++++++++++-++++++-+-- ; (0.0857)
のように,+と-マークの組み合わせで枝が表現されている.従ってau=0.959でrank=1の枝はitem=5だから
5 ++++++++++++++-++++++---- ;
と表現され,これは(NC_001325,...,NC_002503,NC_001643,...,NC_004025)の生物群についた枝である.これ以降に書かれている枝のリスト
# common edges: 44 44 25 1 2 1234567890123456789012345 26 +++++++++++++++++++++++-- ; (1.0000) 27 ++++++++++++++----------- ; (1.0000) 28 +++++++++---------------- ; (1.0000) 29 +++++++------------------ ; (1.0000) 30 ++++++------------------- ; (1.0000) 31 +++++-------------------- ; (1.0000) 32 ++++--------------------- ; (1.0000) 33 ++----------------------- ; (1.0000) 34 +------------------------ ; (1.0000) 35 -+----------------------- ; (1.0000) 36 --++--------------------- ; (1.0000) 37 --+---------------------- ; (1.0000) 38 ---+--------------------- ; (1.0000) 39 ----+-------------------- ; (1.0000) 40 -----+------------------- ; (1.0000) 41 ------+------------------ ; (1.0000) 42 -------++---------------- ; (1.0000) 43 -------+----------------- ; (1.0000) 44 --------+---------------- ; (1.0000) 45 ---------+++++----------- ; (1.0000) 46 ---------+--------------- ; (1.0000) 47 ----------++++----------- ; (1.0000) 48 ----------+-------------- ; (1.0000) 49 -----------+++----------- ; (1.0000) 50 -----------+------------- ; (1.0000) 51 ------------++----------- ; (1.0000) 52 ------------+------------ ; (1.0000) 53 -------------+----------- ; (1.0000) 54 --------------+---------- ; (1.0000) 55 ---------------++++++---- ; (1.0000) 56 ---------------+++++----- ; (1.0000) 57 ---------------++++------ ; (1.0000) 58 ---------------+++------- ; (1.0000) 59 ---------------++-------- ; (1.0000) 60 ---------------+--------- ; (1.0000) 61 ----------------+-------- ; (1.0000) 62 -----------------+------- ; (1.0000) 63 ------------------+------ ; (1.0000) 64 -------------------+----- ; (1.0000) 65 --------------------+---- ; (1.0000) 66 ---------------------+--- ; (1.0000) 67 ----------------------+-- ; (1.0000) 68 -----------------------+- ; (1.0000) 69 ------------------------+ ; (1.0000)
は105通りの候補系統樹すべてに含まれていたので確率値=1.00である.
4. cmd431の出力はなんとも読みにくい.そこでwork-conc-clade.txtに従って生物群にラベルを付けて結果を整理する.まずwork-conc-user.tplもしくはwork-conc-ass.logから系統樹リストの一覧の部分だけを切り取ってwork-conc-tree-list1.txtを作る.そして
repnamegb -sv -p 4 -l work-conc-clade.txt work-conc-tree-list1.txt > work-conc-tree-list2.txt
と入力してwork-conc-tree-list2.txtを書き出す.そして
sortlines -c 3 -f work-conc-tree.log work-conc-tree-list2.txt > work-conc-tree-sort.txt
と入力するとwork-conc-tree-sort.txtが書き出される.これは系統樹のラベルをwork-conc-clade.txtで定義したものに置き換えた上でauの大きい順に系統樹を並び替えたものである.この上位16個
(((((A,B),C),D),E),F); 7 <- 7 (0.0095) (((((A,B),C),E),D),F); 14 <- 14 (0.0095) ((((A,B),C),(D,E)),F); 12 <- 12 (0.0095) ((((A,B),(C,D)),E),F); 16 <- 16 (0.0095) ((((A,(C,D)),B),E),F); 24 <- 24 (0.0095) (((A,(B,(C,D))),E),F); 17 <- 17 (0.0095) ((((A,C),B),(D,E)),F); 6 <- 6 (0.0095) (((A,(C,(D,E))),B),F); 63 <- 63 (0.0095) (((((A,C),B),D),E),F); 1 <- 1 (0.0095) ((((A,B),(C,E)),D),F); 20 <- 20 (0.0095) (((((A,C),B),E),D),F); 5 <- 5 (0.0095) ((((A,(C,E)),B),D),F); 62 <- 62 (0.0095) (((A,(B,C)),(D,E)),F); 10 <- 10 (0.0095) (((A,B),(C,(D,E))),F); 29 <- 29 (0.0095) ((((A,(B,C)),E),D),F); 9 <- 9 (0.0095) ((A,(B,(C,(D,E)))),F); 35 <- 35 (0.0095)
が系統樹の信頼集合である.
4.cmd451の出力も読みにくい.そこでwork-conc-clade.txtに従って生物群にラベルを付けて結果を整理する.まずwork-conc-ass.logから枝リストの一覧の部分だけを切り取ってwork-conc-edge-list1.txtを作る.そして今度は手作業でラベル付けを行いwork-conc-edge-list2.txtを作成する.#のついた行はコメントなので無視される.そして,
sortlines -c 3 -f work-conc-edge.log work-conc-edge-list2.txt > work-conc-edge-sort.txt
で結果をソートしてwork-conc-edge-sort.txtを作成する.この上位13個の枝をみると
5 ++++++++++++++-++++++---- ; (A,B) 19 +++++++++++++++++++++---- ; (A,B,C) 1 ++++++++++++++++++++++--- ; (A,B,C,D) 6 +++++++++++++++++++++-+-- ; (A,B,C,E) 9 ---------------------++-- ; (D,E) 20 +++++++++++++++------+--- ; (A,C,D) 7 --------------+------+--- ; (C,D) 17 --------------++++++++--- ; (B,C,D) 16 +++++++++++++++-------+-- ; (A,C,E) 12 +++++++++++++++------++-- ; (A,C,D,E) 3 +++++++++++++++---------- ; (A,C) 14 --------------+-------+-- ; (C,E) 23 --------------+------++-- ; (C,D,E)となっていてこれらが枝の信頼集合である.
5.TreeViewという系統樹をグラフィックで表示するソフトを使ってwork-conc-user.phyの系統樹のいくつかを表示してみよう(ただしeduにはtreeviewはインストールされていないのでこれはやらなくてもよい).ここでwork-conc-user.phyはwork-conc-user.treをauの大きい順に並べ替えたものである.
AUを一番大きくするのはitem=7である.この系統樹はwork-conc-sort-root-1.pdf(もしくは根無しで書くとwork-conc-sort-unroot-1.pdf)である.この系統樹が真実であるとすると,これまでの生物学の知識に反する結果であるので,新たな発見をしたことになる.つまり,これまではウサギとネズミは近縁で一つのグループを作ると考えられていたからである.ところがAUの順位で14番目(根有り,根無し)と16番目(根有り,根無し)の系統樹をみると,ここでは古典的な結果が示唆されていて,これらの系統樹もAUが0.05以上になっているので棄却されない.従って,新たな発見をしたとは必ずしもいえないわけだ.
自分で分析した結果に対しても同様の考察をしてみよ.