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, :header
  42
  43   # Method to initialize a Sam object.
  44   def initialize(io = nil)
  45     @io     = io
  46     @header = {}
  47
  48     parse_header
  49   end
  50
  51   def each
  52     @io.each_line do |line|
  53       unless line[0] == '@'
  54         entry = parse_alignment(line.chomp)
  55
  56         yield entry if block_given?
  57       end
  58     end
  59   end
  60
  61   private
  62
  63   # Method to parse the header section of a SAM file.
  64   # Each header line should match:
  65   # /^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ or /^@CO\t.*/.
  66   # Tags containing lowercase letters are reserved for end users.
  67   def parse_header
  68     @io.each_line do |line|
  69       if line =~ /^@([A-Za-z][A-Za-z])/
  70         line.chomp!
  71
  72         tag = $1
  73
  74         case tag
  75         when 'HD' then subparse_header(line)
  76         when 'SQ' then subparse_sequence(line)
  77         when 'RG' then subparse_read_group(line)
  78         when 'PG' then subparse_program(line)
  79         when 'CO' then subparse_comment(line)
  80         else
  81           raise SamError, "Unknown header tag: #{tag}"
  82         end
  83       else
  84         @io.rewind
  85         break
  86       end
  87     end
  88
  89     return @header.empty? ? nil : @header
  90   end
  91
  92   # Method to subparse header lines.
  93   def subparse_header(line)
  94     hash   = {}
  95     fields = line.split("\t")
  96
  97     if fields[1] =~ /^VN:([0-9]+\.[0-9]+)$/
  98       hash[:VN] = $1.to_f
  99     else
 100       raise SamError, "Bad version number: #{fields[1]}"
 101     end
 102
 103     if fields.size > 2
 104       if fields[2] =~ /^SO:(unknown|unsorted|queryname|coordinate)$/
 105         hash[:SO] = $1
 106       else
 107         raise SamError, "Bad sort order: #{fields[2]}"
 108       end
 109     end
 110
 111     @header[:HD] = hash
 112   end
 113
 114   # Method to subparse sequence lines.
 115   def subparse_sequence(line)
 116     @header[:SQ] = Hash.new unless @header[:SQ].is_a? Hash
 117     hash = {}
 118
 119     fields = line.split("\t")
 120
 121     if fields[1] =~ /^SN:([!-)+-<>-~][!-~]*)$/
 122       seq_name = $1.to_sym
 123     else
 124       raise SamError, "Bad sequence name: #{fields[1]}"
 125     end
 126
 127     if fields[2] =~ /^LN:(\d+)$/
 128       hash[:LN] = $1.to_i
 129     else
 130       raise SamError, "Bad sequence length: #{fields[2]}"
 131     end
 132
 133     (3 ... fields.size).each do |i|
 134       if fields[i] =~ /^(AS|M5|SP|UR):([ -~]+)$/
 135         hash[$1.to_sym] = $2
 136       else
 137         raise SamError, "Bad sequence tag: #{fields[i]}"
 138       end
 139     end
 140
 141     @header[:SQ][:SN] = Hash.new unless @header[:SQ][:SN].is_a? Hash
 142
 143     if @header[:SQ][:SN].has_key? seq_name
 144       raise SamError, "Non-unique sequence name: #{seq_name}"
 145     else
 146       @header[:SQ][:SN][seq_name] = hash
 147     end
 148   end
 149
 150   # Method to subparse read group lines.
 151   def subparse_read_group(line)
 152     @header[:RG] = Hash.new unless @header[:RG].is_a? Hash
 153     hash = {}
 154
 155     fields = line.split("\t")
 156
 157     if fields[1] =~ /^ID:([ -~]+)$/
 158       id = $1.to_sym
 159     else
 160       raise SamError, "Bad read group identifier: #{fields[1]}"
 161     end
 162
 163     (2 ... fields.size).each do |i|
 164       if fields[i] =~ /^(CN|DS|DT|FO|KS|LB|PG|PI|PL|PU|SM):([ -~]+)$/
 165         hash[$1.to_sym] = $2
 166       else
 167         raise SamError, "Bad read group tag: #{fields[i]}"
 168       end
 169     end
 170
 171     if hash.has_key? :FO
 172       unless hash[:FO] =~ /^\*|[ACMGRSVTWYHKDBN]+$/
 173         raise SamError, "Bad flow order: #{hash[:FO]}"
 174       end
 175     end
 176
 177     if hash.has_key? :PL
 178       unless hash[:PL] =~ /^(CAPILLARY|LS454|ILLUMINA|SOLID|HELICOS|IONTORRENT|PACBIO)$/
 179         raise SamError, "Bad platform: #{hash[:PL]}"
 180       end
 181     end
 182
 183     @header[:RG][:ID] = Hash.new unless @header[:RG][:ID].is_a? Hash
 184
 185     if @header[:RG][:ID].has_key? id
 186       raise SamError, "Non-unique read group identifier: #{id}"
 187     else
 188       @header[:RG][:ID][id] = hash
 189     end
 190   end
 191
 192   # Method to subparse program lines.
 193   def subparse_program(line)
 194     @header[:PG] = Hash.new unless @header[:PG].is_a? Hash
 195     hash = {}
 196
 197     fields = line.split("\t")
 198
 199     if fields[1] =~ /^ID:([ -~]+)$/
 200       id = $1.to_sym
 201     else
 202       raise SamError, "Bad program record identifier: #{fields[1]}"
 203     end
 204
 205     (2 ... fields.size).each do |i|
 206       if fields[i] =~ /^(PN|CL|PP|VN):([ -~]+)$/
 207         hash[$1.to_sym] = $2
 208       else
 209         raise SamError, "Bad program record tag: #{fields[i]}"
 210       end
 211     end
 212
 213     @header[:PG][:ID] = Hash.new unless @header[:PG][:ID].is_a? Hash
 214
 215     if @header[:PG][:ID].has_key? id
 216       raise SamError, "Non-unique program record identifier: #{id}"
 217     else
 218       @header[:PG][:ID][id] = hash
 219     end
 220   end
 221
 222   # Method to subparse comment lines.
 223   def subparse_comment(line)
 224     @header[:CO] = Array.new unless @header[:CO].is_a? Array
 225
 226     if line =~ /^@CO\t(.+)/
 227       @header[:CO] << $1
 228     else
 229       raise SamError, "Bad comment line: #{line}"
 230     end
 231   end
 232
 233   # Method to subparse alignment lines.
 234   def parse_alignment(line)
 235     fields = line.split("\t")
 236
 237     raise SamError, "Bad number of fields: #{fields.size} < 11" if fields.size < 11
 238
 239     qname = fields[0]
 240     flag  = fields[1].to_i
 241     rname = fields[2]
 242     pos   = fields[3].to_i
 243     mapq  = fields[4].to_i
 244     cigar = fields[5]
 245     rnext = fields[6]
 246     pnext = fields[7].to_i
 247     tlen  = fields[8].to_i
 248     seq   = fields[9]
 249     qual  = fields[10]
 250
 251     check_qname(qname)
 252     check_flag(flag)
 253     check_rname(rname)
 254     check_pos(pos)
 255     check_mapq(mapq)
 256     check_cigar(cigar, seq)
 257     check_rnext(rnext)
 258     check_pnext(pnext)
 259     check_tlen(tlen)
 260     check_seq(seq)
 261     check_qual(qual)
 262
 263     entry = {}
 264     entry[:QNAME] = qname
 265     entry[:FLAG]  = flag
 266     entry[:RNAME] = rname
 267     entry[:POS]   = pos
 268     entry[:MAPQ]  = mapq
 269     entry[:CIGAR] = cigar
 270     entry[:RNEXT] = rnext
 271     entry[:PNEXT] = pnext
 272     entry[:TLEN]  = tlen
 273     entry[:SEQ]   = Seq.new(qname, seq)
 274     entry[:QUAL]  = qual
 275
 276     entry
 277   end
 278
 279   # Method to check qname.
 280   def check_qname(qname)
 281     raise SamError, "Bad qname: #{qname}" unless qname =~ /^[!-?A-~]{1,255}$/
 282   end
 283
 284   # Method to check flag.
 285   def check_flag(flag)
 286     raise SamError, "Bad flag: #{flag}" unless (0 .. 2**16 - 1).include? flag
 287   end
 288
 289   # Method to check if rname, when not '*' and
 290   # @SQ header lines are present, is located in
 291   # the header hash.
 292   def check_rname(rname)
 293     raise SamError, "Bad rname: #{rname}" unless rname =~ /^(\*|[!-()+-<>-~][!-~]*)$/
 294
 295     unless @header.empty? or rname == '*'
 296       unless @header[:SQ][:SN].has_key? rname.to_sym
 297         raise SamError, "rname not found in header hash: #{rname}"
 298       end
 299     end
 300   end
 301
 302   # Method to check pos.
 303   def check_pos(pos)
 304     raise SamError, "Bad pos: #{pos}" unless (0 .. 2**29 - 1).include? pos
 305   end
 306
 307   # Method to check mapq.
 308   def check_mapq(mapq)
 309     raise SamError, "Bad mapq: #{mapq}" unless (0 .. 2**8 - 1).include? mapq
 310   end
 311
 312   # Method to check cigar.
 313   def check_cigar(cigar, seq)
 314     raise SamError, "Bad cigar: #{cigar}" unless cigar =~ /^(\*|([0-9]+[MIDNSHPX=])+)$/
 315
 316     # Check cigar hard clipping only at ends.
 317     if cigar.gsub(/^[0-9]+H|[0-9]+H$/, "").match('H')
 318       raise SamError, "Bad cigar with internal H: #{cigar}"
 319     end
 320
 321     # Check cigar soft clipping only at ends or H.
 322     if cigar.gsub(/^[0-9]+H|[0-9]+H$/, "").gsub(/^[0-9]+S|[0-9]+S$/, "").match('S')
 323       raise SamError, "Bad cigar with internal S: #{cigar}"
 324     end
 325
 326     # Check cigar length matches sequence length.
 327     unless cigar == '*' or seq == '*'
 328       cigar_len = 0
 329
 330       cigar.scan(/([0-9]+)([MIDNSHPX=])/).each do |len, op|
 331         cigar_len += len.to_i if op =~ /[MISX=]/
 332       end
 333
 334       if cigar_len != seq.length
 335         raise SamError, "cigar and sequence length mismatch: #{cigar_len} != #{seq.length}"
 336       end
 337     end
 338   end
 339
 340   # Method to check if rnext, when not '*' or '='
 341   # and @SQ header lines are present, is located
 342   # in the header hash.
 343   def check_rnext(rnext)
 344     raise SamError, "Bad rnext: #{rnext}" unless rnext =~ /^(\*|=|[!-()+-<>-~][!-~]*)$/
 345
 346     unless @header.empty? or rnext == '*' or rnext == '='
 347       unless @header[:SQ][:SN].has_key? rnext.to_sym
 348         raise SamError, "rnext not found in header hash: #{rnext}"
 349       end
 350     end
 351   end
 352
 353   # Method to check pnext.
 354   def check_pnext(pnext)
 355     raise SamError, "Bad pnext: #{pnext}" unless (0 .. 2**29 - 1).include? pnext
 356   end
 357
 358   # Method to check tlen.
 359   def check_tlen(tlen)
 360     raise SamError, "Bad tlen: #{tlen}" unless (-2**29 + 1 .. 2**29 - 1).include? tlen
 361   end
 362
 363   # Method to check seq.
 364   def check_seq(seq)
 365     raise SamError, "Bad seq: #{seq}" unless seq  =~ /^(\*|[A-Za-z=.]+)$/
 366   end
 367
 368   # Method to check qual.
 369   def check_qual(qual)
 370     raise SamError, "Bad qual: #{qual}" unless qual =~ /^[!-~]+$/
 371   end
 372 end
 373
 374
 375 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 376
 377
 378 __END__
 379