LoginSignup
0
0
@mhiromi

遺伝研スパコンでANIを使う

Last updated at Posted at 2023-11-21

Average Nucleotide Identity (ANI)をスパコンのクラスタ環境下で使う方法を記す。ANIの実装はいくつかがあるが、オリジナルのソースコードは古すぎて使えず、condaのpyaniが最終的に使えた。https://github.com/widdowquinn/pyani

しかしpyaniも罠だらけなので注意。pyaniは内部でマルチプロセシングを利用しており、デフォルトだとavailableなCPUとメモリを全部使って計算してしまう系の迷惑ツールである。従って、スパコンで利用するには、使用するリソースの指定がとても重要である。このページの内容は試行錯誤の結果なので、もっとスマートなやり方があるかもしれないことに留意。

miniconda3のインストール

スパコンの説明を参照。https://sc.ddbj.nig.ac.jp/software/python/

pyaniのドキュメントの罠

使用上の注意

  • GitHub上には開発版の0.3と安定版の0.2がある。どうやらドキュメントが、0.3と書いてあるけれども0.2の内容だったり、0.2のドキュメントのリンクがどこだかよく分からなかったりする。
  • 0.2でも実装されているオプションの説明が、0.2のドキュメントにはなく0.3にちょっと書いてあることもあるので、google検索などで辿り着いた場合は注意が必要。ヘルプを呼び出した方が分かりやすかったりする。
  • そもそもドキュメントやらに書いていることが間違っていて、エラーの内容から正解（？）に辿り着けることも...。

0.2のREADMEのリンクはこちら https://github.com/widdowquinn/pyani/blob/master/README_v_0_2_x.md

pyaniのインストール

自分のcondaはpython 3.9だったので、それに合わせてインストール

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
conda create --name pyani_env python=3.9.1 -y

を実行。
この後、

conda activate pyani_env

でpyaniが使えるようになる。

重要！リソース使用量の見積り

pyaniの起動時とaverage_nucleotide_identity.pyのオプションで制御する。

(1) condaの起動に関しては、インタラクティブノードで使用する場合は、

export OMP_NUM_THREADS=1

を.bashrcに書いておく。これはUGEのジョブスクリプトに書いてもよい（ジョブスクリプトの例は後述）。

(2) 加えて、average_nucleotide_identity.pyでは後述のように、--workersオプションを指定する（UGEを使用する場合には、-pe def_slotとセットで）。

そうしないと、インタラクティブノードのログイン時やジョブスクリプトで宣言したリソースよりも多い量をあるだけ使って実行してしまって、ノードに負荷をかけてしまい怒られが発生する。例えば、UGEで-l s_vmem=32G,mem_req=32Gでジョブを投げたときに、実際にはそのノードにあった844GBのメモリを全部使ってジョブが落とされていた。怖すぎる。

ちなみに、qreportによるとaverage_nucleotide_identity.pyの実行には、メモリは最大でも4.3GB程度の使用であったので、安全のため多めに見積もったとしても、-l s_vmem=8G,mem_req=8Gでよさそう。

ANIのオプション

手抜きドキュメントよりもaverage_nucleotide_identity.py -hのヘルプ画面のほうがマシなときもある。しかしこのヘルプを出すのにもちょっと反応が重くて、せっかく出てきたヘルプの記述だけではよく分からないこともある...。

(pyani_env) $ average_nucleotide_identity.py -h
usage: average_nucleotide_identity.py [-h] [--version] -o OUTDIRNAME -i INDIRNAME [-v] [-f] [-s FRAGSIZE]
                   [-l LOGFILE] [--skip_nucmer] [--skip_blastn] [--noclobber]
                   [--nocompress] [-g] [--gformat GFORMAT] [--gmethod {mpl,seaborn}]
                   [--labels LABELS] [--classes CLASSES]
                   [-m {ANIm,ANIb,ANIblastall,TETRA}]
                   [--scheduler {multiprocessing,SGE}] [--workers WORKERS]
                   [--SGEgroupsize SGEGROUPSIZE] [--SGEargs SGEARGS] [--maxmatch]
                   [--nucmer_exe NUCMER_EXE] [--filter_exe FILTER_EXE]
                   [--blastn_exe BLASTN_EXE] [--makeblastdb_exe MAKEBLASTDB_EXE]
                   [--blastall_exe BLASTALL_EXE] [--formatdb_exe FORMATDB_EXE]
                   [--write_excel] [--rerender] [--subsample SUBSAMPLE] [--seed SEED]
                   [--jobprefix JOBPREFIX]
optional arguments:
 -h, --help      show this help message and exit
 --version       show program’s version number and exit
 -o OUTDIRNAME, --outdir OUTDIRNAME
            Output directory (required)
 -i INDIRNAME, --indir INDIRNAME
            Input directory name (required)
 -v, --verbose     Give verbose output
 -f, --force      Force file overwriting
 -s FRAGSIZE, --fragsize FRAGSIZE
            Sequence fragment size for ANIb (default 1020)
 -l LOGFILE, --logfile LOGFILE
            Logfile location
 --skip_nucmer     Skip NUCmer runs, for testing (e.g. if output already present)
 --skip_blastn     Skip BLASTN runs, for testing (e.g. if output already present)
 --noclobber      Don’t nuke existing files
 --nocompress     Don’t compress/delete the comparison output
 -g, --graphics    Generate heatmap of ANI
 --gformat GFORMAT   Graphics output format(s) [pdf|png|jpg|svg] (default pdf,png,eps meaning three
            file formats)
 --gmethod {mpl,seaborn}
            Graphics output method (default mpl)
 --labels LABELS    Path to file containing sequence labels
 --classes CLASSES   Path to file containing sequence classes
 -m {ANIm,ANIb,ANIblastall,TETRA}, --method {ANIm,ANIb,ANIblastall,TETRA}
            ANI method (default ANIm)
 --scheduler {multiprocessing,SGE}
            Job scheduler (default multiprocessing, i.e. locally)
 --workers WORKERS   Number of worker processes for multiprocessing (default zero, meaning use all
            available cores)
 --SGEgroupsize SGEGROUPSIZE
            Number of jobs to place in an SGE array group (default 10000)
 --SGEargs SGEARGS   Additional arguments for qsub
 --maxmatch      Override MUMmer to allow all NUCmer matches
 --nucmer_exe NUCMER_EXE
            Path to NUCmer executable
 --filter_exe FILTER_EXE
            Path to delta-filter executable
 --blastn_exe BLASTN_EXE
            Path to BLASTN+ executable
 --makeblastdb_exe MAKEBLASTDB_EXE
            Path to BLAST+ makeblastdb executable
 --blastall_exe BLASTALL_EXE
            Path to BLASTALL executable
 --formatdb_exe FORMATDB_EXE
            Path to BLAST formatdb executable
 --write_excel     Write Excel format output tables
 --rerender      Rerender graphics output without recalculation
 --subsample SUBSAMPLE
            Subsample a percentage [0-1] or specific number (1-n) of input sequences
 --seed SEED      Set random seed for reproducible subsampling.
 --jobprefix JOBPREFIX
            Prefix for SGE jobs (default ANI).

UGEでの実行例

最終的に使用しなかったオプションもあるので、うまくいった例を記す。

ジョブスクリプトを投げるノードで、conda activateをしてからジョブを投げる。

conda activate pyani_env
qsub ani.sh

####  ani.shの例 ####
#!/bin/sh
#$ -S /bin/sh
#$ -cwd

## メモリ8GBあれば起動に足りる
#$ -l s_vmem=8G 
#$ -l mem_req=8G

## 合計メモリ量次第ではmedium (今回は64GB)
#$ -l medium

## ANIの--workersのオプションと数を同じにすると良さそう。
#$ -pe def_slot 8

## 現在の自分の環境(conda)をUGEへ引き継ぐ
#$ -V     

## 利用するスレッド数を1にしておく
export OMP_NUM_THREADS=1

## システム側/opt/のcondaのパスが優先されることがあったため、絶対パスで指定
/home/hoge/miniconda3/bin/conda activate /home/hoge/miniconda3/envs/pyani_env

## condaのエラーが出るときは、infoを出せるか確認してみる
## /home/hoge/miniconda3/bin/conda info

## ANIをBlastで計算しPDFで出力 
average_nucleotide_identity.py -v -f -i . -m ANIb -o ANIb_output -g --gformat pdf　--classes Target_ani_data/Target_genomes.txt --labels Target_ani_data/Target_genomes.txt　 --workers 8

## ANIをMUMerで計算しPDFで出力
average_nucleotide_identity.py -v -f -i . -m ANIm -o ANIm_output -g --gformat pdf　--classes Target_ani_data/Target_genomes.txt --labels Target_ani_data/Target_genomes.txt --workers 8

average_nucleotide_identity.pyのオプション

--workers: 重要なオプション。UGEの-pe def_slotと同じ数を指定しておくと良さそう。

他のオプションでよく分からなかったところ。
--classes と --labels: デフォルトでは、fastaのファイル名(NC_002696.faのNC_002696の部分）が配列名として図に出力される。fastaファイル名そのままの出力で構わない時には、このオプションは指定不要。
もしもラベルを自分で指定したい場合は、1列目がファイル名, 2列目が表示させたい文字列のタブ区切りのファイルを作成して、classes, labelsでこの作成したファイルを指定する。例えばこんな感じ。

NC_002696	C.crescentus
NC_014100	C.segnis
NC_011916	C.crescentus
NC_010338	Caulobacter sp.

ちなみに、githubにはclassesやlabelsのサンプルファイル(labels.tab)が置いてあるが、こちらは3列で（1列目がhashになっている）、エラーメッセージを吐いて終了する。ソースコード上は、行毎にtabでsplitして、1列目は無視して2列目と3列目をdictに突っ込む仕様らしいが、なぜか3列の時に失敗して2列の時に使用できるみたい。なんで？

それ以外のオプションは公式どおり。

  • "-i ." カレントディレクトリの全入力ファイルを利用して計算する。
  • "-o" 出力ディレクトリ名の指定
  • "-f" 出力の上書き許可
  • "-m" 計算に使うプログラムの指定。BLASTであればANIb, MUMerであればANIm
  • "-gformat pdf" 結果はPDFで出力。他にPNGなど指定可

ANIを計算済みで、plotだけし直したいときは--rerenderオプションを使う。

$ average_nucleotide_identity.py -v -f -i . -m ANIm -o ANIm_output -g --gformat pdf --rerender

その他検討したオプション

--scheduler SGE, --SGEGroupsize
最終的にSGE/UGE系のオプションは使わなかった。子プロセスをSGEのアレイジョブに投げるっぽいんだけど、ここまでで十分カオスだったので、力尽きた...

--SGEargs
説明がないのでとりあえず、

--SGEargs "-l s_vmem=8G,mem_req=8G -pe def_slot 4"

みたいに複数のオプションをダブルクォートして指定したら、

average_nucleotide_identity.py: error: argument --SGEargs: expected one argument

というエラーが出たので使うのをやめた。 ソースを読めば分かるかもしれない。

NCBI taxonomy IDを指定すると自動で配列を取ってくるツールもあるが、バッチサイズなどを変えてもDLが途中で止まってしまうので、諦めた。

$ genbank_get_genomes_by_taxon.py -o target -v -t taxonid -l target_dl.log --email xxx@xxx.ac.jp -f --retries 400 --batchsize 100000

他にもオプションがあるが使いこなせてないので、分かる人いたら教えてください。有名ツールなのに使いこなすの難しい（涙）。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0