yat/utility/NeedlemanWunsch.h

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

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