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 require 'maasha/cigar'
  33
  34 REGEX_HEADER  = Regexp.new(/^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/)
  35 REGEX_COMMENT = Regexp.new(/^@CO\t.*/)
  36
  37 FLAG_MULTI               = 0x1   # Template having multiple fragments in sequencing
  38 FLAG_ALIGNED             = 0x2   # Each fragment properly aligned according to the aligner
  39 FLAG_UNMAPPED            = 0x4   # Fragment unmapped
  40 FLAG_NEXT_UNMAPPED       = 0x8   # Next fragment in the template unmapped
  41 FLAG_REVCOMP             = 0x10  # SEQ being reverse complemented
  42 FLAG_NEXT_REVCOMP        = 0x20  # SEQ of the next fragment in the template being reversed
  43 FLAG_FIRST               = 0x40  # The first fragment in the template
  44 FLAG_LAST                = 0x80  # The last fragment in the template
  45 FLAG_SECONDARY_ALIGNMENT = 0x100 # Secondary alignment
  46 FLAG_QUALITY_FAIL        = 0x200 # Not passing quality controls
  47 FLAG_DUPLICATES          = 0x400 # PCR or optical duplicate
  48
  49 # Error class for all exceptions to do with Genbank.
  50 class SamError < StandardError; end
  51
  52 # Class to parse and write SAM files.
  53 class Sam < Filesys
  54   attr_accessor :io, :header
  55
  56   # Class method to convert a SAM entry
  57   # to a Biopiece record.
  58   def self.to_bp(sam)
  59     bp = {}
  60
  61     bp[:REC_TYPE] = 'SAM'
  62     bp[:Q_ID]   = sam[:QNAME]
  63     bp[:STRAND] = sam[:FLAG].revcomp? ? '-' : '+'
  64     bp[:S_ID]   = sam[:RNAME]
  65     bp[:S_BEG]  = sam[:POS]
  66     bp[:MAPQ]   = sam[:MAPQ]
  67     bp[:CIGAR]  = sam[:CIGAR]
  68
  69     unless sam[:RNEXT] == '*'
  70       bp[:Q_ID2]  = sam[:RNEXT]
  71       bp[:S_BEG2] = sam[:PNEXT]
  72       bp[:TLEN]   = sam[:TLEN]
  73     end
  74
  75     bp[:SEQ]      = sam[:SEQ].seq
  76
  77     unless sam[:SEQ].qual.nil?
  78       bp[:SCORES] = sam[:SEQ].convert_phred2illumina!.qual
  79     end
  80
  81     if sam.has_key? :NM and sam[:NM].to_i > 0
  82       bp[:ALIGN] = self.align_descriptors(sam)
  83     end
  84
  85     bp
  86   end
  87
  88   # Create align descriptors according to the KISS format description:
  89   # http://code.google.com/p/biopieces/wiki/KissFormat
  90   def self.align_descriptors(sam)
  91     offset     = 0
  92     align      = []
  93
  94     # Insertions
  95     sam[:CIGAR].each do |len, op|
  96       if op == 'I'
  97         (0 ... len).each_with_index do |i|
  98           nt = sam[:SEQ].seq[offset + i]
  99
 100           align << [offset + i, "->#{nt}"]
 101         end
 102       end
 103
 104       offset += len
 105     end
 106
 107     offset    = 0
 108     deletions = 0
 109
 110     sam[:MD].scan(/\d+|\^[A-Z]+|[A-Z]+/).each do |m|
 111       if m =~ /\d+/      # Matches
 112         offset += m.to_i
 113       elsif m[0] == '^'  # Deletions
 114         m[1 .. -1].each_char do |nt|
 115           align << [offset, "#{nt}>-"]
 116
 117           deletions += 1
 118           offset    += 1
 119         end
 120       else               # Mismatches
 121         m.each_char do |nt|
 122           nt2 = sam[:SEQ].seq[offset - deletions]
 123
 124           align << [offset, "#{nt}>#{nt2}"]
 125
 126           offset += 1
 127         end
 128       end
 129     end
 130
 131     align.sort_by { |a| a.first }.map { |k, v| "#{k}:#{v}" }.join(",")
 132   end
 133
 134   # Method to initialize a Sam object.
 135   def initialize(io = nil)
 136     @io     = io
 137     @header = {}
 138
 139     parse_header
 140   end
 141
 142   def each
 143     @io.each_line do |line|
 144       unless line[0] == '@'
 145         entry = parse_alignment(line.chomp)
 146
 147         yield entry if block_given?
 148       end
 149     end
 150   end
 151
 152   private
 153
 154   # Method to parse the header section of a SAM file.
 155   # Each header line should match:
 156   # /^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ or /^@CO\t.*/.
 157   # Tags containing lowercase letters are reserved for end users.
 158   def parse_header
 159     @io.each_line do |line|
 160       if line =~ /^@([A-Za-z][A-Za-z])/
 161         line.chomp!
 162
 163         tag = $1
 164
 165         case tag
 166         when 'HD' then subparse_header(line)
 167         when 'SQ' then subparse_sequence(line)
 168         when 'RG' then subparse_read_group(line)
 169         when 'PG' then subparse_program(line)
 170         when 'CO' then subparse_comment(line)
 171         else
 172           raise SamError, "Unknown header tag: #{tag}"
 173         end
 174       else
 175         @io.rewind
 176         break
 177       end
 178     end
 179
 180     return @header.empty? ? nil : @header
 181   end
 182
 183   # Method to subparse header lines.
 184   def subparse_header(line)
 185     hash   = {}
 186     fields = line.split("\t")
 187
 188     if fields[1] =~ /^VN:([0-9]+\.[0-9]+)$/
 189       hash[:VN] = $1.to_f
 190     else
 191       raise SamError, "Bad version number: #{fields[1]}"
 192     end
 193
 194     if fields.size > 2
 195       if fields[2] =~ /^SO:(unknown|unsorted|queryname|coordinate)$/
 196         hash[:SO] = $1
 197       else
 198         raise SamError, "Bad sort order: #{fields[2]}"
 199       end
 200     end
 201
 202     @header[:HD] = hash
 203   end
 204
 205   # Method to subparse sequence lines.
 206   def subparse_sequence(line)
 207     @header[:SQ] = Hash.new unless @header[:SQ].is_a? Hash
 208     hash = {}
 209
 210     fields = line.split("\t")
 211
 212     if fields[1] =~ /^SN:([!-)+-<>-~][!-~]*)$/
 213       seq_name = $1.to_sym
 214     else
 215       raise SamError, "Bad sequence name: #{fields[1]}"
 216     end
 217
 218     if fields[2] =~ /^LN:(\d+)$/
 219       hash[:LN] = $1.to_i
 220     else
 221       raise SamError, "Bad sequence length: #{fields[2]}"
 222     end
 223
 224     (3 ... fields.size).each do |i|
 225       if fields[i] =~ /^(AS|M5|SP|UR):([ -~]+)$/
 226         hash[$1.to_sym] = $2
 227       else
 228         raise SamError, "Bad sequence tag: #{fields[i]}"
 229       end
 230     end
 231
 232     @header[:SQ][:SN] = Hash.new unless @header[:SQ][:SN].is_a? Hash
 233
 234     if @header[:SQ][:SN].has_key? seq_name
 235       raise SamError, "Non-unique sequence name: #{seq_name}"
 236     else
 237       @header[:SQ][:SN][seq_name] = hash
 238     end
 239   end
 240
 241   # Method to subparse read group lines.
 242   def subparse_read_group(line)
 243     @header[:RG] = Hash.new unless @header[:RG].is_a? Hash
 244     hash = {}
 245
 246     fields = line.split("\t")
 247
 248     if fields[1] =~ /^ID:([ -~]+)$/
 249       id = $1.to_sym
 250     else
 251       raise SamError, "Bad read group identifier: #{fields[1]}"
 252     end
 253
 254     (2 ... fields.size).each do |i|
 255       if fields[i] =~ /^(CN|DS|DT|FO|KS|LB|PG|PI|PL|PU|SM):([ -~]+)$/
 256         hash[$1.to_sym] = $2
 257       else
 258         raise SamError, "Bad read group tag: #{fields[i]}"
 259       end
 260     end
 261
 262     if hash.has_key? :FO
 263       unless hash[:FO] =~ /^\*|[ACMGRSVTWYHKDBN]+$/
 264         raise SamError, "Bad flow order: #{hash[:FO]}"
 265       end
 266     end
 267
 268     if hash.has_key? :PL
 269       unless hash[:PL] =~ /^(CAPILLARY|LS454|ILLUMINA|SOLID|HELICOS|IONTORRENT|PACBIO)$/
 270         raise SamError, "Bad platform: #{hash[:PL]}"
 271       end
 272     end
 273
 274     @header[:RG][:ID] = Hash.new unless @header[:RG][:ID].is_a? Hash
 275
 276     if @header[:RG][:ID].has_key? id
 277       raise SamError, "Non-unique read group identifier: #{id}"
 278     else
 279       @header[:RG][:ID][id] = hash
 280     end
 281   end
 282
 283   # Method to subparse program lines.
 284   def subparse_program(line)
 285     @header[:PG] = Hash.new unless @header[:PG].is_a? Hash
 286     hash = {}
 287
 288     fields = line.split("\t")
 289
 290     if fields[1] =~ /^ID:([ -~]+)$/
 291       id = $1.to_sym
 292     else
 293       raise SamError, "Bad program record identifier: #{fields[1]}"
 294     end
 295
 296     (2 ... fields.size).each do |i|
 297       if fields[i] =~ /^(PN|CL|PP|VN):([ -~]+)$/
 298         hash[$1.to_sym] = $2
 299       else
 300         raise SamError, "Bad program record tag: #{fields[i]}"
 301       end
 302     end
 303
 304     @header[:PG][:ID] = Hash.new unless @header[:PG][:ID].is_a? Hash
 305
 306     if @header[:PG][:ID].has_key? id
 307       raise SamError, "Non-unique program record identifier: #{id}"
 308     else
 309       @header[:PG][:ID][id] = hash
 310     end
 311   end
 312
 313   # Method to subparse comment lines.
 314   def subparse_comment(line)
 315     @header[:CO] = Array.new unless @header[:CO].is_a? Array
 316
 317     if line =~ /^@CO\t(.+)/
 318       @header[:CO] << $1
 319     else
 320       raise SamError, "Bad comment line: #{line}"
 321     end
 322   end
 323
 324   # Method to subparse alignment lines.
 325   def parse_alignment(line)
 326     fields = line.split("\t")
 327
 328     raise SamError, "Bad number of fields: #{fields.size} < 11" if fields.size < 11
 329
 330     qname = fields[0]
 331     flag  = fields[1].to_i
 332     rname = fields[2]
 333     pos   = fields[3].to_i
 334     mapq  = fields[4].to_i
 335     cigar = fields[5]
 336     rnext = fields[6]
 337     pnext = fields[7].to_i
 338     tlen  = fields[8].to_i
 339     seq   = fields[9]
 340     qual  = fields[10]
 341
 342     check_qname(qname)
 343     check_flag(flag)
 344     check_rname(rname)
 345     check_pos(pos)
 346     check_mapq(mapq)
 347     check_rnext(rnext)
 348     check_pnext(pnext)
 349     check_tlen(tlen)
 350     check_seq(seq)
 351     check_qual(qual)
 352
 353     entry = {}
 354     entry[:QNAME] = qname
 355     entry[:FLAG]  = Flag.new(flag)
 356     entry[:RNAME] = rname
 357     entry[:POS]   = pos
 358     entry[:MAPQ]  = mapq
 359     entry[:CIGAR] = Cigar.new(cigar)
 360     entry[:RNEXT] = rnext
 361     entry[:PNEXT] = pnext
 362     entry[:TLEN]  = tlen
 363     entry[:SEQ]   = (qual == '*') ? Seq.new(qname, seq) : Seq.new(qname, seq, qual)
 364     entry[:QUAL]  = qual
 365
 366     # Optional fields - where some are really important! HATE HATE HATE SAM!!!
 367
 368     fields[11 .. -1].each do |field|
 369       tag, type, val = field.split(':')
 370
 371       raise SamError, "Non-unique optional tag: #{tag}" if entry.has_key? tag.to_sym
 372
 373       # A [!-~] Printable character
 374
 375       # i [-+]?[0-9]+ Singed 32-bit integer
 376       if type == 'i'
 377         raise SamError, "Bad tag in optional field: #{field}" unless val =~ /^[-+]?[0-9]+$/
 378         val = val.to_i
 379       end
 380
 381       # f [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)? Single-precision floating number
 382       # Z [ !-~]+ Printable string, including space
 383       # H [0-9A-F]+ Byte array in the Hex format
 384       # B [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+ Integer or numeric array
 385
 386       entry[tag.to_sym] = val
 387     end
 388
 389     entry
 390   end
 391
 392   # Method to check qname.
 393   def check_qname(qname)
 394     raise SamError, "Bad qname: #{qname}" unless qname =~ /^[!-?A-~]{1,255}$/
 395   end
 396
 397   # Method to check flag.
 398   def check_flag(flag)
 399     raise SamError, "Bad flag: #{flag}" unless (0 .. 2**16 - 1).include? flag
 400   end
 401
 402   # Method to check if rname, when not '*' and
 403   # @SQ header lines are present, is located in
 404   # the header hash.
 405   def check_rname(rname)
 406     raise SamError, "Bad rname: #{rname}" unless rname =~ /^(\*|[!-()+-<>-~][!-~]*)$/
 407
 408     unless @header.empty? or rname == '*'
 409       unless @header[:SQ][:SN].has_key? rname.to_sym
 410         raise SamError, "rname not found in header hash: #{rname}"
 411       end
 412     end
 413   end
 414
 415   # Method to check pos.
 416   def check_pos(pos)
 417     raise SamError, "Bad pos: #{pos}" unless (0 .. 2**29 - 1).include? pos
 418   end
 419
 420   # Method to check mapq.
 421   def check_mapq(mapq)
 422     raise SamError, "Bad mapq: #{mapq}" unless (0 .. 2**8 - 1).include? mapq
 423   end
 424
 425   # Method to check if rnext, when not '*' or '='
 426   # and @SQ header lines are present, is located
 427   # in the header hash.
 428   def check_rnext(rnext)
 429     raise SamError, "Bad rnext: #{rnext}" unless rnext =~ /^(\*|=|[!-()+-<>-~][!-~]*)$/
 430
 431     unless @header.empty? or rnext == '*' or rnext == '='
 432       unless @header[:SQ][:SN].has_key? rnext.to_sym
 433         raise SamError, "rnext not found in header hash: #{rnext}"
 434       end
 435     end
 436   end
 437
 438   # Method to check pnext.
 439   def check_pnext(pnext)
 440     raise SamError, "Bad pnext: #{pnext}" unless (0 .. 2**29 - 1).include? pnext
 441   end
 442
 443   # Method to check tlen.
 444   def check_tlen(tlen)
 445     raise SamError, "Bad tlen: #{tlen}" unless (-2**29 + 1 .. 2**29 - 1).include? tlen
 446   end
 447
 448   # Method to check seq.
 449   def check_seq(seq)
 450     raise SamError, "Bad seq: #{seq}" unless seq  =~ /^(\*|[A-Za-z=.]+)$/
 451   end
 452
 453   # Method to check qual.
 454   def check_qual(qual)
 455     raise SamError, "Bad qual: #{qual}" unless qual =~ /^[!-~]+$/
 456   end
 457
 458   # Method to deconvolute the SAM flag field.
 459   class Flag
 460     attr_reader :flag
 461
 462     # Method to initialize a Flag object.
 463     def initialize(flag)
 464       @flag = flag
 465     end
 466
 467     # Method to test if template have
 468     # multiple fragments in sequencing.
 469     def multi?
 470       flag & FLAG_MULTI
 471     end
 472
 473     # Method to test if each fragment
 474     # properly aligned according to the aligner.
 475     def aligned?
 476       flag & FLAG_ALIGNED
 477     end
 478
 479     # Method to test if the fragment was unmapped.
 480     def unmapped?
 481       flag & FLAG_UNMAPPED
 482     end
 483
 484     # Method to test if the next fragment was unmapped.
 485     def next_unmapped?
 486       flag & FLAG_NEXT_UNMAPPED
 487     end
 488
 489     # Method to test if the fragment was reverse complemented.
 490     def revcomp?
 491       flag & FLAG_REVCOMP
 492     end
 493
 494     # Method to test if the next fragment was reverse complemented.
 495     def next_revcomp?
 496       flag & FLAG_NEXT_REVCOMP
 497     end
 498
 499     # Method to test if the fragment was first in the template.
 500     def first?
 501       flag & FLAG_FIRST
 502     end
 503
 504     # Method to test if the fragment was last in the template.
 505     def last?
 506       flag & FLAG_LAST
 507     end
 508
 509     # Method to test for secondary alignment.
 510     def secondary_alignment?
 511       flag & FLAG_SECONDARY_ALIGNMENT
 512     end
 513
 514     # Method to test for quality fail.
 515     def quality_fail?
 516       flag & FLAG_QUALITY_FAIL
 517     end
 518
 519     # Method to test for PCR or optical duplicates.
 520     def duplicates?
 521       flag & FLAG_DUPLICATES
 522     end
 523   end
 524 end
 525
 526
 527 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 528
 529
 530 __END__
 531