Ignore:
Timestamp:
Sep 14, 2007, 8:59:46 PM (14 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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/yat/utility/Alignment.cc

    r867 r869  
    2626#include "Alignment.h"
    2727#include "matrix.h"
     28#include "yat/utility/stl_utility.h"
    2829
     30#include <algorithm>
    2931#include <utility>
    3032#include <vector>
     
    3335namespace yat {
    3436namespace utility {
    35 
    36   unsigned int ssearch(std::string first, std::string second)
    37   {
    38     matrix m(first.size(), second.size());
    39     for (size_t i=0; i<first.size(); ++i)
    40       for (size_t j=0; j<second.size(); ++j)
    41         m(i,j) = (first[i]==second[j] ? 1 : 0);
    42     std::vector<std::pair<size_t, size_t> > path;
    43     return static_cast<unsigned int>(NeedlemanWunsch(m,path,0));
    44   }
    45 
    4637
    4738  double NeedlemanWunsch(const utility::matrix& s,
     
    9586  }
    9687
     88
     89  double ssearch(std::string first, std::string second, double gap,
     90                 double open_gap)
     91  {
     92    matrix m(first.size(), second.size());
     93    for (size_t i=0; i<first.size(); ++i)
     94      for (size_t j=0; j<second.size(); ++j)
     95        m(i,j) = (first[i]==second[j] ? 1 : 0);
     96    return SmithWaterman(m, gap, open_gap);
     97  }
     98
     99
     100  double SmithWaterman(const utility::matrix& s,
     101                       double gap, double open_gap)
     102  {
     103    // Calculating S-W matrix
     104    matrix m(s.rows()+2,s.columns()+2);
     105    for (size_t i=2; i<m.rows(); ++i)
     106      for (size_t j=2; j<m.columns(); ++j)
     107        m(i,j) = max(0.0, m(i-1,j-1)+s(i-2,j-2),
     108                     m(i-1,j)-gap-open_gap, m(i-2,j)-2*gap,
     109                     m(i,j-1)-gap-open_gap, m(i,j-2)-2*gap);
     110
     111    return max(m);
     112  }
     113
    97114}}} // of namespace utility, yat, and theplu
Note: See TracChangeset for help on using the changeset viewer.