ゲノム情報をもとに数百の遺伝子を使って系統解析をしようとしたのですが、ファイルの結合などのいいツールが見つからなかったので、Biopythonでそこいらの処理が楽になるようにしました。またConselをつかったAU検定でちょっとつまずいたので、そこについても記録しておきます。

初心者なので、間違いなどのコメント歓迎です。

参考サイト

https://www.fifthdimension.jp/documents/molphytextbook/
参考サイトが十分わかり安いと思うので、しっかりやりたい人はそっちを見て下さい。
Concelの部分だけ、私の環境ではうまく動かなかったので、変更してます。

使用ソフトウエア

RAxML (なにはなくとも）
http://sco.h-its.org/exelixis/software.html
Phylogears (treファイルの編集などにつかう）
https://www.fifthdimension.jp/products/phylogears/
　(Phylogears2にしかないツールを使うので、αリリースだけどそっちをいれる）
FigTree (系統樹ファイルを開くのにつかう）
http://tree.bio.ed.ac.uk/software/figtree/
Consel (AU検定などの検定を一挙にやってくれる）
http://www.sigmath.es.osaka-u.ac.jp/shimo-lab/prog/consel/

よく感謝して引用しよう。

AU検定による信頼性検定

所与の配列データから、ある系統関係が有意に棄却されるかを検定するのだ。
実際は、既に系統樹が得られているときに、先行研究と矛盾したとき等に、自分の系統樹の信頼性をブートストラップ値以外で検証したい時に使うことが多い。けどこれもp値の問題を抱えているのでそのうち信頼区間を出せとかいわれんのかなー？

1. 帰無仮説となる系統樹の作成

概要

まず、帰無仮説となる系統樹を作る。信頼性を検証したい系統仮説（ノード）と並立しない（非互換な）系統仮説を想定して系統樹をつくればよい。

すなわち
(A,(B,(C,D)))
があり
(C,D)
の部分が疑わしいならば、
(B,C)
を固定するなどして、帰無系統樹を作れば良い。

しかし、すべての非互換系統仮説を検討していると大変なので、最初にあり得そうなものを選ぶ。その後、その非互換系統仮説部分を固定(制約）して、RaxMLなどで系統樹を作成する。つまり、制約した部分以外は、普通に系統解析をして、全体の樹形を求め、それに対して、検定をおこなう。

樹形の制約条件は、既出の論文を元に自分で選出しても良いだろうが、
参考サイトの手法では以下の手順でBS解析で出てきた非互換な系統仮説を制約としてつかっている。

普通にRaxMLをする(-f aをつかって最尤系統樹とBS値を同時に計算する）
信頼性を検証したい系統仮説（ノード）を取り出す
BS解析で作られた複数のツリーから、検証したい系統仮説とは非互換な系統仮説を持つ物を取り出し、最頻出のツリーを選出する
得られた非互換最頻出仮説を制約として、もう一度RaxML解析をする
得られたbest tree が第一にチェックしとくべき「帰無仮説となる系統樹」である。

これによって、BS解析で少数だが支持された「あり得るかもしれないが非互換な系統仮説」を帰無仮説として検定を行える。一方、矛盾する先行研究が形態データなどを元にしている場合は自分でその系統仮説をnewick形式でつくり、それを制約としてRAxMLをすればよい。

以下詳細

手順１，普通にRaxMLをする

First_RaxML.sh

# !/bin/sh 
# $ -S /bin/bash
# $ -cwd
# $ -v PATH
raxmlHPC-PTHREADS-SSE3 -f a -m PROTGAMMAWAG -p 12345 -x 12345 -# 100 -s $INPUT_FASTA -n $OUT_PREFIX -T 20

手順２，信頼性を検証したい系統仮説（ノード）を取り出す

RAxML_bootstrap.$OUT_PREFIXファイルにはRaxMLで作られたブートストラップ解析結果すべてが入っている。そのなかから各ノード（系統仮説）をpgsumtree(Phylogears2中のツール）をつかってバラバラにする。

all_node_extract.sh

pgsumtree \
--mode=ALL \ 
RAxML_bootstrap.$OUT_PREFIX \ 
RAxML_allhypotheses.tre

RAxML_allhypotheses.tre をfigtreeなどの描画ソフトで開いて、注目しているノードが入ったツリーを選ぶ。見つかったらpgsplicetreeをつかっってそのツリーを別ファイルに切り出す。

tree_extract.sh

pgsplicetree \
5 \ #取り出したいツリーが何番目のツリーかを指定する 
RAxML_allhypotheses.tre \ 　
RAxML_allhypotheses.5.tre

手順３，並立しない系統仮説をBS解析の全ツリーから探し頻出のものを取り出す。

選んだ系統仮説と非互換の仮説の中でBS解析によって最も支持されたものを見つけ出す。

Unmatch_node_extract.sh

pgsumtree \
--mode=MAJi \ #非互換最頻出仮説探索モード
--treefile=RAxML_allhypotheses.5.tre \ #この系統仮説と非互換のものがほしい！
RAxML_bootstrap.$OUT_PREFIX \ 
RAxML_MAJi_hypotheses.5.tre \　#出力ファイル名

出てきたRAxML_MAJi_hypotheses.5.treが非互換最頻出仮説である。

手順4，非互換最頻出仮説を制約として、もう一度RaxML解析をする

limited_RaxML.sh

raxmlHPC-PTHREADS-SSE3 -f a -m PROTGAMMAWAG -p 12345 -x 12345 -# 100 -s $INPUT_FASTA -n   MAJi_hypotheses.5.tre -T 20 -g RAxML_MAJi_hypothesis5.tre

上記の解析ででてきた
RAxML_bestTree.MAJi_hypotheses.5.tre
が帰無仮説系統樹となる

２. 座位毎の尤度計算

最尤系統樹と帰無仮説系統樹のそれぞれについて座位毎の尤度計算を行う。

手順１，最尤系統樹と帰無仮説系統樹を一つのファイルにする

trees_joint.sh

pgjointree \
RAxML_bestTree.MAJi_hypotheses.5.tre \
RAxML_bestTree.$OUT_PREFIX.tre \
RAxML_forAUtest.tre \

手順２，座位毎の尤度計算を行う

eachsite_likely_RaxML.sh

raxmlHPC-PTHREADS-SSE3 -f G -m PROTGAMMAWAG -s $INPUT_FASTA -n eachsite.outfile -T 20 -z RAxML_forAUtest.tre

RAxML_perSiteLLs.eachsite.outfileが結果ファイル

３. ConselによるAU検定

以下のコマンドでp値のリストがでる。

Consel.sh

seqmt --puzzle RAxML_perSiteLLs.eachsite.outfile #"makermt --puzzle"でいっぺんに処理せずに変換をかます
makermt RAxML_perSiteLLs.eachsite
consel RAxML_perSiteLLs
catpv RAxML_perSiteLLs > consel.out

例えば

consel.out

# reading RAxML_perSiteLLs.pv
# rank item    obs     au     np |     bp     pp     kh     sh    wkh    wsh |
#    1    1  -67.8  0.996  0.995 |  0.994  1.000  0.994  0.994  0.994  0.994 |
#    2    2   67.8  0.004  0.005 |  0.006  4e-30  0.006  0.006  0.006  0.006 |