yat/utility/mSmithWaterman2.cc

Code
Comments
Other
Rev Date Author Line
4353 23 Aug 23 peter 1 // $Id$
4353 23 Aug 23 peter 2
4353 23 Aug 23 peter 3 /*
4353 23 Aug 23 peter 4   Copyright (C) 2023 Peter Johansson
4353 23 Aug 23 peter 5
4353 23 Aug 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4353 23 Aug 23 peter 7
4353 23 Aug 23 peter 8   The yat library is free software; you can redistribute it and/or
4353 23 Aug 23 peter 9   modify it under the terms of the GNU General Public License as
4353 23 Aug 23 peter 10   published by the Free Software Foundation; either version 3 of the
4353 23 Aug 23 peter 11   License, or (at your option) any later version.
4353 23 Aug 23 peter 12
4353 23 Aug 23 peter 13   The yat library is distributed in the hope that it will be useful,
4353 23 Aug 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4353 23 Aug 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4353 23 Aug 23 peter 16   General Public License for more details.
4353 23 Aug 23 peter 17
4353 23 Aug 23 peter 18   You should have received a copy of the GNU General Public License
4353 23 Aug 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4353 23 Aug 23 peter 20 */
4353 23 Aug 23 peter 21
4353 23 Aug 23 peter 22 #include <config.h>
4353 23 Aug 23 peter 23
4353 23 Aug 23 peter 24 #include "mSmithWaterman2.h"
4353 23 Aug 23 peter 25
4353 23 Aug 23 peter 26 #include <cassert>
4353 23 Aug 23 peter 27
4353 23 Aug 23 peter 28 namespace theplu {
4353 23 Aug 23 peter 29 namespace yat {
4353 23 Aug 23 peter 30 namespace utility {
4353 23 Aug 23 peter 31
4353 23 Aug 23 peter 32   mSmithWaterman2::mSmithWaterman2(double gap, double open_gap)
4353 23 Aug 23 peter 33     : mSmithWaterman2(gap, open_gap, gap, open_gap)
4353 23 Aug 23 peter 34   {}
4353 23 Aug 23 peter 35
4353 23 Aug 23 peter 36
4353 23 Aug 23 peter 37   mSmithWaterman2::mSmithWaterman2(double insertion_penalty,
4353 23 Aug 23 peter 38                                    double open_insertion_penalty,
4353 23 Aug 23 peter 39                                    double deletion_penalty,
4353 23 Aug 23 peter 40                                    double open_deletion_penalty)
4353 23 Aug 23 peter 41     : mSmithWaterman2(insertion_penalty, open_insertion_penalty,
4353 23 Aug 23 peter 42                       deletion_penalty, open_deletion_penalty,
4353 23 Aug 23 peter 43                       insertion_penalty, open_insertion_penalty,
4353 23 Aug 23 peter 44                       deletion_penalty, open_deletion_penalty)
4353 23 Aug 23 peter 45   {}
4353 23 Aug 23 peter 46
4353 23 Aug 23 peter 47
4353 23 Aug 23 peter 48   mSmithWaterman2::mSmithWaterman2(double insertion_penalty1,
4353 23 Aug 23 peter 49                                    double open_insertion_penalty1,
4353 23 Aug 23 peter 50                                    double deletion_penalty1,
4353 23 Aug 23 peter 51                                    double open_deletion_penalty1,
4353 23 Aug 23 peter 52                                    double insertion_penalty2,
4353 23 Aug 23 peter 53                                    double open_insertion_penalty2,
4353 23 Aug 23 peter 54                                    double deletion_penalty2,
4353 23 Aug 23 peter 55                                    double open_deletion_penalty2)
4353 23 Aug 23 peter 56   : mSmithWaterman2(insertion_penalty1, open_insertion_penalty1,
4353 23 Aug 23 peter 57                     deletion_penalty1, open_deletion_penalty1,
4353 23 Aug 23 peter 58                     insertion_penalty2, open_insertion_penalty2,
4353 23 Aug 23 peter 59                     deletion_penalty2, open_deletion_penalty2, 0, 0)
4353 23 Aug 23 peter 60   {}
4353 23 Aug 23 peter 61
4353 23 Aug 23 peter 62
4353 23 Aug 23 peter 63   mSmithWaterman2::mSmithWaterman2(double insertion_penalty1,
4353 23 Aug 23 peter 64                                    double open_insertion_penalty1,
4353 23 Aug 23 peter 65                                    double deletion_penalty1,
4353 23 Aug 23 peter 66                                    double open_deletion_penalty1,
4353 23 Aug 23 peter 67                                    double insertion_penalty2,
4353 23 Aug 23 peter 68                                    double open_insertion_penalty2,
4353 23 Aug 23 peter 69                                    double deletion_penalty2,
4353 23 Aug 23 peter 70                                    double open_deletion_penalty2,
4353 23 Aug 23 peter 71                                    double gap_penalty,
4353 23 Aug 23 peter 72                                    double open_gap_penalty)
4353 23 Aug 23 peter 73     : insertion_penalty1_(insertion_penalty1),
4353 23 Aug 23 peter 74       open_insertion_penalty1_(open_insertion_penalty1),
4353 23 Aug 23 peter 75       deletion_penalty1_(deletion_penalty1),
4353 23 Aug 23 peter 76       open_deletion_penalty1_(open_deletion_penalty1),
4353 23 Aug 23 peter 77       insertion_penalty2_(insertion_penalty2),
4353 23 Aug 23 peter 78       open_insertion_penalty2_(open_insertion_penalty2),
4353 23 Aug 23 peter 79       deletion_penalty2_(deletion_penalty2),
4353 23 Aug 23 peter 80       open_deletion_penalty2_(open_deletion_penalty2),
4353 23 Aug 23 peter 81       gap_penalty_(gap_penalty), open_gap_penalty_(open_gap_penalty)
4353 23 Aug 23 peter 82   {
4353 23 Aug 23 peter 83   }
4353 23 Aug 23 peter 84
4353 23 Aug 23 peter 85
4353 23 Aug 23 peter 86   const mSmithWaterman2::Result& mSmithWaterman2::first(void) const
4353 23 Aug 23 peter 87   {
4353 23 Aug 23 peter 88     return first_;
4353 23 Aug 23 peter 89   }
4353 23 Aug 23 peter 90
4353 23 Aug 23 peter 91
4353 23 Aug 23 peter 92   const mSmithWaterman2::Result& mSmithWaterman2::second(void) const
4353 23 Aug 23 peter 93   {
4353 23 Aug 23 peter 94     return second_;
4353 23 Aug 23 peter 95   }
4353 23 Aug 23 peter 96
4353 23 Aug 23 peter 97
4353 23 Aug 23 peter 98   size_t mSmithWaterman2::gap(void) const
4353 23 Aug 23 peter 99   {
4353 23 Aug 23 peter 100     return gap_;
4353 23 Aug 23 peter 101   }
4353 23 Aug 23 peter 102
4353 23 Aug 23 peter 103
4353 23 Aug 23 peter 104   double mSmithWaterman2::score(void) const
4353 23 Aug 23 peter 105   {
4353 23 Aug 23 peter 106     return score_;
4353 23 Aug 23 peter 107   }
4353 23 Aug 23 peter 108
4353 23 Aug 23 peter 109
4353 23 Aug 23 peter 110   // mSmithWaterman2::Result
4353 23 Aug 23 peter 111
4353 23 Aug 23 peter 112   size_t mSmithWaterman2::Result::position(void) const
4353 23 Aug 23 peter 113   {
4353 23 Aug 23 peter 114     return position_;
4353 23 Aug 23 peter 115   }
4353 23 Aug 23 peter 116
4353 23 Aug 23 peter 117
4353 23 Aug 23 peter 118   const Aligner::Cigar& mSmithWaterman2::Result::cigar(void) const
4353 23 Aug 23 peter 119   {
4353 23 Aug 23 peter 120     return cigar_;
4353 23 Aug 23 peter 121   }
4353 23 Aug 23 peter 122
4353 23 Aug 23 peter 123
4353 23 Aug 23 peter 124   const Aligner::Cigar& mSmithWaterman2::Result::cigar2(void) const
4353 23 Aug 23 peter 125   {
4353 23 Aug 23 peter 126     return cigar2_;
4353 23 Aug 23 peter 127   }
4353 23 Aug 23 peter 128
4353 23 Aug 23 peter 129
4353 23 Aug 23 peter 130   double mSmithWaterman2::Result::score(void) const
4353 23 Aug 23 peter 131   {
4353 23 Aug 23 peter 132     return score_;
4353 23 Aug 23 peter 133   }
4353 23 Aug 23 peter 134
4353 23 Aug 23 peter 135
4353 23 Aug 23 peter 136   // Aligner 1
4353 23 Aug 23 peter 137   mSmithWaterman2::Aligner1::Aligner1(double xinsertion, double insertion,
4353 23 Aug 23 peter 138                                       double xdeletion, double deletion)
4353 23 Aug 23 peter 139     : mAligner<mSmithWaterman2::Aligner1>(mAlignerBase::Affine(xinsertion,
4353 23 Aug 23 peter 140                                                                insertion),
4353 23 Aug 23 peter 141                                           mAlignerBase::Affine(xdeletion,
4353 23 Aug 23 peter 142                                                                deletion))
4353 23 Aug 23 peter 143   {}
4353 23 Aug 23 peter 144
4353 23 Aug 23 peter 145
4353 23 Aug 23 peter 146   const std::vector<mAlignerBase::Alignment>&
4353 23 Aug 23 peter 147   mSmithWaterman2::Aligner1::col_best(void) const
4353 23 Aug 23 peter 148   {
4353 23 Aug 23 peter 149     return col_best_;
4353 23 Aug 23 peter 150   }
4353 23 Aug 23 peter 151
4353 23 Aug 23 peter 152
4353 23 Aug 23 peter 153   double mSmithWaterman2::Aligner1::init_row(size_t q) const
4353 23 Aug 23 peter 154   {
4353 23 Aug 23 peter 155     return 0;
4353 23 Aug 23 peter 156   }
4353 23 Aug 23 peter 157
4353 23 Aug 23 peter 158
4353 23 Aug 23 peter 159   double mSmithWaterman2::Aligner1::init_col(size_t r) const
4353 23 Aug 23 peter 160   {
4353 23 Aug 23 peter 161     return 0;
4353 23 Aug 23 peter 162   }
4353 23 Aug 23 peter 163
4353 23 Aug 23 peter 164
4353 23 Aug 23 peter 165   double mSmithWaterman2::Aligner1::init(size_t r, size_t q) const
4353 23 Aug 23 peter 166   {
4353 23 Aug 23 peter 167     return 0;
4353 23 Aug 23 peter 168   }
4353 23 Aug 23 peter 169
4353 23 Aug 23 peter 170
4353 23 Aug 23 peter 171   void
4353 23 Aug 23 peter 172   mSmithWaterman2::Aligner1::update(const std::vector<Alignment>& alignments,
4353 23 Aug 23 peter 173                                     size_t r, size_t rsize)
4353 23 Aug 23 peter 174   {
4353 23 Aug 23 peter 175     col_best_.resize(alignments.size());
4353 23 Aug 23 peter 176     for (size_t i=0; i<col_best_.size(); ++i)
4353 23 Aug 23 peter 177       col_best_[i] = std::max(col_best_[i], alignments[i]);
4353 23 Aug 23 peter 178   }
4353 23 Aug 23 peter 179
4353 23 Aug 23 peter 180
4353 23 Aug 23 peter 181   // Aligner 2
4353 23 Aug 23 peter 182   mSmithWaterman2::Aligner2::Aligner2(const std::vector<Alignment>& alignments,
4353 23 Aug 23 peter 183                                       double insertion_penalty,
4353 23 Aug 23 peter 184                                       double open_insertion_penalty,
4353 23 Aug 23 peter 185                                       double deletion_penalty,
4353 23 Aug 23 peter 186                                       double open_deletion_penalty)
4353 23 Aug 23 peter 187     : mAligner<mSmithWaterman2::Aligner2>(Affine(insertion_penalty,
4353 23 Aug 23 peter 188                                                  open_insertion_penalty),
4353 23 Aug 23 peter 189                                           Affine(deletion_penalty,
4353 23 Aug 23 peter 190                                                  open_deletion_penalty)),
4353 23 Aug 23 peter 191       alignments_(alignments)
4353 23 Aug 23 peter 192   {}
4353 23 Aug 23 peter 193
4353 23 Aug 23 peter 194
4353 23 Aug 23 peter 195   double mSmithWaterman2::Aligner2::init_row(size_t q) const
4353 23 Aug 23 peter 196   {
4353 23 Aug 23 peter 197     assert(q < alignments_.size());
4353 23 Aug 23 peter 198     return alignments_[q].score();
4353 23 Aug 23 peter 199   }
4353 23 Aug 23 peter 200
4353 23 Aug 23 peter 201
4353 23 Aug 23 peter 202   double mSmithWaterman2::Aligner2::init_col(size_t r) const
4353 23 Aug 23 peter 203   {
4353 23 Aug 23 peter 204     assert(alignments_.size());
4353 23 Aug 23 peter 205     return alignments_[0].score();
4353 23 Aug 23 peter 206   }
4353 23 Aug 23 peter 207
4353 23 Aug 23 peter 208
4353 23 Aug 23 peter 209   double mSmithWaterman2::Aligner2::init(size_t r, size_t q) const
4353 23 Aug 23 peter 210   {
4353 23 Aug 23 peter 211     assert(1+q < alignments_.size());
4353 23 Aug 23 peter 212     return alignments_[1+q].score();
4353 23 Aug 23 peter 213   }
4353 23 Aug 23 peter 214
4353 23 Aug 23 peter 215
4353 23 Aug 23 peter 216   void
4353 23 Aug 23 peter 217   mSmithWaterman2::Aligner2::update(const std::vector<Alignment>& alignments,
4353 23 Aug 23 peter 218                                     size_t r, size_t rsize)
4353 23 Aug 23 peter 219   {
4353 23 Aug 23 peter 220     auto best = std::max_element(alignments.begin(), alignments.end());
4353 23 Aug 23 peter 221
4353 23 Aug 23 peter 222     if (best_ < *best) {
4353 23 Aug 23 peter 223       best_ = *best;
4353 23 Aug 23 peter 224       best_.cigar().push_back(BAM_CSOFT_CLIP, (alignments.end() - best) - 1);
4353 23 Aug 23 peter 225     }
4353 23 Aug 23 peter 226   }
4353 23 Aug 23 peter 227
4353 23 Aug 23 peter 228 }}}