はじめに
分量が多くなってきたので分割しました。Singularityについては遺伝研スパコンで Singularity とか CWL を使う上での tipsをご参照ください。ここではCWLについて書きます。
CWL (Common Workflow Language) は名前の通り解析や定型作業のワークフローの手順を記述したもので、仮想コンテナを利用することで実行環境によらず再現性のある結果が得られることが大きなメリットだと考えています。
CWL についてより詳しく知りたい場合にはCWL Start Guide JP を、CWL を書いてみたいという人には@tm_tn氏の雑に始める CWL!をおすすめします。特に、
- 自分が作成したワークフローを他の人にもつかってもらいたい
- クラウドを含めた様々な環境でワークフローを実行したい
といった目的には合致すると思います。また、しばらく時間がたってから解析をやり直そうとしたらいつの間にか実行環境が変わっていてツールが動かない、といった悩みからも解放されるはずです。
ここでは、遺伝研スパコンを使って CWL を実行するための環境構築および実際にワークフローを実行する手順を紹介します。
実行環境の構築
遺伝研スパコンで CWL で書かれたワークフローを実行するには下記が必要です。
-
Python
次に挙げる cwltool をインストール・実行するために必要です。公式サイトによると Python 2.7 もサポートしているようですが、Python 3.5 以降をおすすめします。遺伝研スパコンには/usr/local/pkg/python
以下に Python 2.7 系と 3 系の両方がインストールされていますがこれを用いるよりも、pyenv や anyenv を使ってホームディレクトリに好きなバージョンを指定してインストールしてしまうのが良いと思います。Anaconda or Miniconda でも可です。これらのいずれも使用したことがないのであれば、Miniconda を導入しておくと後の手順が簡単です。
以下の手順ではpip
またはconda
コマンドが使用できることを前提とします。 -
cwltool
CWL の実行エンジンです。実行エンジンは他にもいくつかありますが標準実装 (reference implementation) とされています。インストールは、
pip install cwltool
で可能です。Anaconda/Miniconda を使用している場合には、
conda install cwltool
でもインストールできます。(2019/12/18追記 anaconda/miniconda のバージョンによっては conda を使用した場合 conflict が発生することがあるので、pip を使うことを推奨)
-
Singularity
CWL は通常 Docker コンテナと組み合わせて動作しますが、遺伝研スパコンでは Docker が使えないためかわりに Singularity を用います。スパコンにはインストール済みですので特に準備は必要ありません。CWL 実行時には自動で必要な Docker イメージを取得して Singularity のコンテナイメージへの変換が行われます。 -
Node.js
JavaScriptの実行エンジンです。一部の CWL ファイルは部分的に JavaScript を使っているため必要になります。インストールには nodeenv を用いるか、Anaconda/Miniconda を使用している場合にはconda
コマンドを使用するのが簡単です。conda
コマンドを使用する場合のインストール方法は以下の通りです。なお、遺伝研スパコン以外の環境で Docker を使って CWL を実行する場合には、実行時に Node.js の Docker イメージを自動的に取得して用いるためインストールは不要です。
conda install nodejs
実行してみる
例として CWL で書かれたイネの変異解析のためのワークフローを実行してみます。これは、Rice Annotation Project Database (RAP-DB) で公開されている標準解析パイプラインを CWL 化したものです。イネ用に作られたものですが、最後のバリアントのアノテーションを行うステップの手前までは、二倍体であれば他の生物種にも適用することができます。
ワークフローファイルの取得
遺伝研スパコンにログインして、qlogin
コマンドでインタラクティブノードに移動します。Singularity のコンテナイメージを生成する際に比較的大量のメモリを必要としますので、次のように要求メモリを指定してログインします(あるいはコマンド実行時にqsub
コマンドを使ってメモリ指定してバッチジョブとして投入します)。
qlogin -l mem_req=16G,s_vmem=16G
ログインできたら適当なディレクトリに移動して、ワークフローのファイルを取得します。
$ git clone https://github.com/nigyta/rice_reseq.git
$ cd rice_reseq/
$ ls
README.md docs scripts test tools workflows
tools
ディレクトリにはワークフロー中の各実行コマンドの動作を定義したファイルが含まれています。workflows
ディレクトリにはそれらを組み合わせてワークフローの一連の処理を定義したファイルが含まれます。ワークフローを実行する際には基本的にはworkflows
ディレクトリ内のファイルを用いることになります。公開されている CWL ワークフローの多くも同じようなディレクトリ構造になっています。
$ ls workflows/
bam2vcf.cwl prepare_reference.cwl rapdb-pipeline_wo_tasuke.cwl snpeff_all.cwl
fastq2bam.cwl rapdb-pipeline.cwl read_preprocessing.cwl test_samtools.cwl
実行手順
以下では参照ゲノム配列の indexing を行うprepare_reference.cwl
と、参照配列へのリード配列のアラインメントを行うfastq2bam.cwl
を実行する例を示します。
参照ゲノム配列はRAPDBからダウンロードします。また、リード配列はtest
ディレクトリに含まれている動作確認用のファイルを用います。
-
参照ゲノム配列の取得
同梱のスクリプトファイルを使って RAP-DB から必要な参照データを取得し
reference
というディレクトリに保存しておきます。$ ./scripts/download_reference_files.sh $ ls -lh reference/ total 1138824 -rw-r--r-- 1 tanizawa staff 365M 3 11 11:31 IRGSP-1.0_genome_M_C_unanchored.fa -rw-r--r-- 1 tanizawa staff 142M 3 11 11:31 RAP-DB_MSU_concatenated_for_snpEff.gtf -rw-r--r-- 1 tanizawa staff 46M 3 11 11:31 RAP-DB_MSU_concatenated_protein.fa
使用するのは一番上の参照ゲノム配列の FASTA ファイルだけです。残りの2つは snpEff を使用して変異のアノテーションを行うときに使用するためのものなので今回は必要ありません。
-
動作確認用データについて
$ cd test/ $ ls rapdb-pipeline.test.job.yaml read1.fq.gz read2.fq.gz
今回は
read1.fq.gz
とread2.fq.gz
を使って動作確認を行います。rapdb-pipeline.test.job.yaml
は次項以下で説明する`ジョブファイル`と呼ばれるもので、今回は使用しません。 -
参照ゲノム配列の indexing
prepare_reference.cwl
を用いて参照ゲノム配列の indexing を行います。後ほど、このファイルの中を覗いてどのような処理が行われているかみてみましょう。まずは、-h
をつけて cwltool で実行しヘルプを確認します。$ cwltool ../workflows/prepare_reference.cwl -h usage: ../workflows/prepare_reference.cwl [-h] --reference REFERENCE [job_order] positional arguments: job_order Job input json file optional arguments: -h, --help show this help message and exit --reference REFERENCE FASTA file for reference genome
job_order
という名称で示されるオプションは、'ジョブファイル' と呼ばれるワークフローへの入力ファイルやパラメータを記述したファイルを指定するためのものです(ジョブの実行順序?と勘違いしてしまったが、おそらくワークフローへの入力を指示するためのファイルという意味だと思われる)。ジョブファイルを用いて実行する方法は次項(手順4)で紹介するので、ここではコマンドラインから引数を与える方法で実行します。必須の引数は参照ゲノムの FASTA ファイルを指定する--reference
だけです。cwltool --singularity --outdir index ../workflows/prepare_reference.cwl --reference ../reference/IRGSP-1.0_genome_M_C_unanchored.fa
**遺伝研スパコンで CWL を実行する上で一番のポイントとなるのが
--singularity
オプションを指定して Docker のかわりに Singularity を使うようにすることです。**この一点を除けば他の環境で CWL を実行するのと違いはありません。
なお、--outdir
オプションを指定して結果ファイルの出力先としてindex
という名称を指定しました。何も指定しない場合にはカレントディレクトリに結果ファイルが生成されます。出力先のディレクトリは自動で生成されます。終了すると
index
ディレクトリに元のFASTAファイルのコピーとともに各種 index ファイルが出力されています。また、カレントディレクトリに Singularity のコンテナイメージファイル (*.img
) が生成されています。先に述べたように Singularity のイメージを生成する際に大量のメモリが必要となりますので、"Out of memory" で失敗してしまった場合には、qlogin 時の要求メモリを増やして実行してください。同じディレクトリから実行する場合にはイメージファイルは再利用されるため実行時間が短縮され、メモリも少量で済みます (qlogin 時の-l mem_req=16G,s_vmem=16G
は不要になります)。 -
リード配列のアラインメント
続いて
fastq2bam.cwl
を使ってリード配列のアラインメントを行います。ここでは前項とは異なりジョブファイルを利用して実行する方法を紹介します。まずは cwltool の--make-template
オプションを使いジョブファイルの雛形を生成します。$ cwltool --make-template ../workflows/fastq2bam.cwl > fastq2bam.job.yaml
生成されたジョブファイル
fastq2bam.job.yaml
には必要な入力パラメータを記載する項目があります。これを編集し、下記のようにしました。threads: 4 # default value of type "int". reference: # type "File" class: File path: index/IRGSP-1.0_genome_M_C_unanchored.fa outprefix: TEST # type "string" fastq2: # type "File" class: File path: read2.fq.gz fastq1: # type "File" class: File path: read1.fq.gz
おもな入力項目は参照配列の FASTA ファイルとアラインメントするリード配列のファイルのパスです。入力ファイルのパスは、絶対パスあるいはジョブファイルのある場所から見た相対パスで記載します。
outprefix
はサンプル名を識別するための文字列、threads
に並列実行時のスレッド数です。スレッド数に4を指定しているため、インタラクティブノードでそのまま実行すると他のユーザーの迷惑になります。バッチジョブのスクリプトファイルを作って
qsub
コマンドでジョブを投入しましょう。
バッチファイル (cwl_test.sh
の中身は以下の通りです)#$ -S /bin/bash #$ -pe def_slot 4 #$ -cwd #$ -l mem_req=8G,s_vmem=8G cwltool --singularity --outdir result ../workflows/fastq2bam.cwl fastq2bam.job.yaml
バッチファイルの準備ができたらジョブ投入 (
qsub cwl_test.sh
)しましょう。ジョブが完了すると
result
ディレクトリに Picard で重複除去処理 (Markduplicates) を行った後の BAM ファイルとそれに関連するファイルが結果ファイルとして出力されます。以上で実行方法の紹介を終わります。興味があれば README を参考に
rapdb-pipeline.cwl
を使ってアラインメントと変異の同定・アノテーションを含めた一連の工程をコマンド一発で実行してみてください。
ワークフローの中身を見てみよう
自分で CWL を書かず他の人が書いた CWL を利用するだけであっても、どのような処理が行われているか大雑把に理解しておくことは大切です。workflows
ディレクトリ内のファイルを見るとワークフローの中でどのような処理が行われているかを確認することができます。以下は、prepare_reference.cwl の中身です。
#!/usr/bin/env cwl-runner
cwlVersion: v1.0
class: Workflow
requirements:
StepInputExpressionRequirement: {}
InlineJavascriptRequirement: {}
inputs:
reference:
type: File
doc: FASTA file for reference genome
steps:
bwa_index:
run: ../tools/bwa-index.cwl
in:
reference: reference
out: [amb, ann, bwt, pac, sa, log]
create_dict:
run: ../tools/picard-CreateSequenceDictionary.cwl
in:
reference: reference
out: [dict]
create_faidx:
run: ../tools/samtools-faidx.cwl
in:
fasta: reference
out:
[fai]
bind_fasta_index:
run: ../tools/util-bindFastaIndex.cwl
in:
fasta: reference
fai: create_faidx/fai
amb: bwa_index/amb
ann: bwa_index/ann
bwt: bwa_index/bwt
pac: bwa_index/pac
sa: bwa_index/sa
dict: create_dict/dict
out:
[fasta_with_index]
outputs:
fasta_with_index:
type: File
outputSource: bind_fasta_index/fasta_with_index
重要なのはinputs
、steps
、outputs
の部分です。
inputs
にはワークフローへの入力が記述されています。ここでは参照配列の FASTA ファイルが唯一の入力項目です。
steps
の部分ではワークフロー内での各処理が記述されています。
bwa_index
create_dict
create_faidx
-
bind_fasta_index
step 1〜3についてはこの手の解析に慣れている方はだいたい想像がつくと思います。順に、BWA のインデックスファイルの生成、GATK で用いる .dict ファイルの生成、FASTA ファイルのインデックス (.fai) の生成をそれぞれ行っています。各項目中のrun
にはそのステップでの実際の処理内容を定義した CWL ファイルへの参照が、in
とout
には各ステップでの入出力がそれぞれ記載されています。このようにsteps
の部分では、複数の CWL ファイルに分かれて記述された処理をin
とout
を介して繋げることでワークフロー内の一連の処理内容が定義されいます。
なお、step 4 のbind_fasta_index
については他のワークフローと連結させるために便宜上含めたもので、1〜3 の出力結果を受け取ってひとまとめの結果として返すためのものです (FASTA とそのインデックスである .fai のように常に同じディレクトリに置いて使用する必要があるファイルを扱うようなケースに対応させるためです)。
ワークフローの結果ファイルはoutputs
の部分で定義されています。ここでは step 4 のbind_fasta_index
の結果を最終的な結果として出力しています。outputs
に記載の無いファイルは中間ファイルとして実行後に破棄されます。
[補足] ジョブの実行順について
上記 step 1〜3 をよく見るとこの3つは互いに独立している (i.e. あるステップの出力が別のステップの入力になっているわけではない) ことがわかります。このような場合 cwltool ではこの3ステップについては必ずしも記載順に実行されるわけではありません。それに対し、step 4 は step 1〜3 の出力結果を入力としているため、必ずこの3つのステップが完了した後に実行されます。
また、実行時に--parallel
オプションを指定すると step 1〜3 に関しては並列で実行されます。("[experimental]" という注意書きがなされていることと、ログが入り乱れてしまうことにご注意ください)
cwltool --singularity --parallel --outdir index ../workflows/prepare_reference.cwl --reference ref/IRGSP-1.0_genome_M_C_unanchored.fa
利用可能なワークフロー
(書きかけです)
DDBJのレポジトリで公開されているヒトゲノム解析パイプライン
Pitagora Network が提供するPitagora Workflow
次世代シークエンサーDRY解析教本 改訂第2版で紹介された解析手順を CWL 化した DAT2-cwl
その他
Singularityでは動かないワークフロー
Dockerfileからコンテナイメージをビルドする場合
Docker内部のファイル・ディレクトリを変更する処理が含まれるワークフロー
キャッシュの利用
--cachedir
オプションを指定する。
中間ファイルが保持される。
失敗時にresumeできる
udocker の利用
インストールは簡単。
実行は可能。Singularity を使うよりメモリが多く必要だった。
Dockerfileを使うようなケースでは Singularity と同様に失敗する。
遺伝研スパコンでは Singularity が使えるので udocker を使用する特別な必要はないと思われる。