VCFから遺伝研のスパコンでBEAST2（系統樹推定）を動かすまで。

遺伝研のスパコンでBEAST2を動かしたい！

かなり適当ですが悪戦苦闘したら一応動きました。
（7/28、VCF変換を間違えていたため訂正）

BEAST2とは？

なんかベイズとMCMCで祖先推定したり系統樹を書いてくれたりするすごいやつらしいです（まだ触ってから1週間しかたっていないので全然わかりません）。
http://www.beast2.org/
BEAGLEというライブラリとセットで動かします（VCFのimputationをしてくれる同名のソフトとは関係ない）。
このソフトの利点は、DIYABCのようなarlsumstat（Arlequin）を計算に用いるソフトにおける「大量のSNPsが書かれたarpファイルがあるときは動かない」という制限がないことです（これでかなりの時間を無駄にしました）。

手元に必要なもの

・Linux（またはWindowsとWSL、自分はMSLで動かしました。）
・Python3とgitコマンド
・VCFファイル（またはNEXUSファイル、こちらのほうが望ましい）
・表現型や年月のデータ（BEATUtiに登録するときに使います。）
・メモリつよつよのパソコン（自分は32GBです）
（ちなみにGPUを使うモードもありますがスパコンでは試していません。）

下準備（VCFファイルがある方）

まず、Python3でVCFをNEXUSにします。
いろいろ試したり、異常のあるVCFファイルと格闘していたりしてかなり手間取りましたがこれが最善だと思います。
ちなみにIndelsは無視されます。
追記：2021/7/20
このソフトだとVCFはbiallelic（変異が1塩基に対して1つだけ）にしないといけません。先にVCFtoolsで（なぜか染色体の途中で止まった）BCFtoolsで

bgzip YOUR_Multiallelic.vcf
tabix -f -p vcf YOUR_Multiallelic.vcf.gz
bcftools view -m2 -M2 YOUR_Multiallelic.vcf.gz -O v --threads 8 > ${YOURVCF}.vcf.vcf

訂正しました、前回の場合だとBEAST2が読み込んでくれません。（7/28）

YOURVCF=test #vcfの拡張子のない部分
git clone https://github.com/edgardomortiz/vcf2phylip
python ./vcf2phylip/vcf2phylip.py -i ${YOURVCF}.vcf -n -f
sed "s/*/N/g" ${YOURVCF}.min4.nexus > All_biallelic_PASS.nexus

BEAUtiでNEXUSファイルを専用のXMLに加工する

これは手元のパソコンでできます。自分はexeファイルで行いましたが、WSLのGUIが使える方はWSLでもできます。
手元のパソコンのメモリを最大限に生かすために、「./beast/bin/beauti」を以下のように書き換えます。nanoでもサクラエディタでもいいようです。
こちらからBEAUtiとBEASTがダウンロードできます。Windows版はこれで動きました。
http://www.beast2.org/

beatui

if [ -n "$BEAST_EXTRA_LIBS" ]; then 
  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx8g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
else    
  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx8g -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
fi
#これを
if [ -n "$BEAST_EXTRA_LIBS" ]; then 
  "$JAVA" -Dlauncher.wait.for.exit=true -Xmx230g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
else    
  "$JAVA" -Dlauncher.wait.for.exit=true -Xmx230g -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
fi
#このようにする

これで自分のメモリを最大限にしてBEAUtiの処理ができます。
Windowsの場合はbatファイルのJavaに「-Xmx30g」のオプションをつけてbatファイルから実行すると動くはずです。

./beast/bin/beauti #beautiの実行

BEAUtiのGUIの使い方は基本的には公式チュートリアルを見てください。
https://www.beast2.org/tutorials/

ここではXMLをYOUR.xmlにします。

スパコンへのBEASTのインストール

意外とここで苦労しました。現在のバージョンではGithubからそのまま持ってくると動かないみたいです。そのため、2.6.4バージョンをインストールして、アップデートする手法をとります。

wget https://github.com/CompEvol/beast2/releases/download/v2.6.4/BEAST.v2.6.4.Linux.tgz
tar -zxvf BEAST.v2.6.4.Linux.tgz

スパコンへのBEAGLEのインストール（高速化）

高速化はこれを参考にしました。
https://www.beast2.org/performance-suggestions/
https://github.com/beagle-dev/beagle-lib/wiki/LinuxInstallInstructions
高速化するにはBEAGLEというライブラリを使います。これがあるかどうかで速度がだいぶ変わります。

git clone --depth=1 https://github.com/beagle-dev/beagle-lib.git
cd beagle-lib
./autogen.sh
./configure --prefix=$HOME
make install

sed -i -e "s/TreeLikelihood/ThreadedTreeLikelihood/g" YOUR.xml #こうすると並列化できる？

BEASTのスクリプト改変

BEAUtiと同様にBEASTもjavaを使って動くため、これも同様に改変します。
変更方法はnanoでもいいですが、自分はWinSCPの機能を使いました。
250GBにしたらメモリを使いすぎたエラーが出たので230GBにしました。

beast.sh

if [ -n "$BEAST_EXTRA_LIBS" ]; then 
  "$JAVA" -Dlauncher.wait.for.exit=true -Xmx220g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
else    
  "$JAVA" -Dlauncher.wait.for.exit=true -Xmx220g -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
fi

実行ファイル

すべてが終わった後は実行するだけです。必要回数MCMCを回してから終わるかははっきり言ってわかりません……
参考までにshファイルを作りました。これをジョブとして投入すればいいようです。

beast.sh

#$ -S /bin/sh
#$ -cwd
#$ -V
#$ -j y
#$ -l s_vmem=256G,mem_req=256G
#$ -l d_rt=192:00:00
#$ -l s_rt=192:00:00
#上のオプションで8日間動かせます。
export MALLOC_ARENA_MAX=2 #JAVAを遺伝研に使うときのおまじない
export PATH=$HOME/lib:$PATH #不安なのでここにもパスを入れます。
export LD_LIBRARY_PATH=$HOME/lib:$LD_LIBRARY_PATH #BEASTがBEAGLEを見つけられるようにパスを入れます。
./beast/bin/packagemanager  -update 
./beast/bin/packagemanager   -add BDSKY #パッケージを使う時はこんな感じです。
./beast/bin/beast  -beagle_SSE  -overwrite YOUR.xml

これで以下のようなジョブの実行ファイルが生成されたら成功です。

実行ファイル

Installation aborted: BDSKY is already installed.

                        BEAST v2.6.5, 2002-2021
             Bayesian Evolutionary Analysis Sampling Trees
                       Designed and developed by
 Remco Bouckaert, Alexei J. Drummond, Andrew Rambaut & Marc A. Suchard

##略
         Sample      posterior     likelihood          prior
              0 -2.595493066E8 -2.595494123E8         1.4007 --

愚痴

バイオ（とくに集団情報を使うもの）のソフトって遺伝子の入力がVCFじゃなくて独自の規格持っているもの多すぎて難しいです……
あとAsperaが自分の環境ではうまくインストールできませんでした（たしかエラー19）。
使える方は方法を教えてくださるとありがたいです。