Help us understand the problem. What is going on with this article?

続) 遺伝研スパコンで 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 系の両方がインストールされていますがこれを用いるよりも、pyenvanyenv を使ってホームディレクトリに好きなバージョンを指定してインストールしてしまうのが良いと思います。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 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を実行します。

以下ではtestディレクトリに含まれている動作確認用のファイルを使って動かしてみます。

$ cd test/
$ ls
README.md  download_reference_files.sh  rapdb-pipeline.test.job.yaml  read1.fq.gz  read2.fq.gz

まずは、同梱のスクリプトファイルを使って RAP-DB から必要な参照データを取得しrefというディレクトリに保存します。

$ ./download_reference_files.sh
$ ls -lh ref/
total 1139976
-rw-r--r--  1 nigyta  staff   365M 12 12 09:30 IRGSP-1.0_genome_M_C_unanchored.fa
-rw-r--r--  1 nigyta  staff   142M 12 12 09:30 RAP-DB_MSU_concatenated_for_snpEff.gtf
-rw-r--r--  1 nigyta  staff    46M 12 12 09:30 RAP-DB_MSU_concatenated_protein.fa

使用するのは一番上の参照ゲノム配列の FASTA ファイルだけです。残りの2つは snpEff を使用して変異のアノテーションを行うときに使用するためのものなので今回は必要ありません。

実行その1

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 という名称で示されるオプションは、'ジョブファイル' と呼ばれるワークフローへの入力ファイルやパラメータを記述したファイルを指定するためのものです(ジョブの実行順序?と勘違いしてしまったが、おそらくワークフローへの入力を指示するためのファイルという意味だと思われる)。ジョブファイルを用いて実行する方法は実行その2で紹介するので、ここではコマンドラインから引数を与える方法で実行します。必須の引数は参照ゲノムの FASTA ファイルを指定する--referenceだけです。

cwltool --singularity --outdir index ../workflows/prepare_reference.cwl --reference ref/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は不要になります)。

実行その2

続いて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

重要なのはinputsstepsoutputsの部分です。
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 ファイルへの参照が、inoutには各ステップでの入出力がそれぞれ記載されています。このようにstepsの部分では、複数の CWL ファイルに分かれて記述された処理をinoutを介して繋げることでワークフロー内の一連の処理内容が定義されいます。
なお、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 を使用する特別な必要はないと思われる。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした