Changeset 869 for trunk/yat


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.

Location:
trunk/yat/utility
Files:
3 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
  • trunk/yat/utility/Alignment.h

    r867 r869  
    2727*/
    2828
     29#include <string>
    2930#include <utility>
    3031#include <vector>
     
    6162
    6263  /**
    63      SSearch does a rigorous Smith-Waterman (Needleman-Wunsch with no gap
    64      cost) search for similarity between a sequnces. For long sequences this
    65      may be very expensive (both in time and space) and BLAST or FASTA is
    66      preferable.
     64     \brief Local alignment following the Smith-Waterman
     65     algorithm.
    6766
    68      @return number of matches
     67     The original paper can be found here:
     68     http://gel.ym.edu.tw/~chc/AB_papers/03.pdf
     69
     70     Instead of looking at each sequence in its entirety the S-W algorithm
     71     compares segemnt of all possible lengths (LOCAL alignment) and chooses
     72     whichever maximises the similarity measure.
     73
     74     \param s score matrix in which element $(i,j)$ is the score between
     75     element $i$ in first sequence vs element $j$ in second sequence.
     76     \param gap cost for having a gap (insertion or deletion)
     77     \param open_gap cost for open up a gap in sequence, in other words, for a
     78     gap of length $l$ the total cost is $open_gap + l*gap$.
    6979   */
    70   unsigned int ssearch(std::string first, std::string second);
     80  double SmithWaterman(const utility::matrix& s,
     81                       double gap, double open_gap);
     82
     83  /**
     84     SSearch does a rigorous Smith-Waterman search for similarity between
     85     sequnces. For long sequences this may be very expensive (both in time and
     86     space) and BLAST or FASTA is preferable.
     87
     88     @return score
     89   */
     90  double ssearch(std::string first, std::string second, double gap,
     91                 double open_gap);
    7192
    7293}}} // of namespace utility, yat, and theplu
  • trunk/yat/utility/stl_utility.h

    r865 r869  
    3535///
    3636
     37#include <algorithm>
    3738#include <ostream>
    3839#include <string>
     
    5556namespace yat {
    5657namespace utility {
     58
     59  /**
     60     \return max of values
     61   */
     62  template <typename T>
     63  T max(const T& a, const T& b, const T& c)
     64  {
     65    return std::max(std::max(a,b),c);
     66  }
     67
     68
     69  /**
     70     \return max of values
     71   */
     72  template <typename T>
     73  T max(const T& a, const T& b, const T& c, const T& d)
     74  {
     75    return std::max(std::max(a,b), std::max(c,d));
     76  }
     77
     78
     79  /**
     80     \return max of values
     81   */
     82  template <typename T>
     83  T max(const T& a, const T& b, const T& c, const T& d, const T& e)
     84  {
     85    return std::max(max(a,b,c,d), e);
     86  }
     87
     88
     89  /**
     90     \return max of values
     91   */
     92  template <typename T>
     93  T max(const T& a, const T& b, const T& c, const T& d, const T& e, const T& f)
     94  {
     95    return std::max(max(a,b,c,d), std::max(e,f));
     96  }
    5797
    5898
Note: See TracChangeset for help on using the changeset viewer.