シーケンシングのコアファシリティで働くゲノムインフォ系技術員が、Dryひとりでもぼちぼち業務を回せるように日ごろ頼りにしているツール等を紹介します。
Environment Modulesで簡単に環境変数の切り替え (+ツールの管理)
バイオインフォはとにかくツールが多く、その管理だけでも悩ましいです。
コアファシリティともなるとDe novo assemblyから*-Seqまで、様々な解析で使用するツールをザクザクとスパコンにインストールする必要があります。
以前はツールを所定のディレクトリ配下にインストール後、パスを通したり.bashrcに追記したりまたはパスの通ったディレクトリにソフトリンクを張ったり等したかと思うのですが、最近のスパコンやクラスタコンピュータでは環境設定の切り替えにModulesソフトウェアを使うことが増えている？ようです。
Modulesを使うと、簡単にバージョンの異なるツールやライブラリのパスを通したり引っぺがしたりできます。
基本の使い方
Modulesのインストールや基本の使い方については、詳しい記事があるのでそちらをどうぞ：Environment modulesのインストール
Modulesがインストールされていると、moduleコマンドでインストールされているツールやライブラリの確認、環境変数のロード/アンロードなどができます。
たとえばツールのインストールに cmake 3.8以上、Python 3.6以上が必要な場合、下記のようにロードすればOKです。
# 現在のcmakeとpythonのバージョン
$ cmake --version && python3 --version
cmake version 2.8.12.2
Python 3.5.0
# モジュールをロード
$ module load cmake/3.11.4 python/3.7.3
# cmake 3.11.4 と Python 3.7.3 が使えるようになった
$ cmake --version && python3 --version
cmake version 3.11.4
CMake suite maintained and supported by Kitware (kitware.com/cmake).
Python 3.7.3
また、うちではOxford Nanoporeのシーケンサーがあるのですが、Nanoporeのツールはベースコーラーをはじめとしてアップデートが頻繁にあり、ツールの管理も大変なことになってきました。
でもModulesを使うとこのように、インストールされているツールとそのバージョンの一覧が得られます。嬉しい！
$ module avail
----------------------------- /path/to/modules/modulefiles/ -----------------------------
bonito/0.0.1 guppy/3.2.1 guppy/3.3.3 megalodon/1.0.0
guppy/2.3.7 guppy/3.2.2 guppy/3.4.3 pigz/2.4
guppy/3.0.3 guppy/3.2.4 medaka/0.11.1 seqkit/0.11.0
guppy/3.1.5 guppy/3.3.0 medaka/0.6.5 tree/1.6.0
Modulefileの作成と追加
ツールの追加は、まずはツールをインストールしたのち、modulefile(設定ファイル)にパスなどの情報を書きます。
このmodulefileを$MODULEPATHで指定されたパス以下におけばOK。
modulefile - 普通(?)のツールの場合
/path/to/modulefiles以下にツール名のディレクトリを作成し、その下にバージョン名で保存した設定ファイルを置きます。
treeコマンドをスパコンにインストールした際のmodulefileの中身はこんな感じ：
#%Module1.0
set appname [lrange [split [module-info name] {/}] 0 0]
set appversion [lrange [split [module-info name] {/}] 1 1]
set apphome /path/to/modules/$appname/$appversion
prepend-path PATH $apphome/bin
prepend-path LD_LIBRARY_PATH $apphome/lib
prepend-path MANPATH $apphome/share/man
$PATHだけでなく$MANPATHなどの環境変数も宣言できます。
modulefileは Tcl(Tool Command Language) という言語で書くのだそうです。
modulefile - Pythonのvenvを使用するツールの場合
このごろは仮想環境を作成して解析を行うツールも増えてきました。NanoporeのポリッシングツールMedakaや、ベースコーラーのBonitoなどがそれです。
PythonのvenvをModulesでロード/アンロードする場合は、ifを使って下記のような書き方ができます。activate/deactivate をそれぞれ load/unload に対応させています。
NanoporeのBonitoというベースコーラーをインストールした際のmodulefileはこんな感じ:
#%Module1.0
set appname [lrange [split [module-info name] {/}] 0 0]
set appversion [lrange [split [module-info name] {/}] 1 1]
if { [ module-info mode load ] } {
puts stdout "source /path/to/modules/bonito/0.0.1/venv3/bin/activate;"
} elseif { [ module-info mode remove ] } {
puts stdout "deactivate;"
}
すると、下記のように load/unload でvenvの activate/deactivate ができます。嬉しい！
$ module load bonito/0.0.1
(venv3) $ module list
Currently Loaded Modulefiles:
1) /bonito/0.0.1
(venv3) $ module unload bonito/0.0.1
$ module list
No Modulefiles Currently Loaded.
Modulesを使うとツールのバージョン切り替えはもとより、「あのツールインストールしてたっけ？」というときもmodule availするだけなのでとても楽になります。
(おまけ) インストールしたツールをラボメンバーへおしらせ
ラボメンバーにインストールしたツールの情報を共有する用に、SlackのWorkflowを作りました。
私以外は皆Wetなので、「何をするツール？」「うちのサービスにどう使えそう？」といった情報が入ったメッセージをポストするようにしています。
Workflowを起動すると下記のようなフォームが現れるので、質問に答えてSubmitすれば、
1. ツール名
2. 何をするツールか
3. どの計算サーバにインストールされているか
4. コアファシリティのサービスにどのように利用されうるか
といった情報とともに、ラボのSlackに新規ツールインストールのお知らせがポストされます。かしこいねSlack!
ガッと作ったのでお直しする予定。
seqkit - FASTA/Qの前処理・後処理を楽にするtoolkit
FASTA/Qの前処理、後処理・・・以前はcut, grep, sed, awkなどを組み合わせたワンライナーを書いたり、書き捨てのスクリプトで対応していました。
seqkitにはそれらの大抵の操作ができるコマンドがそろっています。
例えばこんな感じ：
Usage:
seqkit [command]
Available Commands:
amplicon retrieve amplicon (or specific region around it) via primer(s)
bam monitoring and online histograms of BAM record features
common find common sequences of multiple files by id/name/sequence
concat concatenate sequences with same ID from multiple files
convert convert FASTQ quality encoding between Sanger, Solexa and Illumina
duplicate duplicate sequences N times
faidx create FASTA index file and extract subsequence
fish look for short sequences in larger sequences using local alignment
... ...
私はよくアセンブル前のリードのランダムサンプリングにsampleを使ったり、アセンブル結果をstatsで確認したりします。
さらに凝ったことも出来て、プロット用にGC含量をタブ区切りテキストで出力することもできます。
（下記は公式からの引用）
$ zcat hairpin.fa.gz \
| seqkit sliding -s 5 -W 30 \
| seqkit fx2tab -n -g
cel-let-7_sliding:1-30 50.00
cel-let-7_sliding:6-35 46.67
cel-let-7_sliding:11-40 43.33
cel-let-7_sliding:16-45 36.67
cel-let-7_sliding:21-50 33.33
cel-let-7_sliding:26-55 40.00
...
現在(v0.11.0)、35個のコマンドがあります。最近コマンドが増えてさらに多くのことが出来るようになりました。
seqkit versionでアップデートの確認もできるのがすごく嬉しいです。
Nextflow - パイプライン構築のお供に
内部向けの解析サービスとしてNanoporeデータのDe novo assemblyを提供しているのですが、そこで使う解析パイプラインをNextflowで構築しました。NextflowはGroovyという言語を拡張したものなんだそうです。
もともとはBashスクリプトで書いていたのですが、若干トリッキーな書き方が必要だったり、2週間前の自分のコードが読めなくなったりと、自分のBashスクリプティング力の低さで日々の生活が辛くなり、ワークフローシステムに手を出しました。
そもそもNextflowなどのワークフローシステムはデータ解析の再現性や自動化、共有といった点を重視しているのかなと思うのですが、私の場合は「こんなことを楽に書けないかな〜」という思いつきを（Bashスクリプトよりは）簡単に実装できたので今のところNextflowでパイプラインを構築しています。
Nextflowについては他の方がアドベントで書いてくださる予定なので、ここでは私のお気に入りをいくつか紹介したいと思います。
個人的に気に入っているのは、パイプライン実行時にエラーが返された場合にジョブを再投入する際に使う errorStrategy, maxRetriesというオプションです。
設定ファイルで下記のように書きます：
process {
executor = 'slurm'
queue = 'compute'
errorStrategy = { task.exitStatus in [1, 137..143] ? 'retry' : 'terminate' }
maxRetries = 2
}
この場合、exit status 1, 137...143 が返された場合は再度ジョブを投入（最大2回）、それ以外の場合はジョブを終了します。
また、withNameでプロセスごとに異なるリソースを割り当てることができ、かつ前述したmaxRetriesとtask.attemptを組み合わせて、リトライのたびにメモリを増やしていくこともできます。
process {
executor = 'slurm'
queue = 'compute'
errorStrategy = { task.exitStatus in [1, 137..143] ? 'retry' : 'terminate' }
maxRetries = 2
# アセンブルでコケたらメモリを100*retry数にして再実行
withName: 'assembly' {
queue = 'largemem'
memory = { 100.GB * task.attempt }
}
# ポリッシングでコケたらメモリを30*retry数にして再実行
withName: 'polishing' {
queue = 'compute'
memory = { 30.GB * task.attempt }
}
}
ほかにも「なんかこんな機能ありそうだな…」と思ってドキュメントを読むとできたりするので、当分はNextflowに助けてもらおうと思います。
Nextflow歴はまだ半年なので、ここに載せたサンプルコードより良い書き方があれば教えてもらえると嬉しいです。
以上ざっくりですが、シーケンシングのコアファシリティで働くゲノムインフォ系技術員から、仕事を楽にしてくれるツール等の紹介でした。