Changeset 2815


Ignore:
Timestamp:
Aug 28, 2012, 4:38:36 AM (11 years ago)
Author:
Peter
Message:

closes #706. Implement a generalization of SmithWaternan? and NeedlemanWunsch? in a new class utility::Aligner. Class also provide a way to access not only alignment score but also which elements were aligned (to maximize the score). Old function SmithWaterman? did not provide this functionality and is deprecated for new member function Aligner::smith_waterman.

Location:
trunk
Files:
2 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEWS

    r2798 r2815  
    44
    55version 0.9 (released NOT YET)
     6  - utility::SmithWaterman is deprecated (ticket #706)
    67  - yat now depends on boost thread and system library (ticket #568 and #714)
    78  - utility::SmartPtr was marked internal (see r2778)

  • trunk/test/alignment.cc

    r2370 r2815  
    2424#include "Suite.h"
    2525
     26#include "yat/utility/Aligner.h"
    2627#include "yat/utility/Alignment.h"
    2728
    … …  
    112113    for (size_t k=0; k<b.size(); ++k)
    113114      m(j,k) = (a[j]==b[k] ? 1 : -1000);
    114   double score=utility::SmithWaterman(m, 100, 100);
     115  double score = utility::Aligner(100, 100).smith_waterman(m);
    115116  if (score!=21)
    116117    suite.add(false);

  • trunk/yat/utility/Alignment.cc

    r2792 r2815  
    2525
    2626#include "Alignment.h"
     27
     28#include "Aligner.h"
    2729#include "Matrix.h"
    2830#include "stl_utility.h"
    … …  
    3638namespace utility {
    3739
    38   double NeedlemanWunsch(const utility::Matrix& s,
     40  double NeedlemanWunsch(const Matrix& s,
    3941                         std::vector<std::pair<size_t, size_t> >& path,
    4042                         const double gap)
    4143  {
    42     utility::Matrix m(s.rows()+1,s.columns()+1);
    43     // Init upper and left border of matrix
    44     for (size_t i=1; i<m.rows(); i++)
    45       m(i,0)=-i*gap;
    46     for (size_t i=1; i<m.columns(); i++)
    47       m(0,i)=-i*gap;
    48     // choice(i,j) tells us how we came to s(i,j). 1 is diagonal, 2
    49     // vertical, and 3 horizontal,
    50     utility::Matrix choice(m.rows(),m.columns());
     44    // align
     45    Aligner aligner(gap);
     46    double score = aligner.needleman_wunsch(s);
    5147
    52     // Calculating NeedlemanWunsch matrix
    53     for (size_t i=1; i<m.rows(); i++)
    54       for (size_t j=1; j<m.columns(); j++){
    55         if (m(i-1,j-1) + s(i-1,j-1) > m(i-1,j)-gap &&
    56             m(i-1,j-1) + s(i-1,j-1) > m(i,j-1)-gap){
    57           m(i,j)=m(i-1,j-1) + s(i-1,j-1);
    58           choice(i,j)=1;
    59         }
    60         else if (m(i-1,j) > m(i,j-1)){
    61           m(i,j)=m(i-1,j)-gap;
    62           choice(i,j)=2;
    63         }
    64         else{
    65           m(i,j)=m(i,j-1)-gap;
    66           choice(i,j)=3;
    67         }
    68       }
    6948    // Going backwards to find best path
    70     size_t i = m.rows()-1;
    71     size_t j= m.columns()-1;
    72     path.resize(0);
    73     while (i && j){
    74       if (choice(i,j)==1){
    75         i--;
    76         j--;
     49    size_t i = s.rows();
     50    size_t j= s.columns();
     51    path.clear();
     52    while (i && j) {
     53      if (aligner.cigar(i,j)==Aligner::diagonal) {
     54        --i;
     55        --j;
    7756        path.push_back(std::make_pair(i,j));
    7857      }
    79       else if (choice(i,j)==2)
    80         i--;
     58      else if (aligner.cigar(i,j)==Aligner::right)
     59        --j;
    8160      else
    82         j--;
     61        --i;
    8362    }
    84 
    85     return m(m.rows()-1,m.columns()-1);
     63    return score;
    8664  }
    8765
    … …  
    9674      for (size_t j=0; j<second.size(); ++j)
    9775        m(i,j) = (first[i]==second[j] ? 1 : -mismatch);
    98     return SmithWaterman(m, gap, open_gap);
     76    Aligner aligner(gap, open_gap);
     77    return aligner.smith_waterman(m);
    9978  }
    10079
    … …  
    10382                       double gap, double open_gap)
    10483  {
    105     enum direction { none, right, down, diagonal };
    106 
    107     // Calculating S-W matrix
    108     Matrix m(s.rows()+1,s.columns()+1);
    109     Matrix array(m);
    110     for (size_t i=1; i<m.rows(); ++i)
    111       for (size_t j=1; j<m.columns(); ++j){
    112         m(i,j)=0;
    113         array(i,j) = none;
    114         if (m(i-1,j-1)+s(i-1,j-1)>m(i,j)){
    115           m(i,j) = m(i-1,j-1)+s(i-1,j-1);
    116           array(i,j) = diagonal;
    117         }
    118         if (array(i-1,j)==right && m(i-1,j)-gap > m(i,j)){
    119           m(i,j) = m(i-1,j)-gap;
    120           array(i,j) = right;
    121         }
    122         if (array(i-1,j)!=right && m(i-1,j)-gap-open_gap > m(i,j)){
    123           m(i,j) = m(i-1,j)-gap-open_gap;
    124           array(i,j) = right;
    125         }
    126         if (array(i,j-1)==down && m(i,j-1)-gap > m(i,j)){
    127           m(i,j) = m(i,j-1)-gap;
    128           array(i,j) = down;
    129         }
    130         if (array(i,j-1)!=down && m(i,j-1)-gap-open_gap > m(i,j)){
    131           m(i,j) = m(i,j-1)-gap-open_gap;
    132           array(i,j) = down;
    133         }
    134       }
    135     return max(m);
     84    Aligner aligner(gap, open_gap);
     85    return aligner.smith_waterman(s);
    13686  }
    13787

  • trunk/yat/utility/Alignment.h

    r2792 r2815  
    2626  along with yat. If not, see <http://www.gnu.org/licenses/>.
    2727*/
     28
     29#include "deprecate.h"
    2830
    2931#include <string>
    … …  
    5759  /// @return the global maximum alignment score.
    5860  ///
     61  /// \see Aligner
     62  ///
    5963  double NeedlemanWunsch(const utility::Matrix& s,
    6064                         std::vector<std::pair<size_t, size_t> >& path,
    … …  
    8387     \param open_gap cost for open up a gap in sequence, in other words, for a
    8488     gap of length L the total cost is open_gap + L*gap.
     89
     90     \deprecated Provided for backward compatibility with 0.8 API. Use
     91     Aligner::smith_waterman() instead.
    8592   */
    8693  double SmithWaterman(const utility::Matrix& s,
    87                        double gap, double open_gap);
     94                       double gap, double open_gap) YAT_DEPRECATE;
    8895
    8996  /**
    … …  
    99106
    100107     \since parameter mismatch available since yat 0.9
     108
     109     \see Aligner
    101110   */
    102111  double ssearch(std::string first, std::string second, double gap,

  • trunk/yat/utility/Makefile.am

    r2787 r2815  
    2323noinst_LTLIBRARIES += yat/utility/libutility.la
    2424yat_utility_libutility_la_SOURCES = \
     25  yat/utility/Aligner.cc \
    2526  yat/utility/Alignment.cc \
    2627  yat/utility/ColumnStream.cc \
    … …  
    5859
    5960nobase_include_HEADERS += \
     61  $(srcdir)/yat/utility/Aligner.h \
    6062  $(srcdir)/yat/utility/Alignment.h \
    6163  $(srcdir)/yat/utility/ColumnStream.h \
Note: See TracChangeset for help on using the changeset viewer.