yat/utility/NeedlemanWunsch.cc

Code
Comments
Other
Rev Date Author Line
4334 25 Mar 23 peter 1 // $Id$
4334 25 Mar 23 peter 2
4334 25 Mar 23 peter 3 /*
4334 25 Mar 23 peter 4   Copyright (C) 2023 Peter Johansson
4334 25 Mar 23 peter 5
4334 25 Mar 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4334 25 Mar 23 peter 7
4334 25 Mar 23 peter 8   The yat library is free software; you can redistribute it and/or
4334 25 Mar 23 peter 9   modify it under the terms of the GNU General Public License as
4334 25 Mar 23 peter 10   published by the Free Software Foundation; either version 3 of the
4334 25 Mar 23 peter 11   License, or (at your option) any later version.
4334 25 Mar 23 peter 12
4334 25 Mar 23 peter 13   The yat library is distributed in the hope that it will be useful,
4334 25 Mar 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4334 25 Mar 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4334 25 Mar 23 peter 16   General Public License for more details.
4334 25 Mar 23 peter 17
4334 25 Mar 23 peter 18   You should have received a copy of the GNU General Public License
4334 25 Mar 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4334 25 Mar 23 peter 20 */
4334 25 Mar 23 peter 21
4334 25 Mar 23 peter 22 #include <config.h>
4334 25 Mar 23 peter 23
4334 25 Mar 23 peter 24 #include "NeedlemanWunsch.h"
4334 25 Mar 23 peter 25 #include "Cigar.h"
4334 25 Mar 23 peter 26
4334 25 Mar 23 peter 27 #include <cassert>
4334 25 Mar 23 peter 28
4334 25 Mar 23 peter 29 namespace theplu {
4334 25 Mar 23 peter 30 namespace yat {
4334 25 Mar 23 peter 31 namespace utility {
4334 25 Mar 23 peter 32
4334 25 Mar 23 peter 33   NeedlemanWunsch::NeedlemanWunsch(double gap, double open_gap)
4334 25 Mar 23 peter 34     : aligner_(gap, open_gap) {}
4334 25 Mar 23 peter 35
4334 25 Mar 23 peter 36
4334 25 Mar 23 peter 37   NeedlemanWunsch::NeedlemanWunsch(double vertical_gap, double open_vertical_gap,
4334 25 Mar 23 peter 38                                double horizon_gap, double open_horizon_gap)
4334 25 Mar 23 peter 39     : aligner_(vertical_gap, open_vertical_gap, horizon_gap, open_horizon_gap)
4334 25 Mar 23 peter 40   {}
4334 25 Mar 23 peter 41
4334 25 Mar 23 peter 42
4334 25 Mar 23 peter 43   const Aligner::Cigar& NeedlemanWunsch::cigar(void) const
4334 25 Mar 23 peter 44   {
4334 25 Mar 23 peter 45     return cigar_;
4334 25 Mar 23 peter 46   }
4334 25 Mar 23 peter 47
4334 25 Mar 23 peter 48
4334 25 Mar 23 peter 49   double NeedlemanWunsch::operator()(const MatrixBase& dot)
4334 25 Mar 23 peter 50   {
4334 25 Mar 23 peter 51     double s = aligner_.needleman_wunsch(dot);
4334 25 Mar 23 peter 52     auto cigar = aligner_.cigar(dot.rows(), dot.columns());
4334 25 Mar 23 peter 53     uint32_t ref = dot.rows() - cigar.reference_length();
4334 25 Mar 23 peter 54     cigar_.clear();
4334 25 Mar 23 peter 55     cigar_.push_back(BAM_CDEL, ref);
4334 25 Mar 23 peter 56     uint32_t query=0;
4334 25 Mar 23 peter 57     for (size_t i=0; i<cigar.size(); ++i) {
4334 25 Mar 23 peter 58       switch (cigar.op(i)) {
4334 25 Mar 23 peter 59       case (BAM_CMATCH) :
4334 25 Mar 23 peter 60         for (size_t j=0; j<cigar.oplen(i); ++j) {
4334 25 Mar 23 peter 61           cigar_.push_back(dot(ref, query) > 0.0 ? BAM_CEQUAL : BAM_CDIFF);
4334 25 Mar 23 peter 62           ++ref;
4334 25 Mar 23 peter 63           ++query;
4334 25 Mar 23 peter 64         }
4334 25 Mar 23 peter 65         break;
4334 25 Mar 23 peter 66       case (BAM_CDEL) :
4334 25 Mar 23 peter 67         cigar_.push_back(cigar.op(i), cigar.oplen(i));
4334 25 Mar 23 peter 68         ref += cigar.oplen(i);
4334 25 Mar 23 peter 69         break;
4334 25 Mar 23 peter 70       case (BAM_CINS) :
4334 25 Mar 23 peter 71         cigar_.push_back(cigar.op(i), cigar.oplen(i));
4334 25 Mar 23 peter 72         query += cigar.oplen(i);
4334 25 Mar 23 peter 73         break;
4334 25 Mar 23 peter 74       default:
4334 25 Mar 23 peter 75         throw runtime_error("NeedlemanWunsch invalid cigar element");
4334 25 Mar 23 peter 76       }
4334 25 Mar 23 peter 77     }
4334 25 Mar 23 peter 78     assert(cigar_.reference_length() == dot.rows());
4334 25 Mar 23 peter 79     assert(cigar_.query_length() == dot.columns());
4334 25 Mar 23 peter 80     return s;
4334 25 Mar 23 peter 81   }
4334 25 Mar 23 peter 82
4334 25 Mar 23 peter 83 }}} // of namespace utility, yat, and theplu

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