続) 遺伝研スパコンで Singularity とか cwl を使う上での tips

はじめに

分量が多くなってきたので分割しました。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ディレクトリに含まれている動作確認用のファイルを用います。

1. 参照ゲノム配列の取得

   同梱のスクリプトファイルを使って 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 ファイルだけです。残りの２つは snpEff を使用して変異のアノテーションを行うときに使用するためのものなので今回は必要ありません。

2. 動作確認用データについて

   $ 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は次項以下で説明する｀ジョブファイル｀と呼ばれるもので、今回は使用しません。

3. 参照ゲノム配列の 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は不要になります)。

4. リード配列のアラインメント

   続いて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の部分ではワークフロー内での各処理が記述されています。
1. bwa_index
2. create_dict
3. create_faidx
4. 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 をよく見るとこの３つは互いに独立している (i.e. あるステップの出力が別のステップの入力になっているわけではない) ことがわかります。このような場合 cwltool ではこの３ステップについては必ずしも記載順に実行されるわけではありません。それに対し、step 4 は step 1〜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
次世代シークエンサーＤＲＹ解析教本　改訂第２版で紹介された解析手順を CWL 化した DAT2-cwl

その他

Singularityでは動かないワークフロー

Dockerfileからコンテナイメージをビルドする場合
Docker内部のファイル・ディレクトリを変更する処理が含まれるワークフロー

キャッシュの利用

--cachedirオプションを指定する。
中間ファイルが保持される。
失敗時にresumeできる

udocker の利用

インストールは簡単。
実行は可能。Singularity を使うよりメモリが多く必要だった。
Dockerfileを使うようなケースでは Singularity と同様に失敗する。
遺伝研スパコンでは Singularity が使えるので udocker を使用する特別な必要はないと思われる。