source: trunk/yat/utility/SmithWaterman.h @ 3205

Last change on this file since 3205 was 3205, checked in by Peter, 8 years ago

implement new class SmithWaterman?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.7 KB
Line 
1#ifndef _theplu_yat_utility_smith_waterman
2#define _theplu_yat_utility_smith_waterman
3
4// $Id: SmithWaterman.h 3205 2014-05-04 02:50:50Z peter $
5
6/*
7  Copyright (C) 2014 Peter Johansson
8
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with yat. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include "Aligner.h"
26#include "Matrix.h"
27
28#include <cstddef> // size_t
29
30namespace theplu {
31namespace yat {
32namespace utility {
33
34  /**
35     \since New in yat 0.12
36   */
37  class SmithWaterman
38  {
39  public:
40    /**
41       \brief Constructor
42
43       Same as SmithWaterman(gap, open_gap, gap, open_gap)
44     */
45    explicit SmithWaterman(double gap=0, double open_gap=0);
46
47    /**
48       \brief Constructor
49
50       \param vertical_gap Penalty for extending a vertical gap. A
51       vertical gap means consuming a reference element i.e. a
52       deletion in query sequence.
53
54       \param open_vertical_gap Open a deletion. Total
55       cost for a deletion is open_vertical_gap + N *
56       vertical_gap.
57
58       \param horizon_gap Penalty for extending a insertion.
59
60       \param open_horizon_gap Penalty for open an insertion. Total
61       penalty for a insertion is open_horizon_gap + N * horizon_gap.
62     */
63    SmithWaterman(double vertical_gap, double open_vertical_gap,
64                  double horizon_gap, double open_horizon_gap);
65
66    /**
67       The CIGAR here is slightly different compared to CIGAR in
68       Aligner as it uses BAM_CEQUAL and BAM_CDIFF rather than
69       BAM_CMATCH, and it also uses BAM_CSOFT_CLIP and ends of CIGAR
70       to reflect if whole query was aligned.
71
72       \return CIGAR reflecting latest performed alignment
73     */
74    const Aligner::Cigar& cigar(void) const;
75
76    /**
77       \return position i.e. first element in reference that is
78       aligned to the query.
79     */
80    size_t position(void) const;
81
82    /**
83       \return score matrix
84     */
85    const Matrix& score(void) const;
86
87    /**
88       Find best alignment using Smith-Waterman algorithm based on the
89       dot-matrix.
90
91       \return smith-waterman score
92     */
93    double operator()(const Matrix& dot);
94
95    /**
96       Calculate a dot-matrix of size query_end-query_begin X
97       reference_end-reference_begin where the element ij is 1.0 if
98       query_begin[i] == reference_begin[j] and -mismatch otherwise.
99     */
100    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
101    double operator()(RandomAccessIterator1 query_begin,
102                      RandomAccessIterator1 query_end,
103                      RandomAccessIterator2 reference_begin,
104                      RandomAccessIterator2 reference_end,
105                      double mismatch=0);
106
107  private:
108    Aligner aligner_;
109    Aligner::Cigar cigar_;
110    size_t position_;
111    Matrix score_;
112  };
113
114  // template implementaion
115
116  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
117  double SmithWaterman::operator()(RandomAccessIterator1 reference_begin,
118                                   RandomAccessIterator1 reference_end,
119                                   RandomAccessIterator2 query_begin,
120                                   RandomAccessIterator2 query_end,
121                                   double mismatch)
122  {
123    Matrix dot(reference_end-reference_begin, query_end-query_begin, -mismatch);
124    for (size_t i=0; i<dot.rows(); ++i)
125      for (size_t j=0; j<dot.columns(); ++j)
126        if (reference_begin[i] == query_begin[j])
127          dot(i, j) = 1;
128    return this->operator()(dot);
129  }
130
131}}} // of namespace utility, yat, and theplu
132
133#endif
Note: See TracBrowser for help on using the repository browser.