多数の配列についてBLAST検索をする例

本ページでは、多数の配列を、並列にblast検索する例を紹介する

まず、多数の配列が含まれるFASTA形式のファイル(例Trinity.fasta)を分割して適当な本数の配列に別れているファイルを用意します。 

#!/usr/bin/ruby
require 'bio'

ff = Bio::FlatFile.open(nil, ARGV[0])
maxcount = (ARGV[1].to_i or 1000)
maxdir = 100
filecount=0
dircount = 0
count=maxcount
owd = Dir.getwd
while fe = ff.next_entry
  if filecount % maxdir == 0 && count % maxcount == 0
    Dir.chdir(owd)
    dircount += 1
    newdir = (ARGV[0] + "." + dircount.to_s)
    Dir.mkdir(newdir)
    Dir.chdir(newdir)
  end
  if count == maxcount then
    filecount += 1
    of=open(ARGV[0]+"."+filecount.to_s,"w")
    count = 0
  end
  of.print fe.to_s
  count += 1
end

のようなプログラムを利用すれば１つのファイルに指定の本数(無指定時1000本)のファイルが含まれ、 各directoryに100個以内のファイルがあるようにサブディレクトリーが作られます。利用方法は例えば 

$ ruby flatsplit Trinity.fasta 100

のようにするとそれぞれのファイルに100本ずつずつのファイルに分割されます。 このファイル１つをqueryにblastx検索するコマンドが例えば、 

$ /bio/bin/blastx -query Trinity.fasta.1 \
        -db TAIR10_pep_20110103_representative_gene_model_updated \
        -num_alignments 10 -num_descriptions 10 -outfmt 7 \
        -out Trinity.1.blastx-Ath

のとき, SGEのarrayジョブにして1番の部分を変数にすると

#!/bin/bash
#$ -cwd
#$ -S /bin/bash
#$ -q normal
#$ -t 1:100

n=$SGE_TASK_ID
export BLASTDB=/home/nasu

/bio/bin/blastx -query Trinity.fasta.$n \
	-db TAIR10_pep_20110103_representative_gene_model_updated \
	-num_alignments 10 -num_descriptions 10 -outfmt 7 \
	-out Trinity.$n.blastx-Ath

のようなjobscriptができて、これを

$ qsub blastx1.job

などとsubmitすることで100個の入力ファイルを処理する事ができます。

これを、さらに多数のディレクトリについて自動的に実行するのに、上記のようなJobを 自動生成してqsubを発行すれば、一命令で全部できる事になります。

例えばrubyで書くなら、

#!/usr/bin/env ruby

owd = Dir.getwd
(1..11).each do |i|
  jobstart = 1 + 100 * (i-1)
  jobend = i*100
  Dir.chdir(owd)
  Dir.chdir("Trinity.fasta.#{i}")
#  jobpipe=IO.popen("cat", "w") do |jobpipe|
  jobpipe=IO.popen("qsub -N blast#{i}", "w") do |jobpipe|
    jobpipe.puts <<EOS
#!/bin/bash
#\$ -cwd
#\$ -S /bin/bash
#\$ -q normal
#\$ -t #{jobstart}:#{jobend}

n=$SGE_TASK_ID
export BLASTDB=/home/nasu

/bio/bin/blastx -query Trinity.fasta.$n \
	-db TAIR10_pep_20110103_representative_gene_model_updated \
	-num_alignments 10 -num_descriptions 10 -outfmt 7 \
	-out Podo1.$n.blastx-Ath
EOS
  end
end

のようにして各ディレクトリーでqsubを発行する事ができる。

この例では、job scriptにもとづいてrubyのプログラムを作成しているが、 jobの中でdirectory毎に異なる要素をあらかじめ計算し、 here documentに#{jobstart}:#{jobend}などとして埋め込むことによって、 jobの文字列を生成している。

なお、上の例では、jobはファイルに書かずに、直接qsubコマンドにパイプで 渡している。

各ファイルに100本ずつ、各ディレクトリーに100ファイルずつで、ディレクトリーあたり１万本処理に相当する。 上の例では、10万本強で11個のディレクトリーに分けてあるので、(1..11).eachとしているが、ここは本数に依存して返る必要がある。

発展としては、より一般化して、総配列数や総ファイル数にもとづいて、jobを生成する回数や、最後のjobの終わりをより 精度良く制御する事も可能である。

ただし、１回実行すれば良いと言う意味では上記で十分であり、そこまでの 一般化をする必要があるわけではないので、ここでは上記までを示して終わりとする。