yat/utility/mSmithWaterman2.h

Code
Comments
Other
Rev Date Author Line
4353 23 Aug 23 peter 1 #ifndef _theplu_yat_utility_msmith_waterman2
4353 23 Aug 23 peter 2 #define _theplu_yat_utility_msmith_waterman2
4353 23 Aug 23 peter 3 // $Id$
4353 23 Aug 23 peter 4
4353 23 Aug 23 peter 5 /*
4353 23 Aug 23 peter 6   Copyright (C) 2023 Peter Johansson
4353 23 Aug 23 peter 7
4353 23 Aug 23 peter 8   This file is part of the yat library, https://dev.thep.lu.se/yat
4353 23 Aug 23 peter 9
4353 23 Aug 23 peter 10   The yat library is free software; you can redistribute it and/or
4353 23 Aug 23 peter 11   modify it under the terms of the GNU General Public License as
4353 23 Aug 23 peter 12   published by the Free Software Foundation; either version 3 of the
4353 23 Aug 23 peter 13   License, or (at your option) any later version.
4353 23 Aug 23 peter 14
4353 23 Aug 23 peter 15   The yat library is distributed in the hope that it will be useful,
4353 23 Aug 23 peter 16   but WITHOUT ANY WARRANTY; without even the implied warranty of
4353 23 Aug 23 peter 17   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4353 23 Aug 23 peter 18   General Public License for more details.
4353 23 Aug 23 peter 19
4353 23 Aug 23 peter 20   You should have received a copy of the GNU General Public License
4353 23 Aug 23 peter 21   along with yat. If not, see <https://www.gnu.org/licenses/>.
4353 23 Aug 23 peter 22 */
4353 23 Aug 23 peter 23
4353 23 Aug 23 peter 24 #include "Aligner.h"
4353 23 Aug 23 peter 25 #include "mAligner.h"
4353 23 Aug 23 peter 26
4353 23 Aug 23 peter 27 #include "yat_assert.h"
4353 23 Aug 23 peter 28
4353 23 Aug 23 peter 29 namespace theplu {
4353 23 Aug 23 peter 30 namespace yat {
4353 23 Aug 23 peter 31 namespace utility {
4353 23 Aug 23 peter 32
4353 23 Aug 23 peter 33   /**
4353 23 Aug 23 peter 34      An aligner that aligns a query sequence against two reference
4353 23 Aug 23 peter 35      sequences such that the query is built up from the alignment
4353 23 Aug 23 peter 36      against the reference1, a potential insertion, and an alignment
4353 23 Aug 23 peter 37      against the reference2. It uses the Smith-Waterman algorithm
4353 23 Aug 23 peter 38      (SMA) to align the query against the two references and combines
4353 23 Aug 23 peter 39      this such that the SMA score against the first reference plus the
4353 23 Aug 23 peter 40      SMA score against the second reference is maximized, potentially
4353 23 Aug 23 peter 41      penalizing an insertion between the two alignments.
4353 23 Aug 23 peter 42
4353 23 Aug 23 peter 43      \since New in yat 0.21
4353 23 Aug 23 peter 44    */
4353 23 Aug 23 peter 45   class mSmithWaterman2
4353 23 Aug 23 peter 46   {
4353 23 Aug 23 peter 47   public:
4353 23 Aug 23 peter 48     /**
4353 23 Aug 23 peter 49        Equivalent with mSmithWaterman2(gap, open_gap, gap, open_gap)
4353 23 Aug 23 peter 50      */
4353 23 Aug 23 peter 51     mSmithWaterman2(double gap, double open_gap);
4353 23 Aug 23 peter 52
4353 23 Aug 23 peter 53     /**
4353 23 Aug 23 peter 54        Equivalent with
4353 23 Aug 23 peter 55        mSmithWaterman2(insertion_penalty, open_insertion_penalty,
4353 23 Aug 23 peter 56                        deletion_penalty, open_deletion_penalty,
4353 23 Aug 23 peter 57                        insertion_penalty, open_insertion_penalty,
4353 23 Aug 23 peter 58                        deletion_penalty, open_deletion_penalty)
4353 23 Aug 23 peter 59      */
4353 23 Aug 23 peter 60     mSmithWaterman2(double insertion_penalty, double open_insertion_penalty,
4353 23 Aug 23 peter 61                     double deletion_penalty, double open_deletion_penalty);
4353 23 Aug 23 peter 62
4353 23 Aug 23 peter 63     /**
4353 23 Aug 23 peter 64        Equivalent with
4353 23 Aug 23 peter 65        mSmithWaterman2(insertion_penalty1, open_insertion_penalty1,
4353 23 Aug 23 peter 66                        deletion_penalty1, open_deletion_penalty1,
4353 23 Aug 23 peter 67                        insertion_penalty2, open_insertion_penalty2,
4353 23 Aug 23 peter 68                        deletion_penalty2, open_deletion_penalty2,
4353 23 Aug 23 peter 69                        0, 0)
4353 23 Aug 23 peter 70      */
4353 23 Aug 23 peter 71     mSmithWaterman2(double insertion_penalty1, double open_insertion_penalty1,
4353 23 Aug 23 peter 72                     double deletion_penalty1, double open_deletion_penalty1,
4353 23 Aug 23 peter 73                     double insertion_penalty2, double open_insertion_penalty2,
4353 23 Aug 23 peter 74                     double deletion_penalty2, double open_deletion_penalty2);
4353 23 Aug 23 peter 75
4353 23 Aug 23 peter 76     /**
4353 23 Aug 23 peter 77        The alignment score against reference1 is calculated with a
4353 23 Aug 23 peter 78        penalty for k-sized insertions as
4353 23 Aug 23 peter 79        k*insertion_penalty1 + open_insertion_penalty1
4353 23 Aug 23 peter 80
4353 23 Aug 23 peter 81        penalty for k-sized deletions as
4353 23 Aug 23 peter 82        k*deletion_penalty1 + open_deletion_penalty1
4353 23 Aug 23 peter 83
4353 23 Aug 23 peter 84        The alignment score against reference2 is calculated with a
4353 23 Aug 23 peter 85        penalty for k-sized insertions as
4353 23 Aug 23 peter 86        k*insertion_penalty2 + open_insertion_penalty2
4353 23 Aug 23 peter 87
4353 23 Aug 23 peter 88        and penalty for k-sized deletions as
4353 23 Aug 23 peter 89        k*deletion_penalty2 + open_deletion_penalty2
4353 23 Aug 23 peter 90
4353 23 Aug 23 peter 91        The penalty for a k-sized insertion between the alignment
4353 23 Aug 23 peter 92        against reference1 and reference2 is calculated as
4353 23 Aug 23 peter 93        k*h=gap_penalty + open_gap_penalty
4353 23 Aug 23 peter 94      */
4353 23 Aug 23 peter 95     mSmithWaterman2(double insertion_penalty1, double open_insertion_penalty1,
4353 23 Aug 23 peter 96                     double deletion_penalty1, double open_deletion_penalty1,
4353 23 Aug 23 peter 97                     double insertion_penalty2, double open_insertion_penalty2,
4353 23 Aug 23 peter 98                     double deletion_penalty2, double open_deletion_penalty2,
4353 23 Aug 23 peter 99                     double gap_penalty, double open_gap_penalty);
4353 23 Aug 23 peter 100
4353 23 Aug 23 peter 101     /**
4353 23 Aug 23 peter 102        Aligned the query sequence in [\c query_begin, \c query_end)
4353 23 Aug 23 peter 103        against reference1, [reference1_begin, reference1_end), and
4353 23 Aug 23 peter 104        reference2, [reference2_begin, reference2_end). For matches
4353 23 Aug 23 peter 105        (or mismatches) in the alignment the alignment score is defined
4353 23 Aug 23 peter 106        by the kernel function, \c kernel.
4353 23 Aug 23 peter 107
4353 23 Aug 23 peter 108        Type Requirements:
4353 23 Aug 23 peter 109        - \c RandomAccessIterator1 is a \random_access_traversal_iterator
4353 23 Aug 23 peter 110        - \c RandomAccessIterator2 is a \random_access_traversal_iterator
4353 23 Aug 23 peter 111        - \c KernelFunction has an operator() that takes two arguments and
4353 23 Aug 23 peter 112        returns a type convertible to \c double.
4353 23 Aug 23 peter 113        - \c RandomAccessIterator1 is convertible to the first argument
4353 23 Aug 23 peter 114        type of KernelFunction::operator()
4353 23 Aug 23 peter 115        - \c RandomAccessIterator2 is convertible to the second argument
4353 23 Aug 23 peter 116        type of KernelFunction::operator()
4353 23 Aug 23 peter 117      */
4353 23 Aug 23 peter 118     template<typename RandomAccessIterator1, typename RandomAccessIterator2,
4353 23 Aug 23 peter 119              class KernelFunction>
4353 23 Aug 23 peter 120     double operator()(RandomAccessIterator1 reference1_begin,
4353 23 Aug 23 peter 121                       RandomAccessIterator1 reference1_end,
4353 23 Aug 23 peter 122                       RandomAccessIterator1 reference2_begin,
4353 23 Aug 23 peter 123                       RandomAccessIterator1 reference2_end,
4353 23 Aug 23 peter 124                       RandomAccessIterator2 query_begin,
4353 23 Aug 23 peter 125                       RandomAccessIterator2 query_end,
4353 23 Aug 23 peter 126                       const KernelFunction& kernel);
4353 23 Aug 23 peter 127
4353 23 Aug 23 peter 128     /**
4353 23 Aug 23 peter 129        Class that contains the alignment result against one of the
4353 23 Aug 23 peter 130        reference sequences.
4353 23 Aug 23 peter 131      */
4353 23 Aug 23 peter 132     class Result
4353 23 Aug 23 peter 133     {
4353 23 Aug 23 peter 134     public:
4353 23 Aug 23 peter 135       /**
4353 23 Aug 23 peter 136          If the position is p, the reference[p] is the first element
4353 23 Aug 23 peter 137          in the reference included in the alignment, that is, if
4353 23 Aug 23 peter 138          position is zero, the first base in reference is included in
4353 23 Aug 23 peter 139          the alignment.
4353 23 Aug 23 peter 140        */
4353 23 Aug 23 peter 141       size_t position(void) const;
4353 23 Aug 23 peter 142
4353 23 Aug 23 peter 143       /**
4353 23 Aug 23 peter 144          A CIGAR class describing the alignment. The alignment uses
4353 23 Aug 23 peter 145          BAM_CMATCH, 'm'.
4353 23 Aug 23 peter 146
4353 23 Aug 23 peter 147          \see mSmithWaterman::cigar()
4353 23 Aug 23 peter 148        */
4353 23 Aug 23 peter 149       const Aligner::Cigar& cigar(void) const;
4353 23 Aug 23 peter 150
4353 23 Aug 23 peter 151       /**
4353 23 Aug 23 peter 152          Same as cigar(), but uses BAM_CEQUAL ('=') and BAM_CDIFF ('X')
4353 23 Aug 23 peter 153          rather than BAM_CMATCH.
4353 23 Aug 23 peter 154
4353 23 Aug 23 peter 155          \see mSmithWaterman::cigar2()
4353 23 Aug 23 peter 156        */
4353 23 Aug 23 peter 157       const Aligner::Cigar& cigar2(void) const;
4353 23 Aug 23 peter 158
4353 23 Aug 23 peter 159       /**
4353 23 Aug 23 peter 160          \brief alignment score
4353 23 Aug 23 peter 161       */
4353 23 Aug 23 peter 162       double score(void) const;
4353 23 Aug 23 peter 163     private:
4353 23 Aug 23 peter 164       friend class mSmithWaterman2;
4353 23 Aug 23 peter 165       size_t position_;
4353 23 Aug 23 peter 166       Aligner::Cigar cigar_;
4353 23 Aug 23 peter 167       Aligner::Cigar cigar2_;
4353 23 Aug 23 peter 168       double score_;
4353 23 Aug 23 peter 169     };
4353 23 Aug 23 peter 170
4353 23 Aug 23 peter 171     /**
4353 23 Aug 23 peter 172        Object describing the alignment against the first reference.
4353 23 Aug 23 peter 173      */
4353 23 Aug 23 peter 174     const Result& first(void) const;
4353 23 Aug 23 peter 175
4353 23 Aug 23 peter 176     /**
4353 23 Aug 23 peter 177        Object describing the alignment against the second reference.
4353 23 Aug 23 peter 178      */
4353 23 Aug 23 peter 179     const Result& second(void) const;
4353 23 Aug 23 peter 180
4353 23 Aug 23 peter 181     /**
4353 23 Aug 23 peter 182        Size of the insertion in the query sequence between the two
4353 23 Aug 23 peter 183        alignments.
4353 23 Aug 23 peter 184      */
4353 23 Aug 23 peter 185     size_t gap(void) const;
4353 23 Aug 23 peter 186
4353 23 Aug 23 peter 187     /**
4353 23 Aug 23 peter 188        \return the total score, i.e., the sum of the alignment against
4353 23 Aug 23 peter 189        first reference, penalty for intermediate insertion, and
4353 23 Aug 23 peter 190        alignment against second reference.
4353 23 Aug 23 peter 191      */
4353 23 Aug 23 peter 192     double score(void) const;
4353 23 Aug 23 peter 193   private:
4353 23 Aug 23 peter 194     double insertion_penalty1_;
4353 23 Aug 23 peter 195     double open_insertion_penalty1_;
4353 23 Aug 23 peter 196     double deletion_penalty1_;
4353 23 Aug 23 peter 197     double open_deletion_penalty1_;
4353 23 Aug 23 peter 198     double insertion_penalty2_;
4353 23 Aug 23 peter 199     double open_insertion_penalty2_;
4353 23 Aug 23 peter 200     double deletion_penalty2_;
4353 23 Aug 23 peter 201     double open_deletion_penalty2_;
4353 23 Aug 23 peter 202     double gap_penalty_;
4353 23 Aug 23 peter 203     double open_gap_penalty_;
4353 23 Aug 23 peter 204
4353 23 Aug 23 peter 205     /**
4353 23 Aug 23 peter 206        Class used to align against second reference
4353 23 Aug 23 peter 207      */
4353 23 Aug 23 peter 208     class Aligner1 : public mAligner<mSmithWaterman2::Aligner1>
4353 23 Aug 23 peter 209     {
4353 23 Aug 23 peter 210     public:
4353 23 Aug 23 peter 211       Aligner1(double insertion_penalty,
4353 23 Aug 23 peter 212                double open_insertion_penalty,
4353 23 Aug 23 peter 213                double deletion_penalty,
4353 23 Aug 23 peter 214                double open_deletion_penalty);
4353 23 Aug 23 peter 215
4353 23 Aug 23 peter 216       const std::vector<Alignment>& col_best(void) const;
4353 23 Aug 23 peter 217     private:
4353 23 Aug 23 peter 218       std::vector<Alignment> col_best_;
4353 23 Aug 23 peter 219
4353 23 Aug 23 peter 220       friend class mAligner;
4353 23 Aug 23 peter 221       double init_row(size_t q) const;
4353 23 Aug 23 peter 222       double init_col(size_t r) const;
4353 23 Aug 23 peter 223       double init(size_t r, size_t q) const;
4353 23 Aug 23 peter 224       void update(const std::vector<Alignment>& alignments, size_t r,
4353 23 Aug 23 peter 225                   size_t rsize);
4353 23 Aug 23 peter 226     };
4353 23 Aug 23 peter 227
4353 23 Aug 23 peter 228
4353 23 Aug 23 peter 229     /**
4353 23 Aug 23 peter 230        Class used to align against second reference
4353 23 Aug 23 peter 231      */
4353 23 Aug 23 peter 232     class Aligner2 : public mAligner<mSmithWaterman2::Aligner2>
4353 23 Aug 23 peter 233     {
4353 23 Aug 23 peter 234     public:
4353 23 Aug 23 peter 235       Aligner2(const std::vector<Alignment>& alignments,
4353 23 Aug 23 peter 236                double insertion_penalty1,
4353 23 Aug 23 peter 237                double open_insertion_penalty1,
4353 23 Aug 23 peter 238                double deletion_penalty1,
4353 23 Aug 23 peter 239                double open_deletion_penalty1);
4353 23 Aug 23 peter 240
4353 23 Aug 23 peter 241     private:
4353 23 Aug 23 peter 242       const std::vector<Alignment>& alignments_;
4353 23 Aug 23 peter 243       friend class mAligner;
4353 23 Aug 23 peter 244       double init_row(size_t q) const;
4353 23 Aug 23 peter 245       double init_col(size_t r) const;
4353 23 Aug 23 peter 246       double init(size_t r, size_t q) const;
4353 23 Aug 23 peter 247       void update(const std::vector<Alignment>& alignments, size_t r,
4353 23 Aug 23 peter 248                   size_t rsize);
4353 23 Aug 23 peter 249     };
4353 23 Aug 23 peter 250
4353 23 Aug 23 peter 251     Result first_;
4353 23 Aug 23 peter 252     Result second_;
4353 23 Aug 23 peter 253     size_t gap_;
4353 23 Aug 23 peter 254     double score_;
4353 23 Aug 23 peter 255   };
4353 23 Aug 23 peter 256
4353 23 Aug 23 peter 257
4353 23 Aug 23 peter 258   // template implementation
4353 23 Aug 23 peter 259
4353 23 Aug 23 peter 260   template<typename RandomAccessIterator1, typename RandomAccessIterator2,
4353 23 Aug 23 peter 261            class KernelFunction>
4353 23 Aug 23 peter 262   double mSmithWaterman2::operator()(RandomAccessIterator1 reference1_begin,
4353 23 Aug 23 peter 263                                      RandomAccessIterator1 reference1_end,
4353 23 Aug 23 peter 264                                      RandomAccessIterator1 reference2_begin,
4353 23 Aug 23 peter 265                                      RandomAccessIterator1 reference2_end,
4353 23 Aug 23 peter 266                                      RandomAccessIterator2 query_begin,
4353 23 Aug 23 peter 267                                      RandomAccessIterator2 query_end,
4353 23 Aug 23 peter 268                                      const KernelFunction& kernel)
4353 23 Aug 23 peter 269     {
4353 23 Aug 23 peter 270       mSmithWaterman2::Aligner1
4353 23 Aug 23 peter 271         aligner1(insertion_penalty1_, open_insertion_penalty1_,
4353 23 Aug 23 peter 272                  deletion_penalty1_, open_deletion_penalty1_);
4353 23 Aug 23 peter 273
4353 23 Aug 23 peter 274       // align against reference1
4353 23 Aug 23 peter 275       aligner1(reference1_begin, reference1_end,
4353 23 Aug 23 peter 276                query_begin, query_end, kernel);
4353 23 Aug 23 peter 277
4353 23 Aug 23 peter 278       // allow inserted bases between the two alignments
4353 23 Aug 23 peter 279       std::vector<mAlignerBase::Alignment>
4353 23 Aug 23 peter 280         alignments(aligner1.col_best().size());
4353 23 Aug 23 peter 281       alignments[0].score() = aligner1.col_best()[0].score();
4353 23 Aug 23 peter 282       for (size_t i=1; i<alignments.size(); ++i) {
4353 23 Aug 23 peter 283         alignments[i].score() = aligner1.col_best()[i].score();
4353 23 Aug 23 peter 284
4353 23 Aug 23 peter 285         double cand_score = alignments[i-1].score() - gap_penalty_;
4353 23 Aug 23 peter 286         if (!alignments[i-1].cigar().size())
4353 23 Aug 23 peter 287           cand_score -= open_gap_penalty_;
4353 23 Aug 23 peter 288
4353 23 Aug 23 peter 289         if (cand_score > alignments[i].score()) {
4353 23 Aug 23 peter 290           alignments[i].score() = cand_score;
4353 23 Aug 23 peter 291           alignments[i].cigar() = alignments[i-1].cigar();
4353 23 Aug 23 peter 292           alignments[i].cigar().push_back(BAM_CINS, 1);
4353 23 Aug 23 peter 293         }
4353 23 Aug 23 peter 294       }
4353 23 Aug 23 peter 295
4353 23 Aug 23 peter 296       mSmithWaterman2::Aligner2
4353 23 Aug 23 peter 297         aligner2(alignments,
4353 23 Aug 23 peter 298                  insertion_penalty2_, open_insertion_penalty2_,
4353 23 Aug 23 peter 299                  deletion_penalty2_, open_deletion_penalty2_);
4353 23 Aug 23 peter 300
4353 23 Aug 23 peter 301       aligner2(reference2_begin, reference2_end,
4353 23 Aug 23 peter 302                query_begin, query_end, kernel);
4353 23 Aug 23 peter 303
4353 23 Aug 23 peter 304       size_t q2 = 0;
4353 23 Aug 23 peter 305       if (aligner2.cigar().size() && aligner2.cigar().op(0) == BAM_CSOFT_CLIP)
4353 23 Aug 23 peter 306         q2 = aligner2.cigar().oplen(0);
4353 23 Aug 23 peter 307
4353 23 Aug 23 peter 308       size_t qsize = std::distance(query_begin, query_end);
4353 23 Aug 23 peter 309
4353 23 Aug 23 peter 310       score_ = aligner2.score();
4353 23 Aug 23 peter 311       second_.score_ = score_ - alignments[q2].score();
4353 23 Aug 23 peter 312       second_.position_ = aligner2.position();
4353 23 Aug 23 peter 313       second_.cigar_ = aligner2.cigar();
4353 23 Aug 23 peter 314       second_.cigar2_ = aligner2.cigar2();
4353 23 Aug 23 peter 315
4353 23 Aug 23 peter 316       gap_ = 0;
4353 23 Aug 23 peter 317       size_t cig_size = alignments[q2].cigar().size();
4353 23 Aug 23 peter 318       if (cig_size)
4353 23 Aug 23 peter 319         gap_ = alignments[q2].cigar().oplen(cig_size-1);
4353 23 Aug 23 peter 320
4353 23 Aug 23 peter 321       size_t q1 = q2 - gap_;
4353 23 Aug 23 peter 322       first_.score_ = aligner1.col_best()[q1].score();
4353 23 Aug 23 peter 323       first_.position_ = aligner1.col_best()[q1].position();
4353 23 Aug 23 peter 324       first_.cigar_ = aligner1.col_best()[q1].cigar();
4353 23 Aug 23 peter 325       first_.cigar_.push_back(BAM_CSOFT_CLIP,
4353 23 Aug 23 peter 326                               qsize - first_.cigar_.query_length());
4353 23 Aug 23 peter 327       first_.cigar2_ =
4353 23 Aug 23 peter 328         cigar2(reference1_begin+first_.position_, reference1_end,
4353 23 Aug 23 peter 329                query_begin, query_end, kernel, first_.cigar_);
4353 23 Aug 23 peter 330
4353 23 Aug 23 peter 331       return score();
4353 23 Aug 23 peter 332     }
4353 23 Aug 23 peter 333
4353 23 Aug 23 peter 334
4353 23 Aug 23 peter 335 }}}
4353 23 Aug 23 peter 336 #endif

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