ABySSむけデータ準備

ABySS向けのデータとしてはIllumina qseq formatでもfastq formatでも構わない。しかしながら、アラインメントのフェーズで各リードのread1とread2がセットで存在しないと非常に大量のメモリーを必要とする。このため、read1とread2がセットになるようにread1とread2のファイルをあらかじめマージすることが望ましい。

qseq.txt formatの場合

Illuminaからの出力はread1, read2が分かれたqseq.txt format とそれをfastq formatにしたものが多く使われる。

qseq.txt formatはs_[レーン番号]_[リード番号(1 or 2)]_[４桁数字]_qseq.txt というようなファイル多数になっている。

1行で1 readのデータである read1とread2のデータは同様に並んでいるものとすると、 read1のファイルから1行、read2のファイルから1行読んで出力する事を繰り返せば、同じ断片のread1とread2が隣に来るファイルができる。例えば以下のようなプログラムmerge12.rbを作成し、

#!/usr/local/bin/ruby
require 'fileutils'
from_dir = ['read1', 'read2']
to_dir = 'read12'
cur_dir=Dir.pwd
FileUtils.mkpath(to_dir)
Dir.chdir(from_dir[0])
files=Dir.glob("*").sort
Dir.chdir(cur_dir)
#p files

files.each do |filename|
  pairfilename = filename.sub(/(s_\d+)_1_/,'\1_2_')
  newfilename = filename.sub(/(s_\d+)_1_/,'\1_')
#      p filename
#      p pairfilename
  open("#{from_dir[0]}/#{filename}") do |if1|
    open("#{from_dir[1]}/#{pairfilename}") do |if2|
      open("#{to_dir}/#{newfilename}", "w") do |of|
        if1.each_line do |line|
          of.puts line
          of.puts if2.gets
        end
      end
    end
  end
end

$ echo ruby merge12.rb | qsub -cwd

ここで、もちろんruby merge12.rbで実行する事もできるが、 "echo"と "|qsub -cwd"ではさむことにより、SGEのjobとして適当なノードで実行させる事を指示している。これにより、現在ログインしているコンピュータのIO負荷が増大するのを避ける。もちろん、nfsサーバーにはしっかり負荷がかかります。

これで、read12というディレクトリにバラバラのファイルが書かれます。このさき、pass filerのreadを選択して単一ファイルにまとめるには

cat read12/* | grep 1$ > ../LIBqf.qseq

fastq formatの場合

fastq formatは基本的に4行で1 readのデータを表すフォーマットであるから、上と同様に交互に４行ずつよんで書き出せば良いかのようにも思われるが、コメント等に対する頑健性があまりにも欠けるやり方なので、区切りを認識して処理するようにする。

このため、biorubyのライブラリーを用いることとする。すると、システムにインストールされているrubyでは自分でいれたライブラリーを使えないこと、ruby-1.8よりruby-1.9の方が速いことを踏まえ、rubyからインストールする。

$ wget http://koke.asrc.kanazawa-u.ac.jp/lpm/repository/ruby.lpm
$ lpm install ruby.lpm
$ wget http://koke.asrc.kanazawa-u.ac.jp/lpm/repository/bioruby.lpm
$ lpm install bioruby.lpm

若干美しくないが、以下のようなプログラムを保存。（例えば~/rubysrc/fastqmerge.rbとして)

#!/bin/env ruby

require 'bio'

ff1 = Bio::FlatFile.open(nil, ARGV[0])
ff2 = Bio::FlatFile.open(nil, ARGV[1])

ff1.each_entry do |fe1|
  fe2 = ff2.next_entry
  puts "@#{fe1.definition}\n#{fe1.sequence_string}\n+#{fe1.definition}\n#{fe1.quality_string}\n"
  puts "@#{fe2.definition}\n#{fe2.sequence_string}\n+#{fe2.definition}\n#{fe2.quality_string}\n"
end

$ ruby ~/rubysrc/fastqmerge.rb /path/to/sequence/files/s_8_1_2_sequence.txt /path/to/sequence/files/s_8_2_2_sequence.txt > s_8_2.sequence.txt

このようにして、ライブラリーごとにread1, read2をマージします。 8GBほどのファイルを二つ結合するのに45分 (約3 MB/s)

入力ファイルがgzip等で圧縮されている場合には、openの所を工夫し、

#!/bin/env ruby

require 'bio'
fh1 = IO.popen("gzip -dc " + ARGV[0])
fh2 = IO.popen("gzip -dc " + ARGV[1])
ff1 = Bio::FlatFile.new(nil, fh1)
ff2 = Bio::FlatFile.new(nil, fh2)

ff1.each_entry do |fe1|
  fe2 = ff2.next_entry
  puts "@#{fe1.definition}\n#{fe1.sequence_string}\n+#{fe1.definition}\n#{fe1.quality_string}\n"
  puts "@#{fe2.definition}\n#{fe2.sequence_string}\n+#{fe2.definition}\n#{fe2.quality_string}\n"
end

ruby ~/rubysrc/fastqgzmerge.rb /path/to/files/s_8_?_4_sequence.txt.gz > s_8_4.sequence.txt

のように発行する事ができる。所要時間は非圧縮とほとんど変わらない。展開はgzipが並行にやってくれるし、rubyの処理部分がgzipよりずっと遅いということ。

$ echo  'ruby ~/rubysrc/fastqgzmerge.rb /path/to/files/s_8_?_5_sequence.txt.gz > s_8_5.sequence.txt' | qsub -cwd

などのようにすればよい。リダイレクト">"させるので、echo のあとのコマンド全体を引用符で囲う。