NGSデータ解析まとめ

サカナ研究者の手探りNGS解析(おもに進化生物学)

Shell scriptでも試してみる

複数FASTA fileのinputから、MrBayesのinfileを順次作成して、順次MrBayesによる解析を実行する。各遺伝子ごとの解析結果から、BUCKyのinput file (*.in)をmbsumで作成する。

#!/usr/bin/sh
# multi_mb.sh

# 0. mrbayes mpiで動く設定にする // まだできていない <-- 不要
#mpirun -np 8 mb-mpi

# for $fname
# 1. .fasta(.aln) --> .nexに変換 & mrbayes block追加 // fasta_to_nex.pl使用
# 2. mrbayes実行 (mpi)
# 3. sedで$fnameの拡張子除く
# 4. mrbayesの結果ファイルを$fnameのdirectoryに移動
# done

mkdir mb_results		# 解析を行うdirectory

for fname in *.aln; do
# 1. sedで$fnameの拡張子除く
fname=`echo "$fname" | sed -E 's/\.aln//g'`		# 拡張子を除く

# 2. .fasta(.aln) --> .nexに変換 & mrbayes block追加 // fasta_to_nex.pl使用
./fasta_to_nex.pl $fname.aln > $fname.nex
mv $fname.nex mb_results

# 3. mrbayes実行 // MPIで実行すると早く動くが、たまにエラーが出て止まる。原因不明 ctrl-cでskipできる(手動だけど・・)
# set quitonerror=no;を.nex fileのmrbayes blockに入れることで上記のエラー停止を回避できる fasta_to_nex.plを修正 <-- 150325
mpirun -np 4 mb ./mb_results/$fname.nex 
#mb ./mb_results/$fname.nex		# mpiでなく実行

# 4. sedで$fnameのいろんな拡張子除く
fname=`echo "$fname" | sed -E 's/\.fasta.+//g'`

# 5. mbsum実行 (n=500 <-- burnin 50,000)
mbsum -n 501 -o ./mb_results/$fname.in ./mb_results/$fname.fasta.blastp.cds.masked.nex.run?.t

# 6. mrbayesの結果ファイルを$fnameのdirectoryに移動
mkdir ./mb_results/$fname
mv ./mb_results/$fname.* ./mb_results/$fname

done

# 7. bucky実行 // shell scriptとしては動かさない alpha valueなどいろいろ調整がある
# bucky ./mb_results/comp*/*.in