Changeset 3395

Timestamp:

Mar 23, 2015, 6:31:15 AM (8 years ago)

Author:

Peter

Message:

closes #830. New function that calculates alignment_score of a bam alignment. Also a new function that returns bams one char in bam's query sequence. Using this function in sequence function. This might slow down the code slightly when building against old libbam, but speed is not a priority in that case.

Location:

trunk

Files:

• test/bam.cc (modified) (2 diffs)
• yat/omic/BamRead.cc (modified) (1 diff)
• yat/omic/BamRead.h (modified) (2 diffs)

trunk/test/bam.cc

@@ -45 +45 @@
 #ifdef YAT_HAVE_LIBBAM
 using namespace omic;
+
+void test_alignment_score(test::Suite& suite, const BamRead& b,
+                          const std::string& ref, double correct)
+{
+  int score = alignment_score(b, ref.begin(), 3, 2, 5);
+  if (!suite.equal(score, correct)) {
+    suite.add(false);
+    suite.err() << "error: alignment_score: " << score
+                << " expected " << correct << "\n";
+    suite.err() << "query: " << b.sequence() << "\n";
+    suite.err() << "cigar: " << b.cigar_str() << "\n";
+    suite.err() << "ref:   " << ref << "\n";
+  }
+}
+
+
+void test_alignment_score(test::Suite& suite, const BamRead& b)
+{
+  std::string reference = b.sequence();
+  test_alignment_score(suite, b, reference, 100);
+  reference[48] = 'N';
+  test_alignment_score(suite, b, reference, 96);
+
+  utility::Aligner::Cigar cig;
+  cig.push_back(BAM_CEQUAL, 48);
+  cig.push_back(BAM_CDIFF, 1);
+  cig.push_back(BAM_CEQUAL, 51);
+  BamRead bam(b);
+  bam.cigar(cig);
+  test_alignment_score(suite, bam, reference, 96);
+
+  cig.clear();
+  cig.push_back(BAM_CSOFT_CLIP, 10);
+  cig.push_back(BAM_CEQUAL, 38);
+  cig.push_back(BAM_CDIFF, 1);
+  cig.push_back(BAM_CEQUAL, 51);
+  bam.cigar(cig);
+  test_alignment_score(suite, bam, reference, 86);
+
+  cig.clear();
+  cig.push_back(BAM_CSOFT_CLIP, 10);
+  cig.push_back(BAM_CEQUAL, 38);
+  cig.push_back(BAM_CDEL, 2);
+  cig.push_back(BAM_CEQUAL, 10);
+  cig.push_back(BAM_CINS, 2);
+  cig.push_back(BAM_CEQUAL, 10);
+  cig.push_back(BAM_CSOFT_CLIP, 2);
+  bam.cigar(cig);
+  test_alignment_score(suite, bam, reference, 58 - (5+2*2) - (5+2*2));
+}
+
 
 void test_aux(test::Suite& suite, const BamRead& b)
@@ -310 +361 @@
                 << " expected " << new_q << "\n";
   }
+  test_alignment_score(suite, bam);
   test_cigar(suite, bam, in.header());
   test_aux(suite, bam);

trunk/yat/omic/BamRead.cc

@@ -340 +340 @@
   {
     std::string result(sequence_length(), ' ');
-#if YAT_HAVE_HTSLIB
     for (size_t i=0; i<result.size(); ++i)
-      result[i] = seq_nt16_str[sequence(i)];
+      result[i] = sequence_str(i);
+    return result;
+  }
+
+
+  char BamRead::sequence_str(size_t i) const
+  {
+#if YAT_HAVE_HTSLIB
+    return seq_nt16_str[sequence(i)];
 #elif HAVE_BAM_NT16_REV_TABLE
-    for (size_t i=0; i<result.size(); ++i)
-      result[i] = bam_nt16_rev_table[sequence(i)];
+    return bam_nt16_rev_table[sequence(i)];
 #else
     std::string table = "=ACMGRSVTWYHKDBN";
-    for (size_t i=0; i<result.size(); ++i)
-      result[i] = table[sequence(i)];
-#endif
-    return result;
+    return table[sequence(i)];
+#endif
   }
 

trunk/yat/omic/BamRead.h

@@ -266 +266 @@
      */
     uint8_t sequence(size_t index) const;
+
+    /**
+       \return A, C, G, T, or N.
+
+       \since New in yat 0.13
+     */
+    char sequence_str(size_t index) const;
 
     /**
@@ -439 +446 @@
   };
 
+
+  /**
+     Calculate alignment score based on \a bam's cigar and sequence
+     and the reference as in [\a ref, ref + bam.end()-bam.pos()).
+
+     For indels the score is calculated as open_gap + len * gap.  For
+     mismatches the score is calculated -mismatch. For matches the
+     score is unity.
+
+     For cigar element BAM_CMATCH the sequence and reference is used
+     to differentiate match from mismatch.
+
+     \since New in yat 0.13
+   */
+  template<typename RandomAccessIterator>
+  double alignment_score(const BamRead& bam, RandomAccessIterator ref,
+                         double mismatch, double gap=0.0, double open_gap=0.0)
+  {
+    double res = 0;
+    size_t q = 0;
+
+    for (uint32_t i = 0; i<bam.core().n_cigar; ++i) {
+      switch (bam.cigar_op(i)) {
+      case (BAM_CMATCH): {
+        RandomAccessIterator end = ref + bam.cigar_oplen(i);
+        for (; ref!=end; ++ref) {
+          if (bam.sequence_str(q) == *ref)
+            res += 1.0;
+          else
+            res -= mismatch;
+          ++q;
+        }
+        break;
+      }
+      case (BAM_CEQUAL): {
+        q += bam.cigar_oplen(i);
+        ref += bam.cigar_oplen(i);
+        res += bam.cigar_oplen(i);
+        break;
+      }
+      case (BAM_CDIFF): {
+        q += bam.cigar_oplen(i);
+        ref += bam.cigar_oplen(i);
+        res -= bam.cigar_oplen(i)*mismatch;
+        break;
+      }
+      case (BAM_CSOFT_CLIP): {
+        q += bam.cigar_oplen(i);
+        break;
+      }
+      case (BAM_CINS): {
+        q += bam.cigar_oplen(i);
+        res -= open_gap + bam.cigar_oplen(i) * gap;
+        break;
+      }
+      case (BAM_CDEL): {
+        ref += bam.cigar_oplen(i);
+        res -= open_gap + bam.cigar_oplen(i) * gap;
+        break;
+      }
+      } // end switch
+    }
+    return res;
+  }
+
 }}}
 #endif