code_ruby/lib/maasha/sam.rb

   1 # Copyright (C) 2007-2011 Martin A. Hansen.
   2
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software
  15 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16
  17 # http://www.gnu.org/copyleft/gpl.html
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # This software is part of the Biopieces framework (www.biopieces.org).
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25 # SAM format version v1.4-r962 - April 17, 2011
  26 #
  27 # http://samtools.sourceforge.net/SAM1.pdf
  28
  29 require 'pp'
  30 require 'maasha/filesys'
  31 require 'maasha/seq'
  32
  33 # Error class for all exceptions to do with Genbank.
  34 class SamError < StandardError; end
  35
  36 REGEX_HEADER  = Regexp.new(/^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/)
  37 REGEX_COMMENT = Regexp.new(/^@CO\t.*/)
  38
  39 # Class to parse and write SAM files.
  40 class Sam < Filesys
  41   attr_accessor :io
  42
  43   # Method to initialize a Sam object.
  44   def initialize(io = nil)
  45     @io          = io
  46     @header_hash = {}
  47   end
  48
  49   # Method to parse the header of a SAM file.
  50   # Each header line should match:
  51   # /^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ or /^@CO\t.*/.
  52   # Tags containing lowercase letters are reserved for end users.
  53   def header
  54     @io.each_line do |line|
  55       if line =~ /^@([A-Za-z][A-Za-z])/
  56         line.chomp!
  57
  58         tag = $1
  59
  60         case tag
  61         when 'HD' then parse_header(line)
  62         when 'SQ' then parse_sequence(line)
  63         when 'RG' then parse_read_group(line)
  64         when 'PG' then parse_program(line)
  65         when 'CO' then parse_comment(line)
  66         else
  67           raise SamError, "Unknown header tag: #{tag}"
  68         end
  69       else
  70         @io.rewind
  71         break
  72       end
  73     end
  74
  75     return @header_hash.empty? ? nil : @header_hash
  76   end
  77
  78   def each
  79     @io.each_line do |line|
  80       unless line[0] == '@'
  81         entry = parse_alignment(line.chomp)
  82
  83         yield entry if block_given?
  84       end
  85     end
  86   end
  87
  88   private
  89
  90   # Method to subparse header lines.
  91   def parse_header(line)
  92     hash   = {}
  93     fields = line.split("\t")
  94
  95     if fields[1] =~ /^VN:([0-9]+\.[0-9]+)$/
  96       hash[:VN] = $1.to_f
  97     else
  98       raise SamError, "Bad version number: #{fields[1]}"
  99     end
 100
 101     if fields.size > 2
 102       if fields[2] =~ /^SO:(unknown|unsorted|queryname|coordinate)$/
 103         hash[:SO] = $1
 104       else
 105         raise SamError, "Bad sort order: #{fields[2]}"
 106       end
 107     end
 108
 109     @header_hash[:HD] = hash
 110   end
 111
 112   # Method to subparse sequence lines.
 113   def parse_sequence(line)
 114     @header_hash[:SQ] = Hash.new unless @header_hash[:SQ].is_a? Hash
 115     hash = {}
 116
 117     fields = line.split("\t")
 118
 119     if fields[1] =~ /^SN:([!-)+-<>-~][!-~]*)$/
 120       seq_name = $1.to_sym
 121     else
 122       raise SamError, "Bad sequence name: #{fields[1]}"
 123     end
 124
 125     if fields[2] =~ /^LN:(\d+)$/
 126       hash[:LN] = $1.to_i
 127     else
 128       raise SamError, "Bad sequence length: #{fields[2]}"
 129     end
 130
 131     (3 ... fields.size).each do |i|
 132       if fields[i] =~ /^(AS|M5|SP|UR):([ -~]+)$/
 133         hash[$1.to_sym] = $2
 134       else
 135         raise SamError, "Bad sequence tag: #{fields[i]}"
 136       end
 137     end
 138
 139     @header_hash[:SQ][:SN] = Hash.new unless @header_hash[:SQ][:SN].is_a? Hash
 140
 141     if @header_hash[:SQ][:SN].has_key? seq_name
 142       raise SamError, "Non-unique sequence name: #{seq_name}"
 143     else
 144       @header_hash[:SQ][:SN][seq_name] = hash
 145     end
 146   end
 147
 148   # Method to subparse read group lines.
 149   def parse_read_group(line)
 150     @header_hash[:RG] = Hash.new unless @header_hash[:RG].is_a? Hash
 151     hash = {}
 152
 153     fields = line.split("\t")
 154
 155     if fields[1] =~ /^ID:([ -~]+)$/
 156       id = $1.to_sym
 157     else
 158       raise SamError, "Bad read group identifier: #{fields[1]}"
 159     end
 160
 161     (2 ... fields.size).each do |i|
 162       if fields[i] =~ /^(CN|DS|DT|FO|KS|LB|PG|PI|PL|PU|SM):([ -~]+)$/
 163         hash[$1.to_sym] = $2
 164       else
 165         raise SamError, "Bad read group tag: #{fields[i]}"
 166       end
 167     end
 168
 169     if hash.has_key? :FO
 170       unless hash[:FO] =~ /^\*|[ACMGRSVTWYHKDBN]+$/
 171         raise SamError, "Bad flow order: #{hash[:FO]}"
 172       end
 173     end
 174
 175     if hash.has_key? :PL
 176       unless hash[:PL] =~ /^(CAPILLARY|LS454|ILLUMINA|SOLID|HELICOS|IONTORRENT|PACBIO)$/
 177         raise SamError, "Bad platform: #{hash[:PL]}"
 178       end
 179     end
 180
 181     @header_hash[:RG][:ID] = Hash.new unless @header_hash[:RG][:ID].is_a? Hash
 182
 183     if @header_hash[:RG][:ID].has_key? id
 184       raise SamError, "Non-unique read group identifier: #{id}"
 185     else
 186       @header_hash[:RG][:ID][id] = hash
 187     end
 188   end
 189
 190   # Method to subparse program lines.
 191   def parse_program(line)
 192     @header_hash[:PG] = Hash.new unless @header_hash[:PG].is_a? Hash
 193     hash = {}
 194
 195     fields = line.split("\t")
 196
 197     if fields[1] =~ /^ID:([ -~]+)$/
 198       id = $1.to_sym
 199     else
 200       raise SamError, "Bad program record identifier: #{fields[1]}"
 201     end
 202
 203     (2 ... fields.size).each do |i|
 204       if fields[i] =~ /^(PN|CL|PP|VN):([ -~]+)$/
 205         hash[$1.to_sym] = $2
 206       else
 207         raise SamError, "Bad program record tag: #{fields[i]}"
 208       end
 209     end
 210
 211     @header_hash[:PG][:ID] = Hash.new unless @header_hash[:PG][:ID].is_a? Hash
 212
 213     if @header_hash[:PG][:ID].has_key? id
 214       raise SamError, "Non-unique program record identifier: #{id}"
 215     else
 216       @header_hash[:PG][:ID][id] = hash
 217     end
 218   end
 219
 220   # Method to subparse comment lines.
 221   def parse_comment(line)
 222     @header_hash[:CO] = Array.new unless @header_hash[:CO].is_a? Array
 223
 224     if line =~ /^@CO\t(.+)/
 225       @header_hash[:CO] << $1
 226     else
 227       raise SamError, "Bad comment line: #{line}"
 228     end
 229   end
 230
 231   # Method to subparse alignment lines.
 232   def parse_alignment(line)
 233     fields = line.split("\t")
 234
 235     raise SamError, "Bad number of fields: #{fields.size} < 11" if fields.size < 11
 236
 237     qname = fields[0]
 238     flag  = fields[1].to_i
 239     rname = fields[2]
 240     pos   = fields[3].to_i
 241     mapq  = fields[4].to_i
 242     cigar = fields[5]
 243     rnext = fields[6]
 244     pnext = fields[7].to_i
 245     tlen  = fields[8].to_i
 246     seq   = fields[9]
 247     qual  = fields[10]
 248
 249     raise SamError, "Bad qname: #{qname}" unless qname =~ /^[!-?A-~]{1,255}$/
 250     raise SamError, "Bad flag: #{flag}"   unless (0 .. 2**16 - 1).include? flag
 251     raise SamError, "Bad rname: #{rname}" unless rname =~ /^(\*|[!-()+-<>-~][!-~]*)$/
 252     raise SamError, "Bad pos: #{pos}"     unless (0 .. 2**29 - 1).include? pos
 253     raise SamError, "Bad mapq: #{mapq}"   unless (0 .. 2**8 - 1).include? mapq
 254     raise SamError, "Bad cigar: #{cigar}" unless cigar =~ /^(\*|([0-9]+[MIDNSHPX=])+)$/
 255     raise SamError, "Bad rnext: #{rnext}" unless rnext =~ /^(\*|=|[!-()+-<>-~][!-~]*)$/
 256     raise SamError, "Bad pnext: #{pnext}" unless (0 .. 2**29 - 1).include? pnext
 257     raise SamError, "Bad tlen: #{tlen}"   unless (-2**29 + 1 .. 2**29 - 1).include? tlen
 258     raise SamError, "Bad seq: #{seq}"     unless seq  =~ /^(\*|[A-Za-z=.]+)$/
 259     raise SamError, "Bad qual: #{qual}"   unless qual =~ /^[!-~]+$/
 260
 261     entry = {}
 262     entry[:QNAME] = qname
 263     entry[:FLAG]  = flag
 264     entry[:RNAME] = rname
 265     entry[:POS]   = pos
 266     entry[:MAPQ]  = mapq
 267     entry[:CIGAR] = cigar
 268     entry[:RNEXT] = rnext
 269     entry[:PNEXT] = pnext
 270     entry[:TLEN]  = tlen
 271     entry[:SEQ]   = seq
 272     entry[:QUAL]  = qual
 273
 274     entry
 275   end
 276 end
 277
 278
 279 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 280
 281
 282 __END__
 283