yat/utility/SmithWaterman.h

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

Generated on Sat Oct 14 02:25:19 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)