2010-07-03

BLASTの結果をSAM形式に変換するスクリプト

BLAST(blastn)の結果を次世代シークエンサで使われるSAM形式に変換するスクリプトです．SAM toolsの中にも，blast2sam.plというperl スクリプトがあり，これのruby移植です．
少し拡張してあって，SAM toolsのblast2sam.plだと配列やquarity valueを出してくれないせいで，tabletなどのSAM/BAMを見られるviewerで読めないことがあります．以下のスクリプトだと，quarity scoreも（適当に）出力するので，tabletで読み込むことが可能です．
使い方は，blastn2sam.rb とすると

% blastn2sam.rb -t 1e-5 input_blast_result > output_sam_file

です．-t の後に示した値以下のE-valueを持つ結果だけsamに変換します．指定しないと全部変換します．
BLASTだとblastnしか対応していませんが，BLATのtblastx相当のオプションでBLAST出力した場合(-t=dnax -q=dnax -out=blast)もparseできます．

#!/usr/bin/env ruby

require 'optparse'

maxeval = 1.0
show_seq = true
opt = OptionParser.new
opt.on("-t", "--threshold [VALUE]", Float){|v| maxeval = v}
opt.parse!(ARGV)

def blast2sam(filename, maxeval, show_seq)
  fp = open(filename, "r")
  query_count = 1
  query_name = ""
  sam = {
    :qname=>'', :flag=>0x00, :rname=>'', :pos=>nil, 
    :mapq=>255, :cigar=>'', :nrnm=>'*', :mpos=>0, 
    :isize=>0, :seq=>'', :qual=>'', :as=>'', :ev=>''}
  cmaux = {:cur_op=>nil, :count=>0} # [0,0,0,'']
  cigar = Array.new
  qlen = 0
  qend = 0
  rname = ""
  query_next = false
  skip_size = 0  # number of spaces before alignment figure
  while line = fp.gets
    if cigar.size !=0 && (line =~ /^Query=/ || line =~ /Score =.*bits.*Expect/) 
      blast_print_sam(sam, cigar, cmaux, qlen - qend) if sam[:ev].to_f <= maxeval
      cigar = Array.new
    end
    if line =~ /^Query= (\S+)/
      query_name = $1
      query_count = 1
    elsif line =~ /\((\S+)\s+letters\)/
      qlen_str = $1
      qlen_str.gsub!(/,/, "")
      qlen = qlen_str.to_i
    elsif line =~ /^>(\S+)/
      rname = $1
    elsif line =~ /Length = (\d+)/
      slen = $1
    elsif line =~ /Score =\s+(\S+) bits.+Expect(\(\d+\))? = (\S+)/
      as = ($1.to_f + 0.499).to_i
      ev = $3
      ev = "1#{ev}" if ev =~ /^e/
      sam[:flag] = 0x49
      sam[:isize] = 0
      sam[:as] = as
      sam[:ev] = ev
      sam[:seq] = ""
      sam[:qual] = ""
      sam[:pos] = 0
      qbeg = 0
      cigar = Array.new
      cmaux = {:cur_op=>nil, :count=>0} # [0,0,0,'']  # @cmaux = (0, 0, 0, '');
      sam[:qname] = "#{query_name}_n#{query_count}" 
      sam[:rname] = rname
      query_count += 1
    elsif line =~ /Strand = (\S+) \/ (\S+)/
      sam[:flag] |= 0x10 if $2 == "Minus"
    elsif line =~ /Query\:\s(\d+)(\s*)(\S+)\s(\d+)/
      start = $1.to_i
      skip_size = 7 + $1.size + $2.size
      q = $3
      if qbeg == 0
        qbeg = start
        cigar << "#{(start-1)}H" if start > 1
      end
      qend = $4.to_i
      if show_seq == true
        sam[:seq] += q.gsub(/-/, "")
      end
      query_next = true
    elsif query_next == true
      query_next = false
      line.chomp!
      align = ""
      0.upto(q.size) do |i|
        align += line[skip_size+i,1] unless q[i,1] == "-"
      end
      align.gsub!(" ", ";")
      sam[:qual] += align
    elsif line =~ /Sbjct\:\s(\d+)\s*(\S+)\s(\d+)/
      s = $2
      if sam[:flag] & 0x10 != 0
        sam[:pos] = $3.to_i
      else
        sam[:pos] = $1.to_i if sam[:pos] == 0
      end
      cigar = aln2cm(cigar, q, s, cmaux)
    end
  end    
  blast_print_sam(sam, cigar, cmaux, qlen - qend) if sam[:ev].to_f <= maxeval
end

def blast_print_sam(sam, cigar, cmaux, qrest)
  cigar << cmaux[:count].to_s + cmaux[:cur_op].to_s
  if sam[:flag] & 0x10 != 0
    cigar = cigar.reverse
    sam[:seq] = sam[:seq].reverse unless sam[:seq].nil?
    sam[:seq].tr!("atgcrymkswATGCRYMKSW","tacgyrkmswTACGYRKMSW")
    sam[:qual] = sam[:qual].reverse unless sam[:qual] == "*"
  end
  sam[:seq] = "*" if sam[:seq].nil?
  sam[:cigar] = cigar.join("")

  puts [sam[:qname],sam[:flag],sam[:rname],sam[:pos], 
        sam[:mapq],sam[:cigar],sam[:nrnm],sam[:mpos], 
        sam[:isize], sam[:seq], sam[:qual], 
        "AS:i:#{sam[:as]}", "EV:Z:#{sam[:ev]}"].join("\t")
end

def aln2cm(cigar, q, s, cmaux)
  l = q.size
  op = nil
  0.upto(l-1) do |i|
    if q[i,1] == "-"
      op = :D
    elsif s[i,1] == "-"
      op = :I
    else
      op = :M
    end
    if cmaux[:cur_op] == op
      cmaux[:count] += 1
    else
      cigar << cmaux[:count].to_s + cmaux[:cur_op].to_s unless cmaux[:cur_op].nil?
      cmaux[:cur_op] = op
      cmaux[:count] = 1
    end
  end
  cigar
end

blast2sam(ARGV[0], maxeval, show_seq)

2010-05-21

次世代シークエンサ(NGS)解析で使われるソフトの簡単なまとめ

research bioinformatics

2010/5/25．EST assemblerのViewerに関して追記
2010/6/9. ちょいちょい追記しています．
2010/7/26. ちょいちょい編集しました．

少し調べたのでまとめを晒してみる．比較的よく使われて沿うなソフトをまとめてみました．既に解析をガシガシやられている先人の方から，こんなソフト使ってるよーとか，そんなん使わん，とか突っ込み歓迎です．あとソフトが見つからないところが有るので，補完コメントいただけると助かります．

主に2つの用途があると思います．

ゲノム配列既知の種にNGSで読んだshort readをmapping
ゲノム配列未知の種からNGSで読んだshort readのESTをアセンブル

(注)NGSの機材のメーカーIllumina/SOLiD/454によっても利用できるソフトウエアが変わります*1が，以下では一緒くたにして書いてあるので注意してください．454は配列が長いので従来のSanger法でのソフトが使えるケースがあります．SOLiDはcolor-space(http://seqanswers.com/forums/showthread.php?t=10 など参照)に対応しているか否かで使えるソフトが分かれます．（みゃー＠ゲノムさん，ご指摘有り難うございます）

順に並べます(ソフトウエアはメーカ謹製のものではなく，アカデミアで自由に使えるものを中心に選んであります）．

ゲノム配列既知の種にNGSで読んだshort readをmapping
- 目標は各short readが，ゲノムのどの位置から出てきたか．手順は，ゲノムへのマップ，形式変換，可視化．
- mappingソフト(BLASTとかBLATだと遅かったりセンシティビティが問題だったり，出力ファイルが軽く数ギガ超えたりするので専用が用いられる)
  - Bowtie(http://bowtie-bio.sourceforge.net/index.shtml)．とりあえず，これ使っとけ，なソフト．BWAに比べると遅いがHadoop対応版であるCrossBowが存在する(http://bowtie-bio.sourceforge.net/crossbow/index.shtml)事や，RNA-seq解析のソフトも用意されている事から使いやすい．
  - MAQ(http://maq.sourceforge.net/)を使った説明が多いよう．でも１世代前のアルゴリズムになっているのと開発が止まっているので，同じ作者による下記のBWAの方を使う
  - BWA(http://bio-bwa.sourceforge.net/bwa.shtml). 汎用，高速．並列化に対応していないのが難点．とりあえず使うには十分．
  - 上記MAQ&BWAの作者による様々なソフトの比較＆まとめ(http://lh3lh3.users.sourceforge.net/NGSalign.shtml)．利点と欠点など良くまとめられている．あんまり並列化に言及が無いのはなぜだろう．
- 形式変換
  - 各mappingソフトから出てきた出力の形式変換
  - GFFとかだと激しく重くて使い勝手悪いので，SAM/BAMもしくは独自形式が利用される．Bowtieは，何も指定しないと独自形式を吐く（オプションでSAMを吐ける）．
  - SAMtools(http://samtools.sourceforge.net/)
- 可視化
  - Gbrowseの次世代シークエンサデータ用チュートリアル．http://gmod.org/wiki/GBrowse_NGS_Tutorial　（takeshi kawashimaさん，ありがとうございます）．
  - ショウジョウバエに対しRNA-Seqの結果をGbrowseを使って可視化した例があります．http://flybase.org/cgi-bin/gbrowse/dmelrnaseq/
  - GbrowseやUSCS genome browserとか今まで使われてたブラウザも使われているっぽいようですが，細かいレベルの塩基の置換とか見るのには適してないような気もします(私見)．ブラウザがインタラクティブに扱える情報量の限界超えているかなぁと．
  - IGV(http://www.broadinstitute.org/igv/home) が使われそうな予感？
  - MapView(http://evolution.sysu.edu.cn/mapview/)：シンプル．
EST assembly (ゲノム配列未知の種からNGSで読んだshort readのESTをアセンブル)
- Assembly
  - BLASTClust. Blastの中に入っているソフト．最短距離法で近い配列を集める．各クラスタに属しているreadは分かるけどalignmentはしてくれない．
  - MIRA3(http://sourceforge.net/apps/mediawiki/mira-assembler/index.php?title=Main_Page). ACE形式にしてViewerでalignmentを見られるっぽい．
  - CAP3...古すぎるよなぁ...
  - ABySS(http://www.bcgsc.ca/platform/bioinfo/software/abyss). MPI使用の並列化も可能。
  - Velvet(http://www.ebi.ac.uk/~zerbino/velvet/).

short readのassemblyが可能．マニュアルでは454対応と書いてあるが，ちょっときつそう．

- 可視化
  - Tablet(http://bioinf.scri.ac.uk/tablet/). ACEファイルだけでなくSOAPdenovoの入力も可能っぽい．
  - GAP5. 使えるかもしれない．要チェック．

様々なソフトへのリンク
- その1: http://seqanswers.com/wiki/Software/list
- その2: http://ngslib.genome.tugraz.at/
三菱スペースソフトの提供している情報。このページよりまとまっているなぁ^^;。http://www.mss-bio.net/ngs/index.html
次世代シークエンサそのもの，何が難しいの？っていう方はこちらのNature Methods 2009年11月のCommentary等参照．http://www.nature.com/nmeth/journal/v6/n11s/full/nmeth.f.268.html

Mapping以後の扱い(RNA-seq等)に関しては，上記では記述していません．

*1:454は400塩基位とSangar法に近い長さの配列を出力するのに対し，IlluminaやSOLiDは現状では100塩基以下の短い連続配列のみが読めます

2010-04-03

2010年第1四半期の1言読書感想文

book