Rabbit Ears

gmapとは

gmapはTopHatなどのshort read処理ツールをSun Grid Engineで並列分散処理させる際の支援ツールです．gmapを使うことで，各ツールのオプションの設定や出力ディレクトリの作成，Sun Grid Engineへのジョブ投入といった一連の処理がコマンドひとつで簡単に実行できるようになります．TopHat，Bowtie，SOAP2の実行に対応しています．

現在gmapはGitHubで公開中です．以下のURL先から入手することができます．利用にはRuby(1.9.1以降)とSun Grid Engine，およびTopHat，Bowtie，SOAP2などのツールがインストールされているクラスタ環境が必要です．
https://github.com/mickey24/gmap

以降，gmapを使った場合のメリットについて，例を交えながら紹介していきたいと思います．

TopHatをコマンドラインから使う場合

以下の設定でTopHatによる解析を行う場合について考えます．

• マッピング先のゲノムのBowtieインデックス：/path/to/genome/mouse/chr1 など (染色体番号：1〜19, X, Y)
• 入力ファイル：/path/to/SRP0123/SRX0123/SRR0123.fastq
• 出力ディレクトリ：./result/1 など(各染色体ごとに結果を保存)

TopHatをコマンドラインから手入力で実行した場合の例を以下に示します．

$ mkdir result
$ tophat -o result/1 /path/to/genome/mouse/chr1 \
  /path/to/SRP0123/SRX0123/SRR0123.fastq && \
tophat -o result/2 /path/to/genome/mouse/chr2 \
  /path/to/SRP0123/SRX0123/SRR0123.fastq && \
tophat -o result/3 /path/to/genome/mouse/chr3 \
  /path/to/SRP0123/SRX0123/SRR0123.fastq && \
...
tophat -o result/19 /path/to/genome/mouse/chr19 \
  /path/to/SRP0123/SRX0123/SRR0123.fastq && \
tophat -o result/X /path/to/genome/mouse/chrX \
  /path/to/SRP0123/SRX0123/SRR0123.fastq && \
tophat -o result/Y /path/to/genome/mouse/chrY \
  /path/to/SRP0123/SRX0123/SRR0123.fastq

これは流石に面倒ですね．同じようなコマンドを染色体の数(1〜19, X, Yで計21)だけ入力しなければならず，またミスもしやすく危険です．for文を使えば以下のようにある程度短く効率的に書けますが，それでもTopHatを使う度に毎回自分の手でfor文を入力するのは大変です．

$ mkdir result
$ for i in `seq 1 19` X Y
do
  tophat -o result/${i} /path/to/genome/mouse/chr${i} \
    /path/to/SRP0123/SRX0123/SRR0123.fastq
done

シェルスクリプトを使う場合

シェルスクリプトは同じようなコマンドを何度も実行する時に使うと非常に便利です．例えば以下のようなシェルスクリプトを書き，mouse.shなどのファイル名で保存します．

input=$1
output=$2

mkdir ${output}

for i in `seq 1 19` X Y
do
  tophat -o ${output}/${i} /path/to/genome/mouse/chr${i} ${input}
done

これを以下のようにコマンドラインから実行します．

$ sh mouse.sh /path/to/SRP0123/SRX0123/SRR0123.fastq result

これでmouseのゲノムに対するTopHatの処理が随分楽に実行できるようになりました．しかし，この場合でも以下のような問題が残ります．

1. ゲノムのBowtieインデックスやTopHatのオプションを変えたい時にスクリプトを書き換えなければならない
2. TopHatを染色体の数だけ逐次実行(上の場合は21回)するのは非常に時間がかかる

1についてはシェルスクリプトのコマンドライン引数でゲノムなどを指定するようにすればある程度解決できますが，今度はスクリプト実行時に毎回Bowtieインデックスのpathのような冗長な引数を入力しなければならなくなるという問題が発生してしまいます．染色体番号の指定も面倒になりそうです．かといってゲノムやオプションごとにスクリプトをコピーして当該箇所だけ書き換えたものを作っていくと，同じようなスクリプトが散在して管理が大変になる怖れがあります．

2については，ひとつのマシンで全染色体に対する処理を行うのではなく，以下のようにクラスタ環境で並列分散処理をすることで解決できます．

Sun Grid Engineを使う場合

Sun Grid Engineはクラスタ環境用のジョブ投入システムです．qsubコマンドを使うことで，シェルスクリプトなどのプログラムをジョブとしてクラスタ環境に投入し，Sun Grid Engineが割り当てたマシンやリソースを使って実行させることができます．

Sun Grid Engineを使った場合のメリットは以下の通りです．

1. 大量のジョブを複数のマシンに投入し，同時に実行させることができる
2. ジョブのスケジューリングやリソース割り当てといった問題はSun Grid Engineが解決してくれる

TopHat，Bowtie，SOAP2は，ひとつの入力ファイルに対する染色体ごとの処理は基本的にそれぞれ独立して実行させることができます．そのため，Sun Grid Engineを使って各染色体ごとに1ジョブずつ投入することで同時に実行させることが可能です．十分な計算リソースがあれば，ひとつのマシンで逐次実行する場合と比べて大幅な高速化が可能です．

しかし，Sun Grid Engineを使うことで新たな問題もいくつか出てきます．

1. 新たにqsubコマンドの引数も適切に指定する必要がある．シェルスクリプトの引数でこれらを指定できるようにした場合，スクリプトのコマンドライン引数の数がさらに増え，スクリプトのコードも複雑化する．
2. PATHなどの環境変数も気を遣う必要がある
3. ツールの引数や環境変数が不適切だった場合，投入されたジョブが異常終了することがある．この場合(特に環境変数が原因の場合)エラーの特定が難しい．

TopHatなどのツールを実行する度に，これらのことに気を遣いながらシェルスクリプトを書くのは非常に大変です．また，前節で挙げた問題1もまだ解決していません．そこでgmapの出番です．

gmapを使う場合

gmapを使った場合のコマンドを以下に示します．

$ gmap -t tophat -g mouse -o result \
  /path/to/SRP0123/SRX0123/SRR0123.fastq

これだけで，mouseの各染色体についてTopHatのジョブをSun Grid Engineで同時実行することができます．各ツールの引数の指定，Sun Grid Engineへのジョブ投入，出力ディレクトリの作成といった作業は全部gmapがやってくれます．

gmapを使った場合の利点を以下に示します．

1. ツールやゲノムの指定が簡単にできる
2. 設定ファイル(.gmap)を元にgmapがSun Grid Engineや各ツールの引数を指定してくれる
3. ジョブ投入や出力ディレクトリの作成といった作業はgmapがやってくれる
4. gmapが設定や与えられた引数が適切かどうかチェックし，ジョブが異常終了することを未然に防いでくれる

gmapを使うには，あらかじめホームディレクトリに.gmapという名前のファイルを設置し，そこに設定を書いておく必要があります．.gmapファイルには主に以下の設定を記述します．

• TopHatなど各ツールのpath
• genome_config: 各ゲノムのBowtieインデックスのpath，染色体番号
• project_config: TopHatなどのツールのオプション(SRPごとに指定)

各ゲノムのBowtieインデックスについての設定例を以下に示します．

genome_config:
 mouse:
  genome_path : /path/to/genome/mouse/chr%s
  chrnum      : 1-19,X,Y
 human:
  genome_path : /path/to/genome/human/chr%s
  chrnum      : 1-22,X,Y

このように設定を記述しておくことで，違うゲノム(例: human)を使ってTopHatを処理したい時も簡単に実行できるようになります．

$ gmap -t tophat -g human -o result \
  /path/to/SRP0123/SRX0123/SRR0123.fastq

tophatなどのツールに与える引数も.gmapファイルに書きます．入力ファイルのSRPの番号ごとに同じ引数を使うことが多いため，各ツールの引数はSRPごとに指定できるようになっています．以下に例を示します．

project_config:
 default:
  tophat    : "--solexa-quals -r 200 --mate-std-dev 50"
  bowtie    : "-X 250 -I 150"
  soap2     : "-x 250 -m 150"
 SRP0123:
  tophat    : "--solexa-quals -r 200 --mate-std-dev 50"
  bowtie    : "-X 250 -I 150"
  soap2     : "-x 250 -m 150"

これらの設定を書いておくだけで，入力ファイルに合わせてTopHatなどのツールに適切な引数が与えられるようになります．

以上

gmapでTopHat，Bowtie，SOAP2を使う大規模計算を少しでもサポートできれば幸いです．