CommonWL for practical use 1

作成したCWL ツール・ワークフローを実際に動かしてみた(その1)

次世代シークエンサーDRY解析教本 改訂第2版(以下、DRY解析教本2)の目次が出版社のページから公開された。 その目次を見ての通り、Level1(準備編)、Level2(実践編)、Level3(応用編)の三編構成で、第1版と大きな構成に変化はない。 大きくは、Level2にメタゲノム解析、アッセンブル解析が加わって、Level3が総取っ替えされたといえば適切だろうか。 ちなみに、この本のtwitterのハッシュタグは #NGS_DAT で、この本に関わるtweetは、このハッシュタグを付けて是非。

しかしながら、大きな変更点が実はそれ以外にある。 Common Workflow Language (CWL)の追加だ。 Level2の一番最後(p331-338)に付録「CWLがあれば、DRY解析はもう怖くない」がある。 CWLについて日本語でコンパクトに解説してくれて、非常に有用である。 しかし、実はすごいのはここだけでない。 Level2で紹介した全てのワークフローが、CWLで書かれたWorkflowとしてこの記事のサポートサイトで公開されている。 つまり、ここにあるCWLをgit cloneしてきて自分にとって必要なように書き換えさえすれば、CWLを使ったNGSデータ解析が割と手間が少なくできる。

残念ながら、HISAT2-Stringtieによるワークフローは、DRY解析教本2には収録されなかった。 そこで、このワークフローをCWLで記述して実際に動かすことを目標として、2019年12月16日(月)に行われた23rd Workflow Meetupにて、自分の課題として取り組んだ。 Meetup当日のlive codingに関しては別に書かれている(「実際に、指先一つで立ち上げる CWL ツール・ワークフロー作成環境を試してみた」)ので、そちらを参照。 このエントリは、その後の顛末を。

まず、live codingで書いてもらったのはHISAT2のmappingの部分からなので、その前のindexを作成する必要がある。 もっともメジャーな生物だとすでにindexを直接ダウンロードできたりするが、いつもそうであるとは限らない。 FASTA形式のゲノム配列ファイルをゲットしてそれに対してindexを作成する。 幸い、pitagora-networkプロジェクトでHISAT2のCWLを作成し公開してくれているので、それをgit cloneしてありがたく利用させていただく。

1
2
3
4
#!/bin/bash
git clone https://github.com/pitagora-network/pitagora-cwl
cwltool pitagora-cwl/tools/hisat2/index/hisat2_index.cwl \
--index_basename org1 --reference_fasta organism1.fa

無事、indexができたので、先日のlive codingで作成していただいたrepoをgit cloneしてくる。 そして、それをつかってHISAT2のmappingを実行する。

1
2
3
4
5
#!/bin/bash
git clone https://github.com/tom-tan/cwl-live-coding-2019-12-16
cwltool cwl-live-coding-2019-12-16/cwl/hisat2-stringtie_wf_se.cwl \
--annotation organism1.gtf --fastq fibroblast.fq.gz \
--hisat2_idx_basedir . --hisat2_idx_basename org1 --nthreads 6

実行中に、Samtools view が並列化されていないことに気が付いたので、pitagora-cwl/tools/samtools/sam2bam/samtools_sam2bam.cwlを以下のように改変して並列化可能なように。

cwlVersion: v1.0
class: CommandLineTool

hints:
  DockerRequirement:
    dockerPull: quay.io/biocontainers/samtools:1.9--h46bd0b3_0

baseCommand: ["samtools", "view", "-b"]

arguments:
  - prefix: -o
    valueFrom: $(runtime.outdir)/$(inputs.output_filename)

inputs:
  output_filename:
    label: "Filename to write output to"
    doc: "Filename to write output to"
    type: string
    default: out.bam
  nthreads:
    type: int
    default: 1
    inputBinding:
      prefix: -@
  input_sam:
    label: "SAM format input file"
    doc: "SAM format input file"
    type: File
    inputBinding:
      position: 50

outputs:
  bamfile:
    type: File
    outputBinding:
      glob: $(inputs.output_filename)

$namespaces:
  s: https://schema.org/
  edam: http://edamontology.org/

s:license: https://spdx.org/licenses/Apache-2.0
s:codeRepository: https://github.com/pitagora-network/pitagora-cwl
s:author:
  - class: s:Person
    s:identifier: https://orcid.org/0000-0003-3777-5945
    s:email: mailto:inutano@gmail.com
    s:name: Tazro Ohta

$schemas:
  - https://schema.org/docs/schema_org_rdfa.html
  - http://edamontology.org/EDAM_1.18.owl

その上で、cwl-live-coding-2019-12-16/cwl/hisat2-stringtie_wf_se.cwlsamtools_sam2bamの部分をローカルなCWLを参照するように、

  samtools_sam2bam:
    run: samtools_sam2bam.cwl

に書き換えて実行すると、sam2bam(samtools view)も並列で動くようになった。

これを使っていくつかのFASTQに対して実行し、複数のGTFを作成することができた。 そこで調子に乗って、stringtie mergeでそれらのGTFを束ねるプログラムもCWLで、ともう一つ作成してもらったCWLを動かそうとした。しかし、

1
2
3
4
#!/bin/bash
cwltool cwl-live-coding-2019-12-16/cwl/stringtie-merge.cwl \
--nthreads 6 --guide organism1.gtf --o_name merged.gtf \
--mergelist_txt mergelist.txt

とすると Error: cannot find file 'Brain.gtf' and /var/lib/cwl/stg4116466b-c5a7-4b66-8aac-76fe5a0ef6d1/mergelist.txt does not look like GFF! というエラーとなった。 どうも、テキストファイルの中に書いたファイルというのがうまくCWL内から見られていないようで。 今日のところはここで力尽きました(続く)。


Written by Hidemasa Bono in NGS_DAT on 木 19 12月 2019.