code_ruby/lib/maasha/sam.rb

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