幾霜::残日録::2004/11/05 (金)

　seqbootでjackknifeでリサンプリングしたデータセットを作って解析すれば現状でできるはずでは。BootPHYMLはseqbootでbootstrapで作成したデータセットをPHYMLに食わせてbootstrapやるんですが、この時seqbootに渡すオプションを書き換えてjackknifeでデータセットを書き出すようにしてしまえばlikelihood jackknifingになるはずです。でもこれってPAUP*でもできませんでした? 何か勘違いしているのかなぁ?

　1st,2nd,3rdコドンを分けて解析し、totalml(MOLPHY)とかで全体の系統樹を作成するような場合、bootstrapはどうするんや? PAUP*の場合はtotalml無しでそういう解析ができるけど、bootstrapやろうとするとその設定ではできませんと文句を言ってくる。えーと、まさか各コドンのデータを用意してそれぞれのbootstrapサンプリングしたデータからtotalmlでまとめて・・・をbootstrapの繰り返し分やるわけ? すげぇめんどいんですけど。スクリプト書けばいいっちゃーいいんですが、誰か書いてないのか? 発現領域ならよくやる解析行程やろ? うーむ、自分で書かないといけないんだろうか・・・。
　ところで、bootstrapの繰り返しのそれぞれでモデル選択ってしなくていいのか? これもよくわからん。っつーか各コドンのデータそれぞれでモデル選択はどうなんやろ。totalmlでまとめる時にはトポロジーと尤度さえあればええんやったっけ? ならモデルが違っても技術的には可能やけど・・・。
　で、モデル選択はともかくとして、bootstrapの手順も問題やなぁ。上記のように各コドンのデータからbootstrapサンプリングしてできた配列それぞれで系統解析してtotalmlでまとめて・・・とやる場合、totalmlでまとめるのが1000回で1000 replicateなのか、bootstrapサンプリングそれぞれ1000回で1000 replicateなのか。後者だとするなら結局1000^3=10億回totalmlでまとめないといけないんですけど。それはつらいなぁ。
　別の実装としては、元データから3塩基単位でbootstrapサンプリングしてから各コドンのデータに分ける。それぞれを解析して、totalmlでまとめる時は組になっていたデータから作られたものどうしをまとめる。これならまとめるのもbootstrapサンプリングも1000回で済むし、なんとなく組になっているものをちゃんと組として扱うことになるのでreasonableな気はする。しかし、問題はそういう風にbootstrapサンプリングするソフトが無いんとちゃうか、っちゅことやな。seqboot(PHYLIP)ってそんなことできへんだよなぁ。・・・と思ってseqbootのドキュメント見たら、あるじゃぁないですか、block-bootstrapping。ならスクリプトさえ書けば何とかなるかな。
　というわけでやはりモデル選択は各コドンで必要か、bootstrapのそれぞれで必要か、bootstrapサンプリングしたデータの各コドンで必要か、その辺が問題やなぁ。Inferring Phylogeniesにそんな話載っとったか? うーむ。コドンそれぞれでは必要だと思うし、bootstrapのそれぞれで必要なら、各コドン別にやる必要はあるでしょうけど、bootstrapのそれぞれで必要かどうかがよーわからん。そのモデルでリサンプリングデータから同じ結果になるかどうかを見るべきなのか、そのデータのベストなモデルで同じ結果になるかどうか見るべきなのか。なんとなく前者のような気はするが、他人に訊かれたら説明できんなぁ。むー。農環研の人にでも訊いてみようか。と書いたら反応あるんとちゃうかな(笑)。

追記 - 17:13:17
　totalmlの扱うllsファイルの形式が分からん。っつーか出力結果の意味も分からん。むぐ。

追記 - 18:30:47
　早速反応が。ありがたやありがたや。

ブーツストラップの精神はそのデータを“母集団”とみなして，そこからリサンプリングすることにあるとぼくは理解している
だから、母集団から選択されたモデルをリサンプリングデータは継承すべきだと。なるほど、確かにそうだ。ああーなんつーか、うまいこと言葉が出てこなかったんですが、良い感じのご説明ありがとうございます。

元データの各コドン位置ごとに分子進化モデルを選ぶのはリーズナブルですが．
ですよねぇ。でも既存の論文見てるとわざわざそこまでやっているのを見たこと無い気がするんですけど。で、簡単にそれができるソフトが無いんですよね。

ブーツストラップのやんちゃ精神は「元データをシバキたおしてばらつかせる」ことを目指しているので，ひょっとしたらパーティションを越えてブーツストラップするのが本来の野生の使い方なのかもしれまへんなあ．パーティションごとにブーツストラップするというようなお座敷風の行儀良さはあかんのんちゃう？
　うーん、そうでしょうか・・・。でもFelsensteinはそれが良いと思っているからblock-bootstrappingのオプションを付けたのでは。しかもドキュメントのその項では最後に
Note also that if you have a DNA sequence data set of an exon of a coding region, you can ensure that equal numbers of first, second, and third coding positions are sampled by using the block bootstrap with B = 3.
と書いてるんですよね。
　で、もしダメだとするなら、各コドン別にbootstrapする場合には上記のように各コドン毎にリサンプリングを行って解析する必要がありますが、問題は一つにまとめる時に組み合わせが膨大になるってことです。1000 replicateやると10億の組み合わせに。考えるだけでも頭が痛くなります。
　でも個人的にはやはりblock-bootstrappingに直感的魅力を感じますね。

　しかし何にしろtotalmlのやっていることと言うか、出力の解釈がよく分からん。これって異なるデータから得られた同じOTU組の複数のMLツリーと尤度を合わせることでデータ全体から推定されるツリーを得るものとちゃうんかなぁ? 最尤法の利点としてそれがよく挙げられていると思うんですけど、これ以外にそういうことができるというソフトを知らないんですが。結局llsファイルもよく分からん。困った。3rdコドンだけでやる? それはマズそうな。情報量減るし復帰置換の可能性が相対的に増えて遠いものでは問題が出てくるかも。

後ほどキンタ君から代理講義料を徴収したろかな．ほほー．
　あーどうぞどうぞ。K田先生は指導教官じゃないし(笑)。

追記 - 18:52:36
　あれ? PAUP*でBremer supportが計算できる? コマンドリファレンスには載ってないような気がするんですが・・・。

　もうMrBayesで突っ走ろっかなー。でもPAUP*と同様にコドン別に解析することはできますが、コドン毎に別のモデルを適用することはできないようで・・・。

追記 - 19:16:33
　ぐは。先ほどPAUP*の結果を見てみたら、MrBayesと全然違う。Lsetでsiterates=partition:by_codonなんてやってコドン別にやった結果。勿論chaesetとcharpartitionは事前に設定してある。これどういう計算しとんのかなぁ。あーソースコードほしぃー。あっても読めへんかもしれへんけど。

追記 - 06日00:46:45
　ただの1000回でいいですかそうですか。良かった(T^T)。というか確かにそう言われてみればそうな気がしてきた(ヲイ)。
　Bremer support indexに関しては最初はこれを読んでも意味が分からなかったんですが、今読み返したら分かりました(ヲイヲイ)。AutoDecayも今軽くソースコード見た感じではすぐ理解できそうだった(Perlで書いてくれてありがとう)のでこれで行こうかと思います。MacCladeは持ってないので・・・。研究室にあるかもしれませんが(笑)。多分あるけど行方不明でしょう。
　PAUP*でratesetってーのはcharpartitionを使うのと結局同じでは(コマンドリファレンスのp90のExample 3か4かの違いだけでやってることは同じではないかと)。でもこれってRates(Lsetのオプション)が各コドンで変わるだけで(=sitespecの場合)、Rmatrix、BaseFreq、Shape、Pinvarまで別個に設定できないですよね。それともそれは最尤法では技術的に無理? Ratesも各コドンに対してgamma補正できて欲しいんですが。勿論、各コドン毎に分けてモデル選択したらgamma補正やinvariant補正は必要無くなる可能性もありますが、それでもRmatrix、BaseFreq、Shapeがコドン毎に設定できないと別々のモデルを適用できていないのではないかと。で、このRates=sitespecの時、PAUP*はbootstrapできないんですね。そういうメッセージが出ます。そこからblock-bootstrappingが出てきているわけです。
　えっ、User's Manual無しでどうしてるのかって? ひたすらコマンドリファレンス読んでますが何か。理解するのに1ヶ月以上かかりましたがこれだけで何とかならなくはないです。この点に関してはPhylogenetic Trees Made Easyもcoding sequenceの場合はこうする、としてNEXUSスクリプトが書かれているのである程度参考になります。
　うちの組長については前に提出した超遅着レポートに書いてありますよ。

　メールが来ていたのでちょっと調べたらやはり話題になっているようですね。個人的にはこれでもうはてなは終わりだなという感じ。アンテナサービスも日記サービスも代替可能なサービスはいくらでもある。しかも今回の個人情報収集は不必要なものでもあるように思われます。不必要な個人情報を収集しようとするということ自体と、その場合にかかる管理コストを分かってないという無能さを露呈したことから、もう方針を撤回してもはてなは衰退に向かうと思います。はてなユーザーにはなまじそういうことに敏感な人間が多いので尚更。以前からはてなはどこか信用しきれないように感じていたので検討はしていましたが、そろそろはてなの利用はやめることにします。まぁ適当に欄を埋めて送信しておくというのも手ですけどね。

　先月は一ヶ月間一日も欠かさずここを更新し続けられたようです。a-Newsを使い始めてからは初めてですね。今月もこの調子で・・・と思ったらいきなり1日に何も書いてない・・・orz。