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

Last change on this file since 3417 was 3417, checked in by Peter, 4 years ago

updating copyright statements

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 4.4 KB
Line 
1#ifndef _theplu_yat_utility_smith_waterman
2#define _theplu_yat_utility_smith_waterman
3
4// $Id: SmithWaterman.h 3417 2015-05-25 01:35:59Z peter $
5
6/*
7  Copyright (C) 2014, 2015 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 <boost/concept_check.hpp>
29
30#include <cstddef> // size_t
31
32namespace theplu {
33namespace yat {
34namespace utility {
35
36  /**
37     \since New in yat 0.12
38   */
39  class SmithWaterman
40  {
41  public:
42    /**
43       \brief Constructor
44
45       Same as SmithWaterman(gap, open_gap, gap, open_gap)
46     */
47    explicit SmithWaterman(double gap=0, double open_gap=0);
48
49    /**
50       \brief Constructor
51
52       \param vertical_gap Penalty for extending a vertical gap. A
53       vertical gap means consuming a reference element i.e. a
54       deletion in query sequence.
55
56       \param open_vertical_gap Open a deletion. Total
57       cost for a deletion is open_vertical_gap + N *
58       vertical_gap.
59
60       \param horizon_gap Penalty for extending a insertion.
61
62       \param open_horizon_gap Penalty for open an insertion. Total
63       penalty for a insertion is open_horizon_gap + N * horizon_gap.
64     */
65    SmithWaterman(double vertical_gap, double open_vertical_gap,
66                  double horizon_gap, double open_horizon_gap);
67
68    /**
69       The CIGAR here is slightly different compared to CIGAR in
70       Aligner as it uses BAM_CEQUAL and BAM_CDIFF rather than
71       BAM_CMATCH, and it also uses BAM_CSOFT_CLIP and ends of CIGAR
72       to reflect if whole query was aligned.
73
74       \return CIGAR reflecting latest performed alignment
75     */
76    const Aligner::Cigar& cigar(void) const;
77
78    /**
79       \return position i.e. first element in reference that is
80       aligned to the query.
81     */
82    size_t position(void) const;
83
84    /**
85       \return score matrix
86     */
87    const Matrix& score(void) const;
88
89    /**
90       Find best alignment using Smith-Waterman algorithm based on the
91       dot-matrix.
92
93       \return smith-waterman score
94     */
95    double operator()(const Matrix& dot);
96
97    /**
98       Calculate a dot-matrix of size query_end-query_begin X
99       reference_end-reference_begin where the element ij is 1.0 if
100       reference_begin[i] == query_begin[j] and -mismatch otherwise.
101
102       Type Requirements:
103       - \c RandomAccessIterator1 is a \readable_iterator
104       - \c RandomAccessIterator1 is a \random_access_traversal_iterator
105       - \c RandomAccessIterator2 is a \readable_iterator
106       - \c RandomAccessIterator2 is a \random_access_traversal_iterator
107     */
108    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
109    double operator()(RandomAccessIterator1 reference_begin,
110                      RandomAccessIterator1 reference_end,
111                      RandomAccessIterator2 query_begin,
112                      RandomAccessIterator2 query_end,
113                      double mismatch=0);
114
115  private:
116    Aligner aligner_;
117    Aligner::Cigar cigar_;
118    size_t position_;
119    Matrix score_;
120  };
121
122  // template implementaion
123
124  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
125  double SmithWaterman::operator()(RandomAccessIterator1 reference_begin,
126                                   RandomAccessIterator1 reference_end,
127                                   RandomAccessIterator2 query_begin,
128                                   RandomAccessIterator2 query_end,
129                                   double mismatch)
130  {
131    using boost_concepts::ReadableIterator;
132    using boost_concepts::RandomAccessTraversal;
133    BOOST_CONCEPT_ASSERT((ReadableIterator<RandomAccessIterator1>));
134    BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator1>));
135    BOOST_CONCEPT_ASSERT((ReadableIterator<RandomAccessIterator2>));
136    BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator2>));
137    Matrix dot(reference_end-reference_begin, query_end-query_begin, -mismatch);
138    for (size_t i=0; i<dot.rows(); ++i)
139      for (size_t j=0; j<dot.columns(); ++j)
140        if (reference_begin[i] == query_begin[j])
141          dot(i, j) = 1;
142    return this->operator()(dot);
143  }
144
145}}} // of namespace utility, yat, and theplu
146
147#endif
Note: See TracBrowser for help on using the repository browser.