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         1.upto(m.size - 1).each do
 115           nt = sam[:SEQ].seq[offset]
 116
 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_rnext(rnext)
 350     check_pnext(pnext)
 351     check_tlen(tlen)
 352     check_seq(seq)
 353     check_qual(qual)
 354
 355     entry = {}
 356     entry[:QNAME] = qname
 357     entry[:FLAG]  = Flag.new(flag)
 358     entry[:RNAME] = rname
 359     entry[:POS]   = pos
 360     entry[:MAPQ]  = mapq
 361     entry[:CIGAR] = Cigar.new(cigar)
 362     entry[:RNEXT] = rnext
 363     entry[:PNEXT] = pnext
 364     entry[:TLEN]  = tlen
 365     entry[:SEQ]   = (qual == '*') ? Seq.new(qname, seq) : Seq.new(qname, seq, qual)
 366     entry[:QUAL]  = qual
 367
 368     # Optional fields - where some are really important! HATE HATE HATE SAM!!!
 369
 370     fields[11 .. -1].each do |field|
 371       tag, type, val = field.split(':')
 372
 373       raise SamError, "Non-unique optional tag: #{tag}" if entry.has_key? tag.to_sym
 374
 375       # A [!-~] Printable character
 376
 377       # i [-+]?[0-9]+ Singed 32-bit integer
 378       if type == 'i'
 379         raise SamError, "Bad tag in optional field: #{field}" unless val =~ /^[-+]?[0-9]+$/
 380         val = val.to_i
 381       end
 382
 383       # f [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)? Single-precision floating number
 384       # Z [ !-~]+ Printable string, including space
 385       # H [0-9A-F]+ Byte array in the Hex format
 386       # B [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+ Integer or numeric array
 387
 388       entry[tag.to_sym] = val
 389     end
 390
 391     entry
 392   end
 393
 394   # Method to check qname.
 395   def check_qname(qname)
 396     raise SamError, "Bad qname: #{qname}" unless qname =~ /^[!-?A-~]{1,255}$/
 397   end
 398
 399   # Method to check flag.
 400   def check_flag(flag)
 401     raise SamError, "Bad flag: #{flag}" unless (0 .. 2**16 - 1).include? flag
 402   end
 403
 404   # Method to check if rname, when not '*' and
 405   # @SQ header lines are present, is located in
 406   # the header hash.
 407   def check_rname(rname)
 408     raise SamError, "Bad rname: #{rname}" unless rname =~ /^(\*|[!-()+-<>-~][!-~]*)$/
 409
 410     unless @header.empty? or rname == '*'
 411       unless @header[:SQ][:SN].has_key? rname.to_sym
 412         raise SamError, "rname not found in header hash: #{rname}"
 413       end
 414     end
 415   end
 416
 417   # Method to check pos.
 418   def check_pos(pos)
 419     raise SamError, "Bad pos: #{pos}" unless (0 .. 2**29 - 1).include? pos
 420   end
 421
 422   # Method to check mapq.
 423   def check_mapq(mapq)
 424     raise SamError, "Bad mapq: #{mapq}" unless (0 .. 2**8 - 1).include? mapq
 425   end
 426
 427   # Method to check if rnext, when not '*' or '='
 428   # and @SQ header lines are present, is located
 429   # in the header hash.
 430   def check_rnext(rnext)
 431     raise SamError, "Bad rnext: #{rnext}" unless rnext =~ /^(\*|=|[!-()+-<>-~][!-~]*)$/
 432
 433     unless @header.empty? or rnext == '*' or rnext == '='
 434       unless @header[:SQ][:SN].has_key? rnext.to_sym
 435         raise SamError, "rnext not found in header hash: #{rnext}"
 436       end
 437     end
 438   end
 439
 440   # Method to check pnext.
 441   def check_pnext(pnext)
 442     raise SamError, "Bad pnext: #{pnext}" unless (0 .. 2**29 - 1).include? pnext
 443   end
 444
 445   # Method to check tlen.
 446   def check_tlen(tlen)
 447     raise SamError, "Bad tlen: #{tlen}" unless (-2**29 + 1 .. 2**29 - 1).include? tlen
 448   end
 449
 450   # Method to check seq.
 451   def check_seq(seq)
 452     raise SamError, "Bad seq: #{seq}" unless seq  =~ /^(\*|[A-Za-z=.]+)$/
 453   end
 454
 455   # Method to check qual.
 456   def check_qual(qual)
 457     raise SamError, "Bad qual: #{qual}" unless qual =~ /^[!-~]+$/
 458   end
 459
 460   # Method to deconvolute the SAM flag field.
 461   class Flag
 462     attr_reader :flag
 463
 464     # Method to initialize a Flag object.
 465     def initialize(flag)
 466       @flag = flag
 467     end
 468
 469     # Method to test if template have
 470     # multiple fragments in sequencing.
 471     def multi?
 472       flag & FLAG_MULTI
 473     end
 474
 475     # Method to test if each fragment
 476     # properly aligned according to the aligner.
 477     def aligned?
 478       flag & FLAG_ALIGNED
 479     end
 480
 481     # Method to test if the fragment was unmapped.
 482     def unmapped?
 483       flag & FLAG_UNMAPPED
 484     end
 485
 486     # Method to test if the next fragment was unmapped.
 487     def next_unmapped?
 488       flag & FLAG_NEXT_UNMAPPED
 489     end
 490
 491     # Method to test if the fragment was reverse complemented.
 492     def revcomp?
 493       flag & FLAG_REVCOMP
 494     end
 495
 496     # Method to test if the next fragment was reverse complemented.
 497     def next_revcomp?
 498       flag & FLAG_NEXT_REVCOMP
 499     end
 500
 501     # Method to test if the fragment was first in the template.
 502     def first?
 503       flag & FLAG_FIRST
 504     end
 505
 506     # Method to test if the fragment was last in the template.
 507     def last?
 508       flag & FLAG_LAST
 509     end
 510
 511     # Method to test for secondary alignment.
 512     def secondary_alignment?
 513       flag & FLAG_SECONDARY_ALIGNMENT
 514     end
 515
 516     # Method to test for quality fail.
 517     def quality_fail?
 518       flag & FLAG_QUALITY_FAIL
 519     end
 520
 521     # Method to test for PCR or optical duplicates.
 522     def duplicates?
 523       flag & FLAG_DUPLICATES
 524     end
 525   end
 526 end
 527
 528
 529 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 530
 531
 532 __END__
 533