]> git.donarmstrong.com Git - biopieces.git/blobdiff - code_ruby/lib/maasha/sam.rb
worked on unit tests for sam.rb
[biopieces.git] / code_ruby / lib / maasha / sam.rb
index 900df867d8a70166f91b9a9cfc5d616e08b7f27f..2a3e5868d36aed6ee0c8fa85ef1b249f6dd0842d 100644 (file)
@@ -253,7 +253,7 @@ class Sam < Filesys
     check_rname(rname)
     check_pos(pos)
     check_mapq(mapq)
-    check_cigar(cigar)
+    check_cigar(cigar, seq)
     check_rnext(rnext)
     check_pnext(pnext)
     check_tlen(tlen)
@@ -306,12 +306,25 @@ class Sam < Filesys
 
   # Method to check mapq.
   def check_mapq(mapq)
-    raise SamError, "Bad mapq: #{mapq}"   unless (0 .. 2**8 - 1).include? mapq
+    raise SamError, "Bad mapq: #{mapq}" unless (0 .. 2**8 - 1).include? mapq
   end
 
   # Method to check cigar.
-  def check_cigar(cigar)
+  def check_cigar(cigar, seq)
     raise SamError, "Bad cigar: #{cigar}" unless cigar =~ /^(\*|([0-9]+[MIDNSHPX=])+)$/
+
+    # Check cigar length matches sequence length.
+    unless cigar == '*' or 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
   end
 
   # Method to check if rnext, when not '*' or '='