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