Snakemakeを試す - ② チュートリアルを実行

OSS bioinformatics Python

明日やろうは馬鹿野郎。

前回の続き。

では準備もできたということで、やっていきしょう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 ドキュメント）

レポート内で画像を表示するとこんな感じ f:id:kimoppy126:20191207164632p:plain

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数分割り当てる)ことが可能になります。

一時ファイルの使用

sortruleによってソートされた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を使ったフローも見てみたい。

次回はもうちょい応用的な使用例を紹介したいと思います。