はじめに

この１年ほどGoogle Cloud Platform (GCP) を使って来たのですが、GCPを使ったゲノム解析のパイプライン構築のベストプラクティスがずっとわからず、色々と調べても満足が行くものがなくて困っていました。ですが、数ヶ月前に「これだ！」と思うフレームワークを見つけて、早速使い始めました。以下のdsubというソフトウェアです。
https://github.com/googlegenomics/dsub

このdsubですが、Google Genomics Pipeline APIというもののラッパーのようです。（基本的には）入力データがGoogle Cloud Storage (GCS)にあるものとして、以下の手順を簡単に実行できるようにするソフトウェアです。

まず、Virtual Machine (VM) を立てる。
立てたVMに、GCSから入力データ（fastqファイルなど）、その他必要なデータ（リファレンスゲノムなど）を転送する。
転送が終わったら、処理（アラインメントなど）をDockerコンテナ上で実行する。
処理が終わったら、成果物（bamファイルなど）を、GCSに転送する。
転送が終わったら、VMを除去する。

事前にVMを立てたり、ソフトウェアのインストールなどをする必要などが全くなくて、タスクの実行に合わせて勝手にVMが立ち上がり、死んでくれるという動作がとても魅力的に感じました。ここでは、テストケースとして、STARを使ってシークエンスデータ（fastqファイル）をアラインメントして、samtoolsを使ってbamファイルを染色体の座標順に並べ替えるプログラムを作成する方法を紹介します。

dsubのタスクの規定

まず、dsubにおいては、以下の二つのファイルにより、「ジョブ」の定義がなされます。

実行環境のDockerイメージ
具体的に実行するシェルスクリプト（またはコマンド）

例えば、STARを実行する際には、例えば以下のようなDockerfileをbuildしてdocker-hubにpushします。

FROM ubuntu:16.04

MAINTAINER Yuichi Shiraishi <friend1ws@gmail.com>

RUN apt-get update && apt-get install -y \
    wget \
    bzip2 \
    make \
    gcc \
    zlib1g-dev

# STAR
RUN wget https://github.com/alexdobin/STAR/archive/2.5.3a.tar.gz && \
    tar xzvf 2.5.3a.tar.gz && \
    mv STAR-2.5.3a/bin/Linux_x86_64_static/STAR /usr/local/bin

# samtools
RUN wget https://github.com/samtools/samtools/releases/download/1.3.1/samtools-1.3.1.tar.bz2 && \
    tar jxvf samtools-1.3.1.tar.bz2 && \
    cd samtools-1.3.1/htslib-1.3.1 && \
    ./configure §§ make && make install && \
    cd ../ && ./configure --without-curses && make && make install

シェルスクリプトとしては、例えば以下のファイルを用意します。

star-alignment.sh

#!/bin/bash

set -o errexit
set -o nounset
set -o xtrace

OUTPUT_PREF=${OUTPUT_DIR}/${SAMPLE}

STAR_OPTION="--runThreadN 6 --outSAMstrandField intronMotif --outSAMunmapped Within --outSAMtype BAM Unsorted"

SAMTOOLS_SORT_OPTION="-@ 6 -m 3G"

/usr/local/bin/STAR --genomeDir ${REFERENCE} --readFilesIn ${INPUT1} ${INPUT2} --outFileNamePrefix ${OUTPUT_PREF}. ${STAR_OPTION}

/usr/local/bin/samtools sort -T ${OUTPUT_PREF}.Aligned.sortedByCoord.out ${SAMTOOLS_SORT_OPTION} ${OUTPUT_PREF}.Aligned.out.bam -O bam > ${OUTPUT_PREF}.Aligned.sortedByCoord.out.bam

/usr/local/bin/samtools index ${OUTPUT_PREF}.Aligned.sortedByCoord.out.bam

rm -rf ${OUTPUT_PREF}.Aligned.out.bam

これで、STAR + samtoolsのタスクの定義ができました。

dsubのタスクの実行

dsubにおけるジョブの実行方法ですが、

シェルスクリプトの引数のリスト（入力・出力データのGCSにおけるパス、その他オプションなど）
dsub実行コマンドファイル

の二つのファイルが必要となります。具体的には、例えば引数のリストは以下のようなファイルを用意します。

--env SAMPLE    --input INPUT1  --input INPUT2  --output-recursive OUTPUT_DIR   --input-recursive REFERENCE
Hela_wt_dox-    gs://genomon_rna_gce/sequence_data/small/Hela_wt_dox-/sequence1.txt     gs://genomon_rna_gce/sequence_data/small/Hela_wt_dox-/sequence2.txt     gs://friend1ws_genomon_rna_gce_test/test171125/star/Hela_wt_dox-       gs://genomon_rna_gce/db/GRCh37.STAR-2.5.2a
MCF-7   gs://genomon_rna_gce/sequence_data/small/MCF-7/sequence1.txt    gs://genomon_rna_gce/sequence_data/small/MCF-7/sequence2.txt    gs://friend1ws_genomon_rna_gce_test/test171125/star/MCF-7      gs://genomon_rna_gce/db/GRCh37.STAR-2.5.2a
DU145_150       gs://genomon_rna_gce/sequence_data/small/DU145_150/sequence1.txt        gs://genomon_rna_gce/sequence_data/small/DU145_150/sequence2.txt        gs://friend1ws_genomon_rna_gce_test/test171125/star/DU145_150  gs://genomon_rna_gce/db/GRCh37.STAR-2.5.2a
K562_150        gs://genomon_rna_gce/sequence_data/small/K562_150/sequence1.txt gs://genomon_rna_gce/sequence_data/small/K562_150/sequence2.txt gs://friend1ws_genomon_rna_gce_test/test171125/star/K562_150   gs://genomon_rna_gce/db/GRCh37.STAR-2.5.2a

–env, –input, –input-recursive, –output, –output-recursiveなどは「どのファイルがGoogle Cloud Storageからの入力、出力か？」を指定するタグを表しています（詳しくはdsubのドキュメントを参照）。そして、具体的なdsub の実行コマンドは以下のような形になります。また、 gs://genomon_rna_gce以下のデータは、公開設定にしているので（2017年11月26日現在）、こちらの入力ファイルはそのまま使うことができる（と思います。もし試した方がいらっしゃって、参照できなかったら、書き込みいただけると嬉しいです）。

dsub \
    --project [your_project] \
    --zones asia-northeast1-a \
    --disk-size 200 \
    --min-cores 6 \
    --min-ram 36 \
    --logging gs://friend1ws_genomon_rna_gce_test/test171125/logging \
    --image friend1ws/star-alignment \
    --tasks star-alignment-tasks.tsv \
    --script star-alignment-script.sh

これを実行することで、サンプルの個数だけのVMが立ち上がり、STARによるアラインメントが実行されます。

Google Genomics Pipeline APIのラッパーのdsubを使って、STARでRNA-seqのアラインメントを実行するのバッチジョブの書き方

はじめに

dsubのタスクの規定

dsubのタスクの実行