X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=code_ruby%2Flib%2Fmaasha%2Fsam.rb;h=8772e6382f20913d2fca3c9076b90475d2824564;hb=a80169a9121e8537f169cd85010d2ceae3a8d4fd;hp=b82c6b1b2742f6f3757d1b90c8a598502ec9844b;hpb=9cbef4726264d5f04d5803a7551510e5304f13f8;p=biopieces.git

diff --git a/code_ruby/lib/maasha/sam.rb b/code_ruby/lib/maasha/sam.rb
index b82c6b1..8772e63 100644
--- a/code_ruby/lib/maasha/sam.rb
+++ b/code_ruby/lib/maasha/sam.rb
@@ -1,4 +1,4 @@
-# Copyright (C) 2007-2011 Martin A. Hansen.
+# Copyright (C) 2007-2013 Martin A. Hansen.
 
 # This program is free software; you can redistribute it and/or
 # modify it under the terms of the GNU General Public License
@@ -29,17 +29,141 @@
 require 'pp'
 require 'maasha/filesys'
 require 'maasha/seq'
-
-# Error class for all exceptions to do with Genbank.
-class SamError < StandardError; end
+require 'maasha/cigar'
 
 REGEX_HEADER  = Regexp.new(/^@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/)
 REGEX_COMMENT = Regexp.new(/^@CO\t.*/)
 
+FLAG_MULTI               = 0x1   # Template having multiple fragments in sequencing
+FLAG_ALIGNED             = 0x2   # Each fragment properly aligned according to the aligner
+FLAG_UNMAPPED            = 0x4   # Fragment unmapped
+FLAG_NEXT_UNMAPPED       = 0x8   # Next fragment in the template unmapped
+FLAG_REVCOMP             = 0x10  # SEQ being reverse complemented
+FLAG_NEXT_REVCOMP        = 0x20  # SEQ of the next fragment in the template being reversed
+FLAG_FIRST               = 0x40  # The first fragment in the template
+FLAG_LAST                = 0x80  # The last fragment in the template
+FLAG_SECONDARY_ALIGNMENT = 0x100 # Secondary alignment
+FLAG_QUALITY_FAIL        = 0x200 # Not passing quality controls
+FLAG_DUPLICATES          = 0x400 # PCR or optical duplicate
+
+# Error class for all exceptions to do with Genbank.
+class SamError < StandardError; end
+
 # Class to parse and write SAM files.
 class Sam < Filesys
   attr_accessor :io, :header
 
+  # Class method to convert a SAM entry
+  # to a Biopiece record.
+  def self.to_bp(sam)
+    bp = {}
+
+    bp[:REC_TYPE] = 'SAM'
+    bp[:Q_ID]     = sam[:QNAME]
+    bp[:STRAND]   = sam[:FLAG].revcomp? ? '-' : '+'
+    bp[:S_ID]     = sam[:RNAME]
+    bp[:S_BEG]    = sam[:POS]
+    bp[:S_END]    = sam[:POS] + sam[:SEQ].length - 1
+    bp[:MAPQ]     = sam[:MAPQ]
+    bp[:CIGAR]    = sam[:CIGAR].to_s
+
+    unless sam[:RNEXT] == '*'
+      bp[:Q_ID2]  = sam[:RNEXT]
+      bp[:S_BEG2] = sam[:PNEXT]
+      bp[:TLEN]   = sam[:TLEN]
+    end
+
+    bp[:SEQ]      = sam[:SEQ].seq
+
+    unless sam[:SEQ].qual.nil?
+      bp[:SCORES] = sam[:SEQ].qual_convert!(:base_33, :base_64).qual
+    end
+
+    if sam[:NM] and sam[:NM].to_i > 0
+      bp[:NM] = sam[:NM]
+      bp[:MD] = sam[:MD]
+      bp[:ALIGN] = self.align_descriptors(sam)
+    end
+
+    bp
+  end
+
+  # Class method to create a new SAM entry from a Biopiece record.
+  # FIXME
+  def self.new_bp(bp)
+    qname = bp[:Q_ID]
+    flag  = 0
+    rname = bp[:S_ID]
+    pos   = bp[:S_BEG]
+    mapq  = bp[:MAPQ]
+    cigar = bp[:CIGAR]
+    rnext = bp[:Q_ID2]  || '*'
+    pnext = bp[:S_BEG2] || 0
+    tlen  = bp[:TLEN]   || 0
+    seq   = bp[:SEQ]
+    qual  = bp[:SCORES] || '*'
+    nm    = "NM:i:#{bp[:NM]}" if bp[:NM]
+    md    = "MD:Z:#{bp[:MD]}" if bp[:MD]
+
+    flag |= FLAG_REVCOMP if bp[:STRAND] == '+'
+
+    if qname && flag && rname && pos && mapq && cigar && rnext && pnext && tlen && seq && qual
+      ary = [qname, flag, rname, pos, mapq, cigar, rnext, pnext, tlen, seq, qual]
+      ary << nm if nm
+      ary << md if md
+      
+      ary.join("\t")
+    end
+  end
+
+  # Create alignment descriptors according to the KISS
+  # format description:
+  # http://code.google.com/p/biopieces/wiki/KissFormat
+  def self.align_descriptors(sam)
+    offset = 0
+    align  = []
+
+    # Insertions
+    sam[:CIGAR].each do |len, op|
+      if op == 'I'
+        (0 ... len).each_with_index do |i|
+          nt = sam[:SEQ].seq[offset + i]
+
+          align << [offset + i, "->#{nt}"]
+        end
+      end
+
+      offset += len
+    end
+
+    offset    = 0
+    deletions = 0
+
+    sam[:MD].scan(/\d+|\^[A-Z]+|[A-Z]+/).each do |m|
+      if m =~ /\d+/      # Matches
+        offset += m.to_i
+      elsif m[0] == '^'  # Deletions
+        m.each_char do |nt|
+          unless nt == '^'
+            align << [offset, "#{nt}>-"]
+            deletions += 1
+            offset    += 1
+          end
+        end
+      else               # Mismatches
+        m.each_char do |nt|
+          nt2 = sam[:SEQ].seq[offset - deletions]
+
+          align << [offset, "#{nt}>#{nt2}"]
+
+          offset += 1
+        end
+      end
+    end
+
+    align.sort_by { |a| a.first }.map { |k, v| "#{k}:#{v}" }.join(",")
+  end
+
   # Method to initialize a Sam object.
   def initialize(io = nil)
     @io     = io
@@ -140,7 +264,7 @@ class Sam < Filesys
 
     @header[:SQ][:SN] = Hash.new unless @header[:SQ][:SN].is_a? Hash
 
-    if @header[:SQ][:SN].has_key? seq_name
+    if @header[:SQ][:SN][seq_name]
       raise SamError, "Non-unique sequence name: #{seq_name}"
     else
       @header[:SQ][:SN][seq_name] = hash
@@ -168,13 +292,13 @@ class Sam < Filesys
       end
     end
 
-    if hash.has_key? :FO
+    if hash[:FO]
       unless hash[:FO] =~ /^\*|[ACMGRSVTWYHKDBN]+$/
         raise SamError, "Bad flow order: #{hash[:FO]}"
       end
     end
 
-    if hash.has_key? :PL
+    if hash[:PL]
       unless hash[:PL] =~ /^(CAPILLARY|LS454|ILLUMINA|SOLID|HELICOS|IONTORRENT|PACBIO)$/
         raise SamError, "Bad platform: #{hash[:PL]}"
       end
@@ -182,7 +306,7 @@ class Sam < Filesys
 
     @header[:RG][:ID] = Hash.new unless @header[:RG][:ID].is_a? Hash
 
-    if @header[:RG][:ID].has_key? id
+    if @header[:RG][:ID][id]
       raise SamError, "Non-unique read group identifier: #{id}"
     else
       @header[:RG][:ID][id] = hash
@@ -212,7 +336,7 @@ class Sam < Filesys
 
     @header[:PG][:ID] = Hash.new unless @header[:PG][:ID].is_a? Hash
 
-    if @header[:PG][:ID].has_key? id
+    if @header[:PG][:ID][id]
       raise SamError, "Non-unique program record identifier: #{id}"
     else
       @header[:PG][:ID][id] = hash
@@ -253,7 +377,6 @@ class Sam < Filesys
     check_rname(rname)
     check_pos(pos)
     check_mapq(mapq)
-    check_cigar(cigar, seq)
     check_rnext(rnext)
     check_pnext(pnext)
     check_tlen(tlen)
@@ -262,17 +385,40 @@ class Sam < Filesys
 
     entry = {}
     entry[:QNAME] = qname
-    entry[:FLAG]  = flag
+    entry[:FLAG]  = Flag.new(flag)
     entry[:RNAME] = rname
     entry[:POS]   = pos
     entry[:MAPQ]  = mapq
-    entry[:CIGAR] = cigar
+    entry[:CIGAR] = Cigar.new(cigar)
     entry[:RNEXT] = rnext
     entry[:PNEXT] = pnext
     entry[:TLEN]  = tlen
-    entry[:SEQ]   = Seq.new(qname, seq)
+    entry[:SEQ]   = (qual == '*') ? Seq.new(seq_name: qname, seq: seq) : Seq.new(seq_name: qname, seq: seq, qual: qual)
     entry[:QUAL]  = qual
 
+    # Optional fields - where some are really important! HATE HATE HATE SAM!!!
+    
+    fields[11 .. -1].each do |field|
+      tag, type, val = field.split(':')
+
+      raise SamError, "Non-unique optional tag: #{tag}" if entry[tag.to_sym]
+
+      # A [!-~] Printable character
+
+      # i [-+]?[0-9]+ Singed 32-bit integer
+      if type == 'i'
+        raise SamError, "Bad tag in optional field: #{field}" unless val =~ /^[-+]?[0-9]+$/
+        val = val.to_i
+      end
+
+      # f [-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)? Single-precision floating number
+      # Z [ !-~]+ Printable string, including space
+      # H [0-9A-F]+ Byte array in the Hex format
+      # B [cCsSiIf](,[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)+ Integer or numeric array
+
+      entry[tag.to_sym] = val
+    end
+
     entry
   end
 
@@ -293,7 +439,7 @@ class Sam < Filesys
     raise SamError, "Bad rname: #{rname}" unless rname =~ /^(\*|[!-()+-<>-~][!-~]*)$/
 
     unless @header.empty? or rname == '*'
-      unless @header[:SQ][:SN].has_key? rname.to_sym
+      unless @header[:SQ][:SN][rname.to_sym]
         raise SamError, "rname not found in header hash: #{rname}"
       end
     end
@@ -309,44 +455,6 @@ class Sam < Filesys
     raise SamError, "Bad mapq: #{mapq}" unless (0 .. 2**8 - 1).include? mapq
   end
 
-  # Method to check cigar string.
-  def check_cigar(cigar, seq)
-    raise SamError, "Bad cigar: #{cigar}" unless cigar =~ /^(\*|([0-9]+[MIDNSHPX=])+)$/
-
-    unless cigar == '*'
-      check_cigar_hard_clip(cigar)
-      check_cigar_soft_clip(cigar)
-      check_cigar_seq_len(cigar, seq) unless seq == '*'
-    end
-  end
-
-  # Method to check cigar hard clipping only at ends.
-  def check_cigar_hard_clip(cigar)
-    if cigar.gsub(/^[0-9]+H|[0-9]+H$/, "").match('H')
-      raise SamError, "Bad cigar with internal H: #{cigar}"
-    end
-  end
-
-  # Method to check cigar soft clipping only at ends or H.
-  def check_cigar_soft_clip(cigar)
-    if cigar.gsub(/^[0-9]+H|[0-9]+H$/, "").gsub(/^[0-9]+S|[0-9]+S$/, "").match('S')
-      raise SamError, "Bad cigar with internal S: #{cigar}"
-    end
-  end
-
-  # Method to check cigar length matches sequence length.
-  def check_cigar_seq_len(cigar, seq)
-    cigar_len = 0
-
-    cigar.scan(/([0-9]+)([MIDNSHPX=])/).each do |len, op|
-      cigar_len += len.to_i if op =~ /[MISX=]/
-    end
-
-    if cigar_len != seq.length
-      raise SamError, "cigar and sequence length mismatch: #{cigar_len} != #{seq.length}"
-    end
-  end
-
   # Method to check if rnext, when not '*' or '='
   # and @SQ header lines are present, is located
   # in the header hash.
@@ -354,7 +462,7 @@ class Sam < Filesys
     raise SamError, "Bad rnext: #{rnext}" unless rnext =~ /^(\*|=|[!-()+-<>-~][!-~]*)$/
 
     unless @header.empty? or rnext == '*' or rnext == '='
-      unless @header[:SQ][:SN].has_key? rnext.to_sym
+      unless @header[:SQ][:SN][rnext.to_sym]
         raise SamError, "rnext not found in header hash: #{rnext}"
       end
     end
@@ -379,6 +487,73 @@ class Sam < Filesys
   def check_qual(qual)
     raise SamError, "Bad qual: #{qual}" unless qual =~ /^[!-~]+$/
   end
+
+  # Method to deconvolute the SAM flag field.
+  class Flag
+    attr_reader :flag
+
+    # Method to initialize a Flag object.
+    def initialize(flag)
+      @flag = flag
+    end
+
+    # Method to test if template have
+    # multiple fragments in sequencing.
+    def multi?
+      (flag & FLAG_MULTI) == 0
+    end
+
+    # Method to test if each fragment
+    # properly aligned according to the aligner.
+    def aligned?
+      (flag & FLAG_ALIGNED) == 0
+    end
+    
+    # Method to test if the fragment was unmapped.
+    def unmapped?
+      (flag & FLAG_UNMAPPED) == 0
+    end
+
+    # Method to test if the next fragment was unmapped.
+    def next_unmapped?
+      (flag & FLAG_NEXT_UNMAPPED) == 0
+    end
+
+    # Method to test if the fragment was reverse complemented.
+    def revcomp?
+      (flag & FLAG_REVCOMP) == 0
+    end
+
+    # Method to test if the next fragment was reverse complemented.
+    def next_revcomp?
+      (flag & FLAG_NEXT_REVCOMP) == 0
+    end
+
+    # Method to test if the fragment was first in the template.
+    def first?
+      (flag & FLAG_FIRST) == 0
+    end
+
+    # Method to test if the fragment was last in the template.
+    def last?
+      (flag & FLAG_LAST) == 0
+    end
+
+    # Method to test for secondary alignment.
+    def secondary_alignment?
+      (flag & FLAG_SECONDARY_ALIGNMENT) == 0
+    end
+
+    # Method to test for quality fail.
+    def quality_fail?
+      (flag & FLAG_QUALITY_FAIL) == 0
+    end
+
+    # Method to test for PCR or optical duplicates.
+    def duplicates?
+      (flag & FLAG_DUPLICATES) == 0
+    end
+  end
 end