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

# はじめに
分量が多くなってきたので分割しました。Singularityについては[遺伝研スパコンで Singularity とか CWL を使う上での tips]()をご参照ください。ここではCWLについて書きます。

CWL (Common Workflow Language)　は名前の通り解析や定型作業のワークフローの手順を記述したもので、**仮想コンテナを利用することで実行環境によらず再現性のある結果が得られる**ことが大きなメリットだと考えています。

CWL についてより詳しく知りたい場合には[CWL Start Guide JP](https://github.com/pitagora-galaxy/cwl/wiki/CWL-Start-Guide-JP) を、CWL を書いてみたいという人には@tm_tn氏の[雑に始める CWL！](https://qiita.com/tm_tn/items/4956f5ca523f7f49f386)をおすすめします。特に、

- 自分が作成したワークフローを他の人にもつかってもらいたい
- クラウドを含めた様々な環境でワークフローを実行したい

といった目的には合致すると思います。また、しばらく時間がたってから解析をやり直そうとしたらいつの間にか実行環境が変わっていてツールが動かない、といった悩みからも解放されるはずです。

ここでは、遺伝研スパコンを使って CWL を実行するための環境構築および実際にワークフローを実行する手順を紹介します。

# 実行環境の構築
遺伝研スパコンで CWL で書かれたワークフローを実行するには下記が必要です。

- Python
次に挙げる cwltool をインストール・実行するために必要です。公式サイトによると Python 2.7 もサポートしているようですが、Python 3.5 以降をおすすめします。遺伝研スパコンには `/usr/local/pkg/python` 以下に Python 2.7 系と 3 系の両方がインストールされていますがこれを用いるよりも、[pyenv](https://github.com/pyenv/pyenv) や [anyenv](https://github.com/anyenv/anyenv) を使ってホームディレクトリに好きなバージョンを指定してインストールしてしまうのが良いと思います。[Anaconda](https://www.anaconda.com) or [Miniconda](https://docs.conda.io/en/latest/miniconda.html) でも可です。これらのいずれも使用したことがないのであれば、Miniconda を導入しておくと後の手順が簡単です。
以下の手順では`pip`または`conda`コマンドが使用できることを前提とします。

- [cwltool](https://github.com/common-workflow-language/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](https://github.com/nodenv/nodenv) を用いるか、Anaconda/Miniconda を使用している場合には`conda`コマンドを使用するのが簡単です。`conda`コマンドを使用する場合のインストール方法は以下の通りです。なお、遺伝研スパコン以外の環境で Docker を使って CWL を実行する場合には、実行時に Node.js の Docker イメージを自動的に取得して用いるためインストールは不要です。

```
conda install nodejs
```
 
# 実行してみる
例として CWL で書かれた[イネの変異解析のためのワークフロー](https://github.com/nigyta/rice_reseq)を実行してみます。これは、Rice Annotation Project Database (RAP-DB) で公開されている[標準解析パイプライン](https://rapdb.dna.affrc.go.jp/genome-wide_variations/Analysis_workflow_for_detection_of_genome-wide_var.html)を 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](https://github.com/nigyta/rice_reseq/blob/master/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のレポジトリで公開されている[ヒトゲノム解析パイプライン](https://github.com/ddbj/human-reseq)
Pitagora Network が提供する[Pitagora Workflow](https://github.com/pitagora-network/pitagora-cwl)
[次世代シークエンサーＤＲＹ解析教本　改訂第２版](https://www.hanmoto.com/bd/isbn/9784780909838)で紹介された解析手順を CWL 化した [DAT2-cwl](https://github.com/pitagora-network/DAT2-cwl) 


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

### キャッシュの利用
`--cachedir`オプションを指定する。
中間ファイルが保持される。
失敗時にresumeできる

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