幾霜::残日録::2008/06/05 (木)

　ぶーとすとらっぱースクリプトをNEXUS用とPHYLIP用で別々にしていたのを、入力ファイル形式判定ルーチンを付けることで統合した。別々だと、メンテナンスコストがかかってしょうがない。

　募集開始してたっけ。会場は東大駒場だっけな。宿泊費が高くつきそうな気がする。なるべく出たいのだが、どうしようか。懐具合と相談かな。えーと、渋谷駅から1kmね。渋谷付近のホテルから歩けば通えるか。学会会場まで電車で通うなんて絶対お断りやから早めに決めなあかんな。

　これも考えてみたらさほど難しいことはなかった。まず指定された外群のOTUを全て含むクレードの中で最も浅いクレードを探す。そのクレード直下の枝とクレード直上の枝が新しい有根系統樹の根点直下の枝となる。クレード直下の枝とそこに繋がっているクレードを出力する有根系統樹に加えたら、逆にクレード直上の枝を遡っていきつつ続けて加えていく。入力系統樹の根点まで到達したら、たどってきたクレード以外を加えて完成。おしまい。

　まぁ実際に実装するのはけっこう面倒なんだが、Treefinderの系統樹読み込みと処理って結構遅いので自前で実装した方が良いかもしれない。

　よそ様に頼ってばかりではいかんと思って考えてみたらさほど難しいことはなかった。

　全ての枝長を1と仮定して全OTU間の樹上の距離行列を得る。その距離行列が一致するなら必ず樹形も同じであり、一致しないなら樹形も必ず異なる。

　また、Symmetric Differenceは距離行列の各マス目の差の絶対値の和を2で割ればよい。ただし、OTU1-OTU2の距離とOTU2-OTU1の距離を両方使った場合は4で割る。最小NNIステップ数はSymmetric Differnceの半分である。

　・・・であっとんのかいな。Semple & Steel (2003)にはきっと書いてあるんやろうけど、難しすぎて死ぬ。簡単なことも書き方が難しくて理解不能なんよねこの本。

　一旦自力で樹形をパースして各枝の番号を得る。得る番号は、PHYLIPのretreeが各ノードに割り当てるのと同じ値になるようにする。具体的には、まず各OTUが1からの値になり、次に最上位のクレードから順に浅いクレードへと番号を振る。指定された外群OTUを全て含む最も浅いクレードの番号を探し、retreeにて有根化する。これを各樹形に対して行えばよい。

　が、こんな面倒なことをしなくともTreefinderで一発であった。MesquiteでもできるがGUIで手作業なので私の好みではない。

　TreefinderにはRerootという関数があるがこれではなく、RedrawTreeを使う。以下のようにコマンドを実行する。

SaveTreeList[
 RedrawTree[
   LoadTreeList["treefile"],
   Outgroup->{"OTU1","OTU2"}
 ],
 "outfile",
 Format->"NEWICK"
]

これにてtreefile内の全系統樹が有根化されてoutfileへ保存される。外群OTUはもっと多数指定しても構わない。

　PAUP*があれば簡単。と言ってもPAUP*は樹形を読むにはその樹形のOTUの形質データが必要なので事前に配列データを読み込ませておく必要があるのが難点ではある。とにかくデータを読み込ませてから以下のコマンドを実行すればよい。

GetTrees File=treefile DupTrees=Eliminate

これで重複は除去された形でメモリ上にロードされるので、SaveTreesすればよい。

　他にもPHYLIPのtreedistを使う方法もある。PHYLIP形式樹形ファイルを与えて全組み合わせの距離を書き出すようにしてやれば樹形間の距離行列が得られるので距離ゼロの組み合わせを除去すればいい・・・が、かなりめんどうくさい。なお、Branch Score Distanceではなく、Symmetric Distanceを出力するように設定する必要がある。

　Treefinderでも一応できなくはない。RobinsonFoulds関数を使う。以下のように用いる。

RobinsonFoulds[LoadTree["treefile1"],LoadTree["treefile2"]]

これで二つの樹形間の距離は得られる。やっぱり距離ゼロの組み合わせを探して除去しなくてはならないのは同じ。しかも1樹形1ファイルにして全対比較しなくてはならないのでtreedistより面倒くさい。

　treedistもTreefinderも、枝長が8.967e-008などという指数表記だとエラーになるようです。枝長を除去してやった方が良さそうですね。

　枝長を除去したファイルを用意してtreedistにて樹形間の距離行列を得、距離ゼロ、即ち同一の樹形を最初のものを残して除去する。といったところでしょうか。まぁそんな難しい処理ではないか。

　PARのppでカレントディレクトリの.plファイルを全て.exeにする。

for %I in (*.pl) do pp -o %~nI.exe %I

%~nIは変数%Iに格納されているファイル名の拡張子を除いた狭義のファイル名部分を取り出す。その他にも色々あるのでfor /?を参照。

　Mersenne TwisterのPure Perl実装モジュール。可搬性に関しては一番安心できそうです。速度があまりに遅くないのならこれを同梱しましょうかね。

追記 - 03:31:02
　Math::Random::MTより若干遅いがものすごい差というわけでもないもよう。ActivePerl内蔵の純正rand()と比べるとかなり遅い。それはMath::Random::MTでも同様。