yat/utility/SmithWaterman.cc

Code
Comments
Other
Rev Date Author Line
3204 03 May 14 peter 1 // $Id$
3204 03 May 14 peter 2
3204 03 May 14 peter 3 /*
4207 26 Aug 22 peter 4   Copyright (C) 2014, 2022 Peter Johansson
3204 03 May 14 peter 5
3204 03 May 14 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
3204 03 May 14 peter 7
3204 03 May 14 peter 8   The yat library is free software; you can redistribute it and/or
3204 03 May 14 peter 9   modify it under the terms of the GNU General Public License as
3204 03 May 14 peter 10   published by the Free Software Foundation; either version 3 of the
3204 03 May 14 peter 11   License, or (at your option) any later version.
3204 03 May 14 peter 12
3204 03 May 14 peter 13   The yat library is distributed in the hope that it will be useful,
3204 03 May 14 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
3204 03 May 14 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3204 03 May 14 peter 16   General Public License for more details.
3204 03 May 14 peter 17
3204 03 May 14 peter 18   You should have received a copy of the GNU General Public License
3204 03 May 14 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
3204 03 May 14 peter 20 */
3204 03 May 14 peter 21
3204 03 May 14 peter 22 #include <config.h>
3204 03 May 14 peter 23
3204 03 May 14 peter 24 #include "SmithWaterman.h"
3212 05 May 14 peter 25 #include "Cigar.h"
3204 03 May 14 peter 26
3204 03 May 14 peter 27 #include <cassert>
3204 03 May 14 peter 28
3204 03 May 14 peter 29 namespace theplu {
3204 03 May 14 peter 30 namespace yat {
3204 03 May 14 peter 31 namespace utility {
3204 03 May 14 peter 32
3205 04 May 14 peter 33   SmithWaterman::SmithWaterman(double gap, double open_gap)
3205 04 May 14 peter 34     : aligner_(gap, open_gap) {}
3205 04 May 14 peter 35
3205 04 May 14 peter 36
3205 04 May 14 peter 37   SmithWaterman::SmithWaterman(double vertical_gap, double open_vertical_gap,
3205 04 May 14 peter 38                                double horizon_gap, double open_horizon_gap)
3205 04 May 14 peter 39     : aligner_(vertical_gap, open_vertical_gap, horizon_gap, open_horizon_gap)
3205 04 May 14 peter 40   {}
3205 04 May 14 peter 41
3205 04 May 14 peter 42
3205 04 May 14 peter 43   const Aligner::Cigar& SmithWaterman::cigar(void) const
3205 04 May 14 peter 44   {
3205 04 May 14 peter 45     return cigar_;
3205 04 May 14 peter 46   }
3205 04 May 14 peter 47
3205 04 May 14 peter 48
3205 04 May 14 peter 49   size_t SmithWaterman::position(void) const
3205 04 May 14 peter 50   {
3205 04 May 14 peter 51     return position_;
3205 04 May 14 peter 52   }
3205 04 May 14 peter 53
3205 04 May 14 peter 54
3205 04 May 14 peter 55   const Matrix& SmithWaterman::score(void) const
3205 04 May 14 peter 56   {
3205 04 May 14 peter 57     return score_;
3205 04 May 14 peter 58   }
3205 04 May 14 peter 59
3205 04 May 14 peter 60
4125 14 Jan 22 peter 61   double SmithWaterman::operator()(const MatrixBase& dot)
3205 04 May 14 peter 62   {
3205 04 May 14 peter 63     score_.resize(dot.rows()+1, dot.columns()+1);
3205 04 May 14 peter 64     aligner_(dot, score_);
3205 04 May 14 peter 65     cigar_.clear();
3205 04 May 14 peter 66     size_t row = 0;
3205 04 May 14 peter 67     size_t column = 0;
3205 04 May 14 peter 68     for (size_t i=0; i<score_.rows(); ++i)
3205 04 May 14 peter 69       for (size_t j=0; j<score_.columns(); ++j)
3205 04 May 14 peter 70         if (score_(i, j) > score_(row, column)) {
3205 04 May 14 peter 71           row = i;
3205 04 May 14 peter 72           column = j;
3205 04 May 14 peter 73         }
3205 04 May 14 peter 74
3205 04 May 14 peter 75     cigar_.push_front(BAM_CSOFT_CLIP, score_.columns()-column-1);
3205 04 May 14 peter 76
3205 04 May 14 peter 77     size_t end_row(row);
3205 04 May 14 peter 78     size_t end_column(column);
3205 04 May 14 peter 79
3205 04 May 14 peter 80     while (aligner_.alignment(row, column)!=Aligner::none) {
3205 04 May 14 peter 81       switch (aligner_.alignment(row, column)) {
3205 04 May 14 peter 82       case(Aligner::diagonal):
3205 04 May 14 peter 83         assert(row);
3205 04 May 14 peter 84         --row;
3205 04 May 14 peter 85         assert(column);
3205 04 May 14 peter 86         --column;
3205 04 May 14 peter 87         if (dot(row, column)>0)
3205 04 May 14 peter 88           cigar_.push_front(BAM_CEQUAL);
3205 04 May 14 peter 89         else
3205 04 May 14 peter 90           cigar_.push_front(BAM_CDIFF);
3205 04 May 14 peter 91         break;
3205 04 May 14 peter 92       case (Aligner::right):
3205 04 May 14 peter 93         assert(column);
3205 04 May 14 peter 94         --column;
3205 04 May 14 peter 95         cigar_.push_front(BAM_CINS);
3205 04 May 14 peter 96         break;
3205 04 May 14 peter 97       case (Aligner::down):
3205 04 May 14 peter 98         assert(row);
3205 04 May 14 peter 99         --row;
3205 04 May 14 peter 100         cigar_.push_front(BAM_CDEL);
3205 04 May 14 peter 101         break;
3205 04 May 14 peter 102       default:;
3205 04 May 14 peter 103       }
3205 04 May 14 peter 104     }
3205 04 May 14 peter 105
3205 04 May 14 peter 106     position_ = row;
3205 04 May 14 peter 107     cigar_.push_front(BAM_CSOFT_CLIP, column);
3205 04 May 14 peter 108
3205 04 May 14 peter 109     return score_(end_row, end_column);
3205 04 May 14 peter 110   }
3205 04 May 14 peter 111
3204 03 May 14 peter 112 }}} // of namespace utility, yat, and theplu

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