はじめに
この１年ほど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
シェルスクリプトとしては、例えば以下のファイルを用意します。
#!/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によるアラインメントが実行されます。