明日やろうは馬鹿野郎。
前回の続き。
では準備もできたということで、やっていきしょうshort tutrial
このチュートリアルは、ざっくりと決められたワークフローを組んでいくことでSnakemakeと仲良くなろうという趣旨のものです。
内容はバイオインフォマティクス分野からで、
bwa
によるアライメント
→ samtools sort
によるソート
→ samtools mpileup
による変異検出
→ 独自ツール(plot-quals.py
)による集計
となっています。
いいですね。ざっくり学べそうですね。
Prerequirement
まずは、データの準備を行います。
今回は公式のレポジトリから必要なデータをダウンロード、解凍します。
$ mkdir snakemake-demo $ cd snakemake-demo $ wget https://github.com/snakemake/snakemake-tutorial-data/archive/v5.4.5.tar.gz $ tar --wildcards -xf v5.4.5.tar.gz --strip 1 "*/data"
Step 1
次にSnakemake
というファイルを作成します。
これがSnakemakeのワークフローを記述するファイルとなります。
次からのステップではここにどんどん追記していきます。
$ touch Snakefile
このファイルを作成すると、dry-run が可能になります。
現時点での実行結果は、Nothing to be done.
です。
何も書いてないから当たり前ですね。
$ snakemake -n Building DAG of jobs... Nothing to be done.
Step 2
tree
コマンドでダウンロードしたdata
ディレクトリの中身を見てみましょう。
$ tree data/ (snakemake) data/ ├── genome.fa ├── genome.fa.amb ├── genome.fa.ann ├── genome.fa.bwt ├── genome.fa.fai ├── genome.fa.pac ├── genome.fa.sa └── samples ├── A.fastq ├── B.fastq └── C.fastq
今回は、data/samples
以下のfastqファイルをdata/genome.fa
にマップし、遺伝子変異を検出するパイプラインを作成します。
まずはbwa
というruleを作成します。
入力ファイルは
data/genome.fa
data/samples/A.fastq
出力ファイルは
mapped/A.bam
実行コマンドは
$ bwa mem {input} | samtools view -Sb - > {output}
です。
コマンドを指定しただけでは不十分で、実行環境にbwa
がインストールされ、PATHが通っている必要があります。
condaでインストール可能なソフトウェアであれば、conda: "envs/mapping.yaml"
みたいに依存パッケージをYAML形式で書いておくことで実行時に環境の作成、必要なパッケージのインストールを自動で行ってくれます。
envs/mapping.yaml
の中身はこんな感じ。
探しに行くcondaのチャンネルと必要なパッケージ名、バージョンを記述します。
channels: - bioconda - conda-forge dependencies: - bwa =0.7.17 - samtools =1.9
上記の内容をSnakefile
に記述するとこうなります。
rule bwa: input: "data/genome.fa", "data/samples/A.fastq" output: "mapped/A.bam" conda: "envs/mapping.yaml" shell: "bwa mem {input} | samtools view -Sb - > {output}"
dry-run(-n
オプション)を行うと、実行前にワークフローの概略を把握することができます。
$ snakemake --use-conda -n mapped/A.bam Building DAG of jobs... Conda environment envs/mapping.yaml will be created. Job counts: count jobs 1 bwa 1 [Thu Dec 5 09:29:51 2019] rule bwa: input: data/genome.fa, data/samples/A.fastq output: mapped/A.bam jobid: 0 Job counts: count jobs 1 bwa 1 This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.
問題なさそうなので実際に実行します。
-n
を外して実行すると、先ほどと違い実際に環境の構築、bwaの実行が行われ、結果が出力されます。
$ snakemake --use-conda -n mapped/A.bam Building DAG of jobs... Creating conda environment envs/mapping.yaml... Downloading and installing remote packages. Environment for envs/mapping.yaml created (location: .snakemake/conda/56da7cc7) Using shell: /bin/bash Provided cores: 1 (use --cores to define parallelism) Rules claiming more threads will be scaled down. Job counts: count jobs 1 bwa 1 [Thu Dec 5 09:35:34 2019] rule bwa: input: data/genome.fa, data/samples/A.fastq output: mapped/A.bam jobid: 0 Activating conda environment: /home/kimoton/snakemake-demo/.snakemake/conda/56da7cc7 [M::bwa_idx_load_from_disk] read 0 ALT contigs [M::process] read 25000 sequences (2525000 bp)... [M::mem_process_seqs] Processed 25000 reads in 1.594 CPU sec, 1.586 real sec [main] Version: 0.7.17-r1188 [main] CMD: bwa mem data/genome.fa data/samples/A.fastq [main] Real time: 1.906 sec; CPU: 1.734 sec [Thu Dec 5 09:35:40 2019] Finished job 0. 1 of 1 steps (100%) done Complete log: /home/kimoton/snakemake-demo/.snakemake/log/2019-12-05T093440.986213.snakemake.log
※ 以降煩わしいので標準出力は省略します。
これが完了すると、mapped/A.bam
が出力されているはずです。
Step 3
Step 2ではbwaの実行対象がA.bam
のみでした。
これをB.bam
, C.bam
についても実行できるようにするために、{sample}
変数を使って入出力ファイルを指定するように変更します。
rule bwa: input: "data/genome.fa", "data/samples/{sample}.fastq" output: "mapped/{sample}.bam" conda: "envs/mapping.yaml" shell: "bwa mem {input} | samtools view -Sb - > {output}"
Step 4
次のrule (sort
) を書いていきます。
入力ファイルは
mapped/{sample}.bam
出力ファイルは
mapped/{sample}.sorted.bam
実行コマンドは
$ samtools sort -o {output} {input}
実行に必要なsamtools
は、Step 2の時点でenvs/mapping.yaml
に記載していたため、特に別途作成する必要はありません。
ということでSnakefile
に下記を追記します。
rule sort: input: "mapped/{sample}.bam" output: "mapped/{sample}.sorted.bam" conda: "envs/mapping.yaml" shell: "samtools sort -o {output} {input}"
先ほど同様にまずはdry-runします。
snakemake --use-conda -n mapped/A.sorted.bam
問題なければ-n
を外して実際にrunします。
snakemake --use-conda mapped/A.sorted.bam
完了すると、mapped/A.sorted.bam
が出力されます。
Step 5
次に変異検出を行うruleを作成します。
このruleでは、サンプルを跨いだjoint callを行うため、AだけでなくB、Cも入力ファイルとします。
そのために、まずSnakefileの最上部で下記を定義します。
samples = ["A", "B", "C"]
expand
関数を用いることで、リスト形式を展開して引数に渡すことができます。
よって、入力ファイルは
- fa="data/genome.fa"
- bam=expand("mapped/{sample}.sorted.bam", sample=samples)
出力ファイルは
- calls/all.vcf
コマンドは
samtools mpileup -g -f {input.fa} {input.bam} | bcftools call -mv - > {output}
となります。
call
という名前でruleを追記します。
rule call: input: fa="data/genome.fa", bam=expand("mapped/{sample}.sorted.bam", sample=samples) output: "calls/all.vcf" conda: "envs/calling.yaml" shell: "samtools mpileup -g -f {input.fa} {input.bam} | " "bcftools call -mv - > {output}"
更に、このステップに必要なパッケージをenvs/calling.yaml
として作成します。
channels: - bioconda - conda-forge dependencies: - bcftools =1.9 - samtools =1.9
Step 6
パイプラインの最後に、変異検出結果の統計的なプロットを出力するようにします。
これには、scripts/plot-quals.py
という下記のようなスクリプトを使用します。
import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from pysam import VariantFile quals = [record.qual for record in VariantFile(snakemake.input[0])] plt.hist(quals) plt.savefig(snakemake.output[0])
このスクリプトを見るとわかるように、
スクリプト内ではsnakemake.input
という形でSnakefileで指定した入力ファイルをリスト形式で受け取ることができます。
このステップに必要なパッケージをenvs/stats.yaml
として作成します。
channels: - bioconda - conda-forge dependencies: - pysam =0.15 - matplotlib =3.1 - python =3.7
確認のため、この時点でもう一度実行してみましょう。
$ snakemake --use-conda plots/quals.svg
無事完了すると、plots/quals.svg
が出力されているはずです。
Step 7
ここまで、Snakemakeを実行する際には、引数にターゲットとなるファイルを指定していましたが、ターゲットのファイルを指定しない場合、SnakemakeはSnakefileの最初の行を実行しようとします。 この性質を使用して、defaultの処理を定義することができます。
最上部にall
という名前のruleを追記します。
rule all: input: "calls/all.vcf", "plots/quals.svg"
ここではoutput
およびcommands
は定義していません。
こうすることにより、
ファイルを指定しなかった場合のターゲットファイルをcalls/all.vcf
とplots/quals.svg
に指定することができます。
Step 8
最後のステップとして、さらに幾つかの工夫をワークフローに適用します。
自動レポーティング
Snakemakeは、HTML形式のレポートを、自動で生成することができます。
$ snakemake --report report.html
このレポートにはruntimeの統計情報であったり、私用しているデータのパスなんかが記載されます。
追加で、ワークフロー内で出力した画像を含めることもできます。
今回の場合、plots/quals.svg
をreport("plots/quals.svg", caption="report/calling.rst")
に書き換えることで、
quals.svg
をレポートに含めることができます。
caption
引数で指定しているreport/calling.rst
には、reStructuredText記法で出力ファイルについての説明を記載します(詳細:reStructuredText入門 — Sphinx 4.0.0+/2c7d64b94 ドキュメント)
レポート内で画像を表示するとこんな感じ
Thread数を増やす
前回の記事で述べたように、Snakemakeの特徴としてワークフロー全体のスケールアップが容易にできることが挙げられます。
スケールアップの際には、--cores
(--jobs
, -j
でも可)で使用するCPUコア数を指定してあげる必要があります。
指定したコア数を使用してSnakefile内のruleをガンガン並列に処理してくれるという寸法です。
デフォルトでは各ruleは単一のCPUで処理される(rule内の処理は並列化されない)のですが、rule内の処理を並列化させたいという場合があるかもしれません。
例えば今回の場合、リードのアライメントに使用しているbwa
はmultithreadでの実行に対応しています。
これを使わない手はありません。
-t
引数でThread数を指定することができるため、下記のように書き換えます。
threads: 8 shell: "bwa mem -t {threads} {input} | samtools view -Sb - > {output}"
このように、Thread数の指定をshell
ディレクティブ内でハードコーディングするのではなく、threadsディレクティブに分けることでrule内の処理を並列化する(コマンドラインで指定したcore数からそのthread数分割り当てる)ことが可能になります。
一時ファイルの使用
sort
ruleによってソートされたBAMファイルが出力されると、ソート前のBAMファイル(bwaファイルの出力ファイル)は不要になります。
このため、ルールbwa
の"mapped/{sample}.bam"
をtemp("mapped/{sample}.bam")
に置き換えて、出力を一時的なものとします。
結果
このチュートリアルで最終的に出来上がるSnakeFileは下記の通りになります。
defaultのall
が実行された際には、
bwa
→ sort
→ call
→ stats
と上から順に行われることになります。
これは別に上から順に実行しているわけではなく、「指定した出力ファイルにどの処理が必要か」というのが再帰的に辿られているだけです。 シンプルで分かりやすいですねー。
samples = ["A", "B"] rule all: input: "calls/all.vcf", "plots/quals.svg" rule bwa: input: "data/genome.fa", "data/samples/{sample}.fastq" output: temp("mapped/{sample}.bam") conda: "envs/mapping.yaml" threads: 8 shell: "bwa mem -t {threads} {input} | samtools view -Sb - > {output}" rule sort: input: "mapped/{sample}.bam" output: "mapped/{sample}.sorted.bam" conda: "envs/mapping.yaml" shell: "samtools sort -o {output} {input}" rule call: input: fa="data/genome.fa", bam=expand("mapped/{sample}.sorted.bam", sample=samples) output: "calls/all.vcf" conda: "envs/calling.yaml" shell: "samtools mpileup -g -f {input.fa} {input.bam} | " "bcftools call -mv - > {output}" rule stats: input: "calls/all.vcf" output: report("plots/quals.svg", caption="report/calling.rst") conda: "envs/stats.yaml" script: "scripts/plot-quals.py"
おわりに
再現性のあるスケーラブルな可搬性のある解析になっていたことが伝わったでしょうか。
個人的には
- ワークフロー、実行環境の定義ともにYAML形式で記述するため、可読性が高い
- 基本的に入出力ファイルを定義してつなぎ合わせるだけなので書き易い
- dry-run機能によりワークフロー概要の把握が可能
辺りが気に入りました。
一方で、レポートのフォーマットがどこまで弄れるのか、ログの出力は扱えるのかなんかが気になりますね。
dockerを使ったフローも見てみたい。
次回はもうちょい応用的な使用例を紹介したいと思います。