Changeset 869 for trunk/yat

Timestamp:

Sep 14, 2007, 8:59:46 PM (15 years ago)

Author:

Peter

Message:

Adding Smith-Waterman local alignment and modified ssearch to use this instead. Also added some convenient max functions in stl_utility.

Location:

trunk/yat/utility

Files:

• Alignment.cc (modified) (3 diffs)
• Alignment.h (modified) (2 diffs)
• stl_utility.h (modified) (2 diffs)

trunk/yat/utility/Alignment.cc

@@ -26 +26 @@
 #include "Alignment.h"
 #include "matrix.h"
+#include "yat/utility/stl_utility.h"
 
+#include <algorithm>
 #include <utility>
 #include <vector>
@@ -33 +35 @@
 namespace yat {
 namespace utility {
-
-  unsigned int ssearch(std::string first, std::string second)
-  {
-    matrix m(first.size(), second.size());
-    for (size_t i=0; i<first.size(); ++i)
-      for (size_t j=0; j<second.size(); ++j)
-        m(i,j) = (first[i]==second[j] ? 1 : 0);
-    std::vector<std::pair<size_t, size_t> > path;
-    return static_cast<unsigned int>(NeedlemanWunsch(m,path,0));
-  }
-
 
   double NeedlemanWunsch(const utility::matrix& s,
@@ -95 +86 @@
   }
 
+
+  double ssearch(std::string first, std::string second, double gap, 
+                 double open_gap)
+  {
+    matrix m(first.size(), second.size());
+    for (size_t i=0; i<first.size(); ++i)
+      for (size_t j=0; j<second.size(); ++j)
+        m(i,j) = (first[i]==second[j] ? 1 : 0);
+    return SmithWaterman(m, gap, open_gap);
+  }
+
+
+  double SmithWaterman(const utility::matrix& s,
+                       double gap, double open_gap)
+  {
+    // Calculating S-W matrix
+    matrix m(s.rows()+2,s.columns()+2);
+    for (size_t i=2; i<m.rows(); ++i) 
+      for (size_t j=2; j<m.columns(); ++j) 
+        m(i,j) = max(0.0, m(i-1,j-1)+s(i-2,j-2), 
+                     m(i-1,j)-gap-open_gap, m(i-2,j)-2*gap, 
+                     m(i,j-1)-gap-open_gap, m(i,j-2)-2*gap);
+
+    return max(m);
+  }
+
 }}} // of namespace utility, yat, and theplu

trunk/yat/utility/Alignment.h

@@ -27 +27 @@
 */
 
+#include <string>
 #include <utility>
 #include <vector>
@@ -61 +62 @@
 
   /**
-     SSearch does a rigorous Smith-Waterman (Needleman-Wunsch with no gap
-     cost) search for similarity between a sequnces. For long sequences this
-     may be very expensive (both in time and space) and BLAST or FASTA is
-     preferable.
+     \brief Local alignment following the Smith-Waterman
+     algorithm. 
 
-     @return number of matches
+     The original paper can be found here:
+     http://gel.ym.edu.tw/~chc/AB_papers/03.pdf
+
+     Instead of looking at each sequence in its entirety the S-W algorithm
+     compares segemnt of all possible lengths (LOCAL alignment) and chooses
+     whichever maximises the similarity measure.
+
+     \param s score matrix in which element $(i,j)$ is the score between
+     element $i$ in first sequence vs element $j$ in second sequence.
+     \param gap cost for having a gap (insertion or deletion)
+     \param open_gap cost for open up a gap in sequence, in other words, for a
+     gap of length $l$ the total cost is $open_gap + l*gap$.
    */
-  unsigned int ssearch(std::string first, std::string second);
+  double SmithWaterman(const utility::matrix& s,
+                       double gap, double open_gap);
+
+  /**
+     SSearch does a rigorous Smith-Waterman search for similarity between 
+     sequnces. For long sequences this may be very expensive (both in time and
+     space) and BLAST or FASTA is preferable.
+
+     @return score
+   */
+  double ssearch(std::string first, std::string second, double gap, 
+                 double open_gap);
 
 }}} // of namespace utility, yat, and theplu

trunk/yat/utility/stl_utility.h

@@ -35 +35 @@
 ///
 
+#include <algorithm>
 #include <ostream>
 #include <string>
@@ -55 +56 @@
 namespace yat {
 namespace utility {
+
+  /**
+     \return max of values
+   */
+  template <typename T>
+  T max(const T& a, const T& b, const T& c)
+  {
+    return std::max(std::max(a,b),c);
+  }
+
+
+  /**
+     \return max of values
+   */
+  template <typename T>
+  T max(const T& a, const T& b, const T& c, const T& d)
+  {
+    return std::max(std::max(a,b), std::max(c,d));
+  }
+
+
+  /**
+     \return max of values
+   */
+  template <typename T>
+  T max(const T& a, const T& b, const T& c, const T& d, const T& e)
+  {
+    return std::max(max(a,b,c,d), e);
+  }
+
+
+  /**
+     \return max of values
+   */
+  template <typename T>
+  T max(const T& a, const T& b, const T& c, const T& d, const T& e, const T& f)
+  {
+    return std::max(max(a,b,c,d), std::max(e,f));
+  }