RAxML via Condor

2007 年 10 月 12 日　改訂
井上 潤


　このページは condor という Supercomputer の管理システムを用いて RAxML の解析を行うための個人的なメモです．
　RAxML の一般的な操作方法に関してはこちらをご覧下さい．
最尤樹推定

[inoue@petal 12tr_53Noto_RM_Fol]$ cat NondRM_Noto53.cmd
initialdir = /home/inoue/53Noto_Fol/RAxML_Noto53/12tr_53Noto_RM_Fol
Rank = kflops

Executable = /usr/common/x86_64-linux/bin/raxmlHPC
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = Noto53_12tr.PHYLIP,partition_12tr

notification = NEVER
arguments = -f d -m GTRMIX -s Noto53_12tr.PHYLIP -q partition_12tr -n NondMultipleOriginal_$(Process)
output = NondMulOri.out$(Process)
error = NondMulOri.error$(Process)
log = NondMulOri.log$(Process)

Queue 100
grep GAMMA-likelihood RAxML_info.MultipleOriginal_*
スクリーンに出てきた表示をテキストファイルとして保存し，Excel で読み込んで尤度で並べ替える．
アミノ酸解析の場合 (PROTMIXMTREVF)，-# 100 は，petal でも phoenix でも途中でとまってしまう．塩基はわからない．
最適な initial rearrangement (-i) とカテゴリー数 (-c) の設定を探す．

最節約法を用いた starting tree の推定までは，condor を用いないで行う．

実際に解析を行った感触では，-i と -c は以下の操作の限りでは，機能しているかどうか不明であった．というのも，-i と -c ではなく，starting tree (ここでは MP tree) にのみ依存して尤度が決まっていると思えるため．つまり，同じ starting tree を用いた場合では，-i Auto と -i 10 あるいは，-c 10 と -c 40 の間で尤度に変化が無かった．


とりあえず，上記の操作を一通り記す．
まず，initial rearrangement の最適な設定を探す．
10 として (-i 10)，5 通りの最尤樹を求める．
[inoue@petal 12tr_53Noto_RM_Fol]$ cat i10RM_Noto53.cmd
initialdir = /home/inoue/53Noto_Fol/RAxML_Noto53/12tr_53Noto_RM_Fol
Rank = kflops

Executable = /usr/common/x86_64-linux/bin/raxmlHPC
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = Noto53_12tr.PHYLIP,partition_12tr,RAxML_parsimonyTree.ST$(Process)

notification = NEVER
arguments = -f d -i 10 -m GTRMIX -s Noto53_12tr.PHYLIP -q partition_12tr -t RAxML_parsimonyTree.ST$(Process) -n i10_$(Process)
output = i10RM.out$(Process)
error = i10RM.error$(Process)
log = i10RM.log$(Process)

Queue 5
rearrangement setting が自動的に設定されるようにして (-i オプションを外す)，5 通りの最尤樹を求める．
[inoue@petal 12tr_53Noto_RM_Fol]$ cat iAutoRM_Noto53.cmd
initialdir = /home/inoue/53Noto_Fol/RAxML_Noto53/12tr_53Noto_RM_Fol
Rank = kflops

Executable = /usr/common/x86_64-linux/bin/raxmlHPC
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = Noto53_12tr.PHYLIP,partition_12tr,RAxML_parsimonyTree.ST$(Process)

notification = NEVER
arguments = -f d -m GTRMIX -s Noto53_12tr.PHYLIP -q partition_12tr -t RAxML_parsimonyTree.ST$(Process) -n iAuto_$(Process)
output = iAutoRM.out$(Process)
error = iAutoRM.error$(Process)
log = iAutoRM.log$(Process)

Queue 5
最適なカテゴリー数を探す．
カテゴリーを 10 (-c 10) として解析する
[inoue@petal 12tr_53Noto_RM_Fol]$ cat c10RM_Noto53.cmd
initialdir = /home/inoue/53Noto_Fol/RAxML_Noto53/12tr_53Noto_RM_Fol
Rank = kflops

Executable = /usr/common/x86_64-linux/bin/raxmlHPC
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = Noto53_12tr.PHYLIP,partition_12tr,RAxML_parsimonyTree.ST$(Process)

notification = NEVER
arguments = -f d -c 10 -m GTRMIX -s Noto53_12tr.PHYLIP -q partition_12tr -t RAxML_parsimonyTree.ST$(Process) -n c10_$(Process)
output = c10RM.out$(Process)
error = c10RM.error$(Process)
log = c10RM.log$(Process)

Queue 5

カテゴリーを 40 (-c 40) として解析する
[inoue@petal 12tr_53Noto_RM_Fol]$ cat c10RM_Noto53.cmd
initialdir = /home/inoue/53Noto_Fol/RAxML_Noto53/12tr_53Noto_RM_Fol
Rank = kflops

Executable = /usr/common/x86_64-linux/bin/raxmlHPC
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = Noto53_12tr.PHYLIP,partition_12tr,RAxML_parsimonyTree.ST$(Process)

notification = NEVER
arguments = -f d -c 10 -m GTRMIX -s Noto53_12tr.PHYLIP -q partition_12tr -t RAxML_parsimonyTree.ST$(Process) -n c10_$(Process)
output = c10RM.out$(Process)
error = c10RM.error$(Process)
log = c10RM.log$(Process)

Queue 5
[inoue@petal 12tr_53Noto_RM_Fol]$ cat c40RM_Noto53.cmd
initialdir = /home/inoue/53Noto_Fol/RAxML_Noto53/12tr_53Noto_RM_Fol
Rank = kflops

Executable = /usr/common/x86_64-linux/bin/raxmlHPC
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = Noto53_12tr.PHYLIP,partition_12tr,RAxML_parsimonyTree.ST$(Process)

notification = NEVER
arguments = -f d -c 40 -m GTRMIX -s Noto53_12tr.PHYLIP -q partition_12tr -t RAxML_parsimonyTree.ST$(Process) -n c40_$(Process)
output = c40RM.out$(Process)
error = c40RM.error$(Process)
log = c40RM.log$(Process)

Queue 5
grep で尤度を比較する．
grep Likelihood RAxML_info.*


Constraint 付きの解析

[inoue@phoenix RM_ML_12tr_Const_40Noto_fol]$ pwd
/home/inoue/b40Noto_fol/ConstAnal_40Noto_fol/
12tr_Const_b40Noto_fol/RM_ML_12tr_Const_40Noto_fol
ブートストラップ解析


ランダムシードナンバーを作成する script

RMBS_com.pl
通常の Bootstrap
[inoue@phoenix RM_BS_123ACtr_40Noto_fol]$ cat RMBS.txt
initialdir = /home/inoue/b40Noto_fol/AC_123trAC_40Noto_fol/RM_BS_123ACtr_40Noto_fol
Rank = kflops
[inoue@petal 12tr_53Noto_RM_Fol]$ cat RMBS.cmd
initialdir = /home/inoue/53Noto_Fol/RAxML_Noto53/12tr_53Noto_RM_Fol
Rank = kflops

Executable = /usr/common/x86_64-linux/bin/raxmlHPC
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = Noto53_12tr.PHYLIP

notification = NEVER
output = RMBS.$(Process)
error = RMBS.$(Process)
log = RMBS.log

arguments =-m GTRCAT -s Noto53_12tr.PHYLIP -b 418260808 -n BOOT0
Queue

arguments =-m GTRCAT -s Noto53_12tr.PHYLIP -b 286988123 -n BOOT1
Queue

 .............

arguments =-m GTRCAT -s Noto53_12tr.PHYLIP -b 935017833 -n BOOT99
Queue

tree ファイルの作成
#nexus
begin trees;

を先頭に

tree mytree =

各行の最初に挿入??て，ファイルの最後に

end;

を加える．
* vi エディターで検索置換を行う場合のコマンド．

1,$s/\n/\r\rtree mytree = /g
paup で consensus tree を行う control file．
ガンマ補正なしの Bootstrap．
infile と同じフォルダにコンパイルしたプログラムを置いておく．
[inoue@petal Ex1stCyt57Act_RMBS]$ pwd
/home/inoue/CybMem_RABS_Act57/Ex1stCyt57Act_RMBS
[inoue@petal Ex1stCyt57Act_RMBS]$ cat Ex1stCybNonG_RMBS.cmd
initialdir = /home/inoue/CybMem_RABS_Act57/Ex1stCyt57Act_RMBS
Rank = kflops

Executable = ./NonG_raxmlHPC
[ここでは gamma 補正を外すために，別途コンパイルしたプログラムを用いている]
Universe = vanilla
requirements = (OpSys =="LINUX" && Arch =="X86_64")

should_transfer_files = YES
when_to_transfer_output = ON_EXIT_OR_EVICT
transfer_input_files = ExCyb1st.PHYLIP

notification = NEVER
output = RMBPCp.$(Process)
error = RMBPCp.$(Process)
log = RMBPCp.log

arguments =-m GTRCAT -s ExCyb1st.PHYLIP -b 418260808 -n BOOT0
Queue

arguments =-m GTRCAT -s ExCyb1st.PHYLIP -b 286988123 -n BOOT1
Queue

.............

arguments =-m GTRCAT -s ExCyb1st.PHYLIP -b 857737186 -n BOOT99
Queue
ソースコードを変更して RAxML を上書き可能にする

Condor は解析の途中で node (?) を変更するため，スイッチしたときにどうも RAxML が outputfile を見つけて解析がストップしてしまいます．60 OTU 程度では問題にならないのですが，120 OTU を超えた解析を行うと，RAxML がストップすることが頻繁に出てきました．

この方法は作者に教えて頂いたのですが，dangerous なのであまりお薦めできない，と言っていました．やるなら名前を変更するにようにコードを書きかえてはと進められましたが，私はそのやり方がわかりません．

The check is in axml.c, function

static void makeFileNames(tree *tr, analdef *adef, int argc)
if(processID == 0)
{
infoFileExists = filexists(infoFileName);
..........
}

if you comment out that one you should be able to overwrite files.

コメントに従って，以下の選択部分を削除してコンパイルしました．いまのところ問題なく解析が動いているようです．得られる結果も問題ないように思えます．
AA rate matrix を書きかえる

AA matrix を Excel ファイルで作成し，スペース区切りのテキストで保存する．BBEDIT で加工して model.c の該当部分を書きかえる．アミノ酸の並び方は PAML と同じ．matrix の要素は，小数点何桁でも良いらしい．それと mtmama にも 0 が含まれているので，0 があっても問題ないようだ．ただし，Pi は 0 があるとまずいと思う．例題
脊椎動物 781 spp の (Nexus) OTU 名を書きかえる perl スクリプト

perl -p ExcludeMam.txt infile > outfile

で実行する．

Vertebrate 781 spp のデータセット (nexus 形式) から mammals 67 spp を抜き出すスクリプト
Vertebrate 781 spp のデータセット (nexus 形式) の OTU 名を，RAxML で読み込み可能なパターンに置き換えるスクリプト