はじめに

※本記事はバイオインフォマティクス Advent Calendar 2019の20日目の記事です。

Snakemakeのインストールと概要
→ Snakemakeを試す - ① インストールと概要 - ばいばいバイオ

チュートリアルの実行
→ Snakemakeを試す - ② チュートリアルを実行 - ばいばいバイオ

Snakemakeでの並列化
→ Snakemakeを試す - ③ 並列化 - ばいばいバイオ

①～③と、Advent Calendar に備え着々と準備を続けてきたのですが、
当日に全然時間がとれず今回は少し薄い内容となっております、どうかお許しを。

Snakemake vs Nextflow

バイオインフォ界隈のワークフローツールはSnakemakeだけでなく、CWL、Galaxy、Nextflowなんかが含まれてきます。多種多様であれがいい、これがいい、そもそもワークフローはいらないだの色々と騒がれています。

少し古い引用ですが、NCBIのハッカソンでこんな議論が繰り広げられていたらしいです。

CWLは労働集約的すぎるため、ほとんどすべてのメンバーに却下された。 CWLの経験を持つ少数の人々は、扱いの困難さ、学習コストの高さを問題点として挙げていた。
SnakemakeはNextflowよりも柔軟性が低いという理由で却下された。多くのユーザーはPythonを使用していると考えていたが、そうではない人もいた。
Nextflowは、入力と出力を管理するだけの定義言語があらゆる言語のAPIで簡単にラップできるようになっていたため、採択された。
Jupyterは、まだすべての言語で扱えるものではない。将来的には魅力的な選択肢となるかもしれない。

当時はなんとも言えない理由でNextflowが採択されたようです。

参考：NCBI Hackathons discussions on Bioinformatics workflow engines

Snakemakeのいいところ

自分はCWL、Galaxyを少しかじった程度しかワークフローを触ったことはなく、pythonでごりごり書くことが多かったのですが、巷ではやり始めたワークフローツールを使えばこれもっと楽にかけるんじゃね。。と思い、最近ではSnakemakeを信仰しています。

個人的に思うSnakemakeを使うべき人は下記の通りです。

Pythonが好き、Pythonでワークフローを記述したい
小難しい文法を勉強したくない（学習コストをかけたくない）
プログラミング初学者が多いチームで使いたい
数TBレベルの大規模処理は扱わない

ほかのワークフローツールと比べて、いろんな意味で軽量なツールかと思います。

GATK best practiceをSnakemakeで

dna-seq-gatk-variant-calling という便利なモノを見つけてしまいました。

これは、Snakemakeの作者であるJohannes Köster氏が作成したGATK best practiceのワークフローです。
彼はBiocondaプロジェクトのfounderや、Rust Bioの著者でもあるみたいです。

ワークフロー概略

Snakemakeのレポーティング機能は、下記のようなワークフロー図を出力してくれます。
少し画質が荒くて恐縮ですが、この図を張り付けておきます。

見てわかる通り、このワークフローはGATK best practiceに忠実なフローとなっています。
ワークフローの具体的な内容は本家を参考にしてください。

入力ファイル

units.tsvとして下記のように定義しています。
PLATFORMとか指定できるようになっていますが、レポーティング用途として使われているのみです。

sample	unit	platform	fq1	fq2
A	1	ILLUMINA	data/reads/a.chr21.1.fq	data/reads/a.chr21.2.fq
B	1	ILLUMINA	data/reads/b.chr21.1.fq	data/reads/b.chr21.2.fq
B	2	ILLUMINA	data/reads/b.chr21.1.fq

https://github.com/snakemake-workflows/dna-seq-gatk-variant-calling/blob/master/units.tsv

設定ファイル

config.yamlにて解析のオプションを設定します。
base recalibration を行うか否か、snv, indelごとのフィルタリング条件、trimmomaticやpicard, GATKのオプションが直感的に指定できるようになっています。

https://github.com/snakemake-workflows/dna-seq-gatk-variant-calling/blob/master/config.yaml

準備

このワークフローではBWAで使用するindex及びGATKが使用する.dictファイルが作成されていることが前提となっています。

それぞれリファレンスゲノムを指定して下記の通り実行することで作成できます。

bwa indexの実行

$ bwa index data/ref/genome.chr21.fa                                                                                         (snakemake)

picard CreateSequenceDictionaryの実行

$ picard CreateSequenceDictionary R=data/ref/genome.chr21.fa  O=data/ref/genome.chr21.dict

実行

Snakemakeには便利なdry-run機能があります。
重い処理のワークフローを実行する前にはまずdry-runして概要が正しいか確認するようにしましょう。

$ snakemake --use-conda -n

実際にワークフローを実行する際は、--coresを指定して並列に処理を走らせることができます。

$ snakemake --use-conda --cores 12

中で使われている便利な機能たち

ログ出力

ログの出力はlogディレクティブに出力ファイル名を指定するだけ。
標準出力が吐き出されます。

明示的なパラメータを含む場合はshellディレクティブ内で指定してあげます。

rule abc:
    input: "input.txt"
    output: "output.txt"
    log: "logs/abc.log"
    shell: "somecommand --log {log} {input} {output}"

とっても簡単です。

参考：log-files

モジュール化

Snakemakeでは、includeディレクティブを使って処理のモジュール化を行うことができます。
ワークフローが大規模になってきた際には積極的に使っていきましょう。

使い方はいたって簡単。

下記のように含めたいSnakefileのパスをincludeディレクティブに指定するだけ。
モジュール内では通常のSnakefileと同じシンタックスが使用できます。

include: "path/to/other.snakefile"

参考：moduralization

wrapper

wrapperディレクティブを使うと、一般的な解析ツールであれば自分でコーディングせずにSnakemakeのワークフローの一部として使用できます。
例えばbwaだとこんな感じ。

rule bwa_mem:
  input:
      ref="data/genome.fa",
      sample=lambda wildcards: config["samples"][wildcards.sample]
  output:
      temp("mapped_reads/{sample}.bam")
  log:
      "logs/bwa_mem/{sample}.log"
  params:
      "-R '@RG\tID:{sample}\tSM:{sample}'"
  threads: 8
  wrapper:
      "0.15.3/bio/bwa/mem"

このディレクティブが指定されると、呼び出し時に要求されたバージョンのラッパーを自動的にダウンロードしてくれます。
--use-condaと組み合わせて使うと、wrapper内で必要なツールの仮想環境へのインストールも自動的に行ってくれます。
割と頻繁に更新されてるし簡単にオプションのカスタムできるし、とっても便利。

まだ数十ツールしかラインナップがないので優秀な皆さん、committerチャンスです。

参考： tool-wrappers