2021 5/5 コードの改行 (\) を除去
2021 5/6 説明を修正
Nextflowは2018年のアップデートでcondaに対応し（リンク）、nextflow側からcondaを呼び出してランできるようになりました。 そこで”Nextflowを使ってバイオインフォマティクスのツールを動かす”第2回の今回は、Nextflowにcondaの環境を作らせてランします。condaの環境作成までNextflowに任せることで、より再現性の高い計算パイプラインを構築することを目指します。
どのようにcondaの環境を指定するのかですが、Nextflow Documenttionでは、configファイルを書いてランする時にオプションで指定するか、”nextflow.config”というファイル名のYML形式のファイルに環境を書いて保存する方法などが説明されています。
詳細はこちら（優先順位あり）
https://www.nextflow.io/docs/latest/config.html
つまりnextflow.configをカレントに保存しておく事で、main.nfを叩いた際にnextflow.configが自動で認識されて設定が読み込まれるということです。
nextflow.configは例えば以下のようなファイルになります。
この設定ファイルにはparamsスコープとprofileスコープが記述してあって、前半のparamsスコープでは、 fastqのディレクトリ、出力ディレクトリ、freebayesのパラメータ（value）等を指定しています。上のようにparams { }で囲むか、params.minc =2などと書きます。これらのパラメータは、前日の記事ではmain.nfに記述していました。
paramsスコープに書かれたパラメータはラン時にオプションとして指定することができます（解説）。つまり、
と書かれているか、
と書かれていれば、デフォルトでは２（value）が$mincに読み込まれます。さらにコードを実行する時にオプションで”nextflow run main.nf --minc 5”と指定すれば、２が５に上書きされ、５（value）が読み込まれます。
後半のprofileスコープではcondaの設定ファイルのパス、dockerのイメージ名（後述）を指定しています。まずconda { }の部分ですが、condaの環境を書いたYMLファイルを読み込んでいます。
指定しているのは次のようなymlファイルです。これも指定されたパス（ここではカレントパス）に置いておく必要があります。
environment.yml
このYMLファイルは、
mamba env export > env.yaml
で作ることができます。condaの仮想環境を作って必要なツールを導入・テストし、問題なければ作って下さい。しかしこの記述だとosが変わる事で依存ライブラリの互換性が代わって環境を再現できない可能性が高くなります。
もっとシンプルに
と書けます。この書き方だと、依存関係に問題がない中で各ツールの最新バージョンが導入されます。
まとめると、condaプロファイルでランするには、メインのmain.nfの他に、設定ファイル：nextflow.config、環境ファイル：environment.ymlの２つのファイルすればいい事になります。
ではコードを書いてみます。
テストコード
今回は６組あるペアエンドfastqについて以下の計算を行う。
- process preprocessing：fastpを使ってペアエンドfastqのアダプタートリミングとクオリティフィルタリングを実行。
- process mapping：前処理したペアエンドfastqを受け取り、minimap2を使ってリファレンスゲノムにマッピング、elprepにパイプしてPCR dulicationをフィルタリングしてcoordinate_sorted.bamを出力。
- process variant_call：ソートされたbamを受け取り、freebayesを使ってvcfファイルを出力。
- process indexing：、vcfを受け取り、圧縮してindexを付ける。
ここでは１つずつバリアントを呼び出していますが、freebayesはジョイントコールに対応しているので、たくさんサンプルがある場合、まとめて呼び出して感度を上げることができます（複数サンプルを参照する事で頻度の少ないバリアントをノイズと区別してコールできる）。このコードはあくまで例と考えて下さい。
メインコード: variant_call_freebayes.nf
（実際のコードは一番下に貼ってあります）
nextflow.config
environment.yml
（実際のコードは一番下に貼ってあります）
環境作成
Nextflow自身にcondaの環境を作らせてランさせるため、利用するツールを予めインストールしておく必要は無くなります。ただしcondaだけは使えるようにしておきます。
conda -h
＜参考＞”-profile conda”だけ指定してランした場合、ymlファイルに書かれているツールだけが動作します。つまり、現在miniconda3/bin/やanaconda3/bin/にパスが通っているコマンドでも、-profile condaをつけてランするとパスは見えなくなります。ymlファイルに書かれいる環境がアクティブになるためです。再現性を確保するためには当然の仕様ですが、注意して下さい。
ディレクトリ構造
現在のパス（$baseDi）から見たディレクトリ構造が下の写真になります。入力データであるfastq.gzファイルはfastq ディレクトリに配置し、fasta形式のリファレンス配列はreferenceディレクトリに配置しています。
prefix はnextflow.configの
に記載されている通り、xxx_1.fastq.gzとxxx_2.fastq.gzになっている必要があります（１組以上あれば何組でも認識されます）。
ラン
準備ができたらnextflowのコードを実行します。-profileでcondaを指定します。
nextflow variant_call_freebayes.nf -profile conda
#freebayesの一部パラメータも指定できるように書いています。指定するなら
nextflow variant_call_freebayes.nf -profile conda --minc 4 --proidy 4
- -profile Choose a configuration profile
condaを指定した場合、まずworkにcondaの環境が作成されます（補足：あまりに依存関係が多いと、依存の解決が延々と続いて終わらないことがあります）
condaの環境作成が終わると計算が開始されます。
4つのプロセス全て正常終了。
executor > local (24)
[4d/45bf1a] process > preprocessing (SRR12340712) [100%] 6 of 6 ✔
[0e/6412dd] process > mapping (SRR12340712) [100%] 6 of 6 ✔
[ce/4b4d4b] process > variant_call (SRR12340712) [100%] 6 of 6 ✔
[02/e65cd3] process > indexing (SRR12340712) [100%] 6 of 6 ✔
Completed at: 04-5月-2021 19:51:43
Duration : 5m 26s
CPU hours : 0.7
Succeeded : 24
出力
流れを動画にしました。condaの環境作成に時間がかかっていますね。
次はdockerです。profileでdockerも指定しています。
このdockerプロファイルが指定されると、kazumax/variant_call:1.0というイメージに含まれているツールを使って計算が行われます。
イメージは予め作っておくか、公開してpullできるようにしておく必要があります。以前、ベースイメージ<kazumax/miniconda:1.0>を作って公開しました(解説)。これをベースに作ります。
#################################################
FROM kazumax/miniconda:1.0
LABEL description="variant-call-freebayes"
#################################################
RUN mamba install -c bioconda -qy samtools minimap2 fastp freebayes elprep2 && mamba clean -ay
#################################################
Dockerfileというファイル名で保存します (空のディレクトリを作ってその中で行う)。
ビルドしましょう。
docker build -t kazumax/variant_call:1.0 .
これで準備ができました。元のパスに戻ってランします。-profileでdockerを指定します。
nextflow variant_call_freebayes.nf -profile docker
- -profile Choose a configuration profile
これでdockerイメージを使ってランされます。指定したイメージ：<TAG>がローカルにも公共サイトにもない場合はエラーになります（サイズを減らす）。
2021 5/5 追記
依存する一部のツール（macosで動作するツール）はcondaを使い、残りのツールはdocker（pure linuxでないとビルドが難しいツール）という風に、２つ以上プロファイルをしてランすることも可能です。-profileでconda,dockerと指定します。
nextflow variant_call_freebayes.nf -profile conda,docker
これでcondaとdockerイメージを使ってランされます。
nextflow.configにはCPUコア数やメモリの使用量の上限も記述できます。例えばnextflow.configに以下を記載しておくと、１プロセスでのCPU使用数最大４つ、メモリは最大7GBになります。
process {
memory='4G'
cpus= 6
scratch = false
withLabel: big_mem_cpus {
cpus = 20
memory = '60G'
}
}
1プロセスでの制限になるので、例えば前処理プロセスが６プロセス動いている時は、36コア使う可能性があります。コードでは withLabel: big_mem_cpus {}が書いてあるので、process内で特別にlabel 'big_mem_cpus' が書いてある場合、そのプロセスのみ、CPUは最大２、メモリは最大60GB の指定になります。複数のプロセスがあって、計算パラメータを一括指定しつつ、ゲノムアセンブリのようなプロセスのみリソースを増やしたい時に便利です。ただし、process内にCPUsやmemoryが記載してある場合、そちらが優先されるようです。
感想
Nextflowは便利ですが、どんな時でもワークフローエンジンを使うのがベストとは限らない点に注意して下さい。計算環境や計算内容によって使い分けるのが良いと思います。
DolphinNextというwebインターフェイスも発表されています。時間があれば紹介します。
https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-020-6714-x
関連
1回目
dockerベースイメージ
簡単な軽量化
コメント
たくさんのファイルをランする先は、ストレージの空き容量に注意してください。わずかに空き容量の余裕があるように見えても、たくさんの中間ファイルが出来るため、思った以上にスペースを使います。
１、variant-call-freebayes.nf （メインコード）
#!/usr/bin/env nextflow
/*
* test code2
* Kazuma Uesaka
*/
Channel
.fromFilePairs( params.reads )
.ifEmpty { error "Cannot find any reads matching: ${params.reads}" }
.set { read_pairs_ch }
// message
println "reads: $params.reads"
println "genome: $params.genome"
// process
process preprocessing {
cpus 2
tag "$pair_id"
publishDir params.outdir1, mode: 'copy'
input:
tuple val(pair_id), path(reads) from read_pairs_ch
output:
set pair_id, "${pair_id}_QT_1.fq.gz" into read1_ch
set pair_id, "${pair_id}_QT_2.fq.gz" into read2_ch
set pair_id, "${pair_id}_fastp*"
"""
fastp -i ${pair_id}_1.fastq.gz -I ${pair_id}_2.fastq.gz -o ${pair_id}_QT_1.fq.gz -O ${pair_id}_QT_2.fq.gz -h ${pair_id}_fastp.html -j ${pair_id}_fastp.json -q 30 -n 5 -t 1 -T 1 -l 20 -w ${task.cpus}
"""
}
process mapping {
tag "$pair_id"
publishDir params.outdir2, mode: 'copy'
cpus 2
input:
path genome from params.genome
tuple val(pair_id), path(reads1) from read1_ch
tuple val(pair_id), path(reads2) from read2_ch
output:
set pair_id, "${pair_id}.bam" into bam_ch
"""
minimap2 -ax sr -t ${task.cpus} ${genome} $reads1 $reads2 | elprep filter /dev/stdin ${pair_id}.bam --mark-duplicates --remove-duplicates --filter-mapping-quality 0 --clean-sam --nr-of-threads ${task.cpus} --sorting-order coordinate --filter-unmapped-reads-strict
"""
}
process variant_call {
tag "$pair_id"
input:
path genome from params.genome
val minc from params.minc
val plo from params.ploidy
tuple val(pair_id), path(bam) from bam_ch
output:
set pair_id, "${pair_id}_freebayes.vcf" into vcf_ch
"""
samtools faidx ${genome}
freebayes -C $minc -u -p $plo -f ${genome} $bam > ${pair_id}_freebayes.vcf
"""
}
process indexing {
tag "$pair_id"
publishDir params.outdir3, mode: 'copy'
input:
tuple val(pair_id), path(vcf) from vcf_ch
output:
set pair_id, "${pair_id}_freebayes.vcf.gz*"
"""
bgzip $vcf
tabix -p vcf ${vcf}.gz
"""
}
２、nextflow.config （設定ファイル）
// parameters
params {
reads = "$baseDir/fastq/*_{1,2}.fastq.gz"
genome = "$baseDir/reference/assembly.fasta"
output = 'results'
outdir1 = "$output/trimmed_fastq/"
outdir2 = "$output/bam/"
outdir3 = "$output/variant-call/"
help = false
minc = 2
ploidy = 1
}
//profile
profiles {
conda {
process.conda = "$baseDir/environment.yml"
}
docker {
process.container = 'kazumax/variant_call:1.0'
docker.enabled = true
}
none {
// don't load any configration (use local $PATH)
}
}
３、environment.yml （condaの環境設定）
name: variant_call
channels:
- bioconda
- conda-forge
- defaults
dependencies:
- fastp
- samtools
- elprep
- minimap2
- freebayes
2021 6/4
test data
fastq/以下の2つのfastqをref/にあるゲノムにマッピングして、バリアントコールするパイプライン（real dataの一部領域のfastqを使用。例外処理もほとんどない簡単なコードです）
git clone git@github.com:kazumaxneo/Nextflow_test_code.git
cd Nextflow_test_code/Github_nextflow_variant_call/
#環境作成
mamba create -n test
conda activate test
mamba install -c bioconda fastp elprep samtools minimap2 freebayes -y
#run
nextflow variant_call.nf
結果はresults以下に保存される。
