系統解析のワークフロー(MrBayes) - ピエール瀧になりたい

Bayes法による微生物のたんぱく質のアミノ酸配列に対しての系統解析手法。

この方がいいという意見あれば教えて下さい。

なんとなく想定してるのは数百オーダーのアミノ酸配列の系統解析をする場合。

①MSA(multiple sequence alignment)データの前処理

最初はとても大事。NCBIから取ってきた配列は大体"|"とか" "とか、イヤなものがたくさん入っている。

MrBayesとかはこの辺の変な文字は結構嫌がるし、biopythonでFastaをパースするときも空白以下は読み込んでくれない。

一番いい解決策は、配列を全部 seq001, seq002 ,......, seqXXXみたいに名前を変更して、変更した名前のファイルを紐付けて保存しておくこと。これなら、最後系統樹を描画する直前に名前を戻してやればいいので楽。

②アラインメント

mafft。オプションは以下のようにいじるといいらしい。

③カラム抽出

trimal。データにもよるが、-gtオプションを0.8〜1.0くらいの間で振って、抽出される配列の長さを見てどれぐらい厳しくするかを決定する。

③モデル推定

AminosanでやるとMrBayesの入力を自動的に出力してくれる。DNAの場合は最適な置換モデル細かくを推定してくれるPartitionFinderの方がいい。でも使い方が割と面倒&計算に時間がかかる&スパコンに入れるのが大変(僕は結局スパコンで使うのは諦めました)。

④系統解析

置換モデルについてはAminosanの出力を利用。

その他のコマンドとしては、

mcmc Starttree=random ngen = 100000 printfreq = 1000 samplefreq = 100 nchains = 4

でmcmcを走らせ、Average standard deviation of split frequenciesの減少が見られなくなるまで続行。ngenはどれだけ時間かかるか予測つかない場合はちょっと少なめに設定して様子を見たりする。

sump burnin = xxx sumt

*ここでxxxは(最終的な世代数)/ Samplefreq数 / 4

事後確率が改善されていない(定常期に入っている)ことを確認。事後確率の散布図が出てくるから、それがランダムな分布をしていれば良し。

PSRFが1.0から離れている場合はサンプルの頻度を上げるか世代数(ngen)を増やす。

sumt conformat = simple Contype = Allcompat burnin = xxx

conformat = simpleにしないとR(ape)で樹形を弄れない。Contype = Allcompatを指定しない場合は、事後確率50%以下の分岐がコンセンサスになる。