Changeset 3205

Timestamp:

May 4, 2014, 4:50:50 AM (9 years ago)

Author:

Peter

Message:

implement new class SmithWaterman?

Location:

trunk

Files:

• test/Makefile.am (modified) (1 diff)
• yat/utility/Aligner.cc (modified) (2 diffs)
• yat/utility/Aligner.h (modified) (1 diff)
• yat/utility/SmithWaterman.cc (modified) (2 diffs)
• yat/utility/SmithWaterman.h (modified) (3 diffs)

trunk/test/Makefile.am

@@ -76 +76 @@
   test/roc.test \
   test/score.test test/segment.test test/smart_ptr.test   \
+  test/smith_waterman.test \
   test/smoother.test test/spearman.test \
   test/split.test test/statistics.test test/stream_redirect.test \

trunk/yat/utility/Aligner.cc

@@ -220 +220 @@
   void Aligner::Cigar::push_back(uint8_t op, uint32_t len)
   {
+    if (len==0)
+      return;
     if (cigar_.empty() || this->op(size()-1)!=op)
       cigar_.push_back(bam_cigar_gen(len, op));
@@ -229 +231 @@
   void Aligner::Cigar::push_front(uint8_t op, uint32_t len)
   {
+    if (len==0)
+      return;
     if (cigar_.empty() || this->op(0)!=op)
       cigar_.push_front(bam_cigar_gen(len, op));

trunk/yat/utility/Aligner.h

@@ -120 +120 @@
        gap, and gap + open_gap otherwise.
 
+       To get the CIGAR to behave as expected, the each row in \a dot
+       corresponds to an element in reference sequence and each row in
+       \a dot corresponds to an element in query reference. If that is
+       transposed, insertions and deletions are swapped in CIGAR.
      */
     void operator()(const Matrix& dot, Matrix& score);

trunk/yat/utility/SmithWaterman.cc

@@ -23 +23 @@
 
 #include "SmithWaterman.h"
+#include "yat/omic/BamRead.h" // for cigar defines
 
 #include <cassert>
@@ -30 +31 @@
 namespace utility {
 
+  SmithWaterman::SmithWaterman(double gap, double open_gap)
+    : aligner_(gap, open_gap) {}
+
+
+  SmithWaterman::SmithWaterman(double vertical_gap, double open_vertical_gap,
+                               double horizon_gap, double open_horizon_gap)
+    : aligner_(vertical_gap, open_vertical_gap, horizon_gap, open_horizon_gap)
+  {}
+
+
+  const Aligner::Cigar& SmithWaterman::cigar(void) const
+  {
+    return cigar_;
+  }
+
+
+  size_t SmithWaterman::position(void) const
+  {
+    return position_;
+  }
+
+
+  const Matrix& SmithWaterman::score(void) const
+  {
+    return score_;
+  }
+
+
+  double SmithWaterman::operator()(const Matrix& dot)
+  {
+    score_.resize(dot.rows()+1, dot.columns()+1);
+    aligner_(dot, score_);
+    cigar_.clear();
+    size_t row = 0;
+    size_t column = 0;
+    for (size_t i=0; i<score_.rows(); ++i)
+      for (size_t j=0; j<score_.columns(); ++j)
+        if (score_(i, j) > score_(row, column)) {
+          row = i;
+          column = j;
+        }
+
+    cigar_.push_front(BAM_CSOFT_CLIP, score_.columns()-column-1);
+
+    size_t end_row(row);
+    size_t end_column(column);
+
+    while (aligner_.alignment(row, column)!=Aligner::none) {
+      switch (aligner_.alignment(row, column)) {
+      case(Aligner::diagonal):
+        assert(row);
+        --row;
+        assert(column);
+        --column;
+        if (dot(row, column)>0)
+          cigar_.push_front(BAM_CEQUAL);
+        else
+          cigar_.push_front(BAM_CDIFF);
+        break;
+      case (Aligner::right):
+        assert(column);
+        --column;
+        cigar_.push_front(BAM_CINS);
+        break;
+      case (Aligner::down):
+        assert(row);
+        --row;
+        cigar_.push_front(BAM_CDEL);
+        break;
+      default:;
+      }
+    }
+
+    position_ = row;
+    cigar_.push_front(BAM_CSOFT_CLIP, column);
+
+    return score_(end_row, end_column);
+  }
+
 }}} // of namespace utility, yat, and theplu

trunk/yat/utility/SmithWaterman.h

@@ -43 +43 @@
        Same as SmithWaterman(gap, open_gap, gap, open_gap)
      */
-    SmithWaterman(double gap=0, double open_gap=0);
+    explicit SmithWaterman(double gap=0, double open_gap=0);
 
     /**
-       \brief constructor
+       \brief Constructor
 
        \param vertical_gap Penalty for extending a vertical gap. A
-       vertical gap means alignment consumes an element in second
-       element and not, in the first element, i.e., an element has
-       been inserted to the second element or equivalently an element
-       has been deleted in first sequence.
+       vertical gap means consuming a reference element i.e. a
+       deletion in query sequence.
 
-       \param open_vertical_gap Penalty for open a vertical gap. Total
-       cost for a vertical gap is open_vertical_gap + N *
+       \param open_vertical_gap Open a deletion. Total
+       cost for a deletion is open_vertical_gap + N *
        vertical_gap.
 
-       \param horizon_gap Penalty for extending a horizontal gap
+       \param horizon_gap Penalty for extending a insertion.
 
-       \param open_horizon_gap Penalty for open a horizontal
-       gap. Total penalty for a horizontal gap is open_horizon_gap + N
-       * horizon_gap.
+       \param open_horizon_gap Penalty for open an insertion. Total
+       penalty for a insertion is open_horizon_gap + N * horizon_gap.
      */
     SmithWaterman(double vertical_gap, double open_vertical_gap,
@@ -68 +65 @@
 
     /**
+       The CIGAR here is slightly different compared to CIGAR in
+       Aligner as it uses BAM_CEQUAL and BAM_CDIFF rather than
+       BAM_CMATCH, and it also uses BAM_CSOFT_CLIP and ends of CIGAR
+       to reflect if whole query was aligned.
+
        \return CIGAR reflecting latest performed alignment
      */
@@ -110 +112 @@
   };
 
+  // template implementaion
+
+  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
+  double SmithWaterman::operator()(RandomAccessIterator1 reference_begin,
+                                   RandomAccessIterator1 reference_end,
+                                   RandomAccessIterator2 query_begin,
+                                   RandomAccessIterator2 query_end,
+                                   double mismatch)
+  {
+    Matrix dot(reference_end-reference_begin, query_end-query_begin, -mismatch);
+    for (size_t i=0; i<dot.rows(); ++i)
+      for (size_t j=0; j<dot.columns(); ++j)
+        if (reference_begin[i] == query_begin[j])
+          dot(i, j) = 1;
+    return this->operator()(dot);
+  }
+
 }}} // of namespace utility, yat, and theplu