yat/utility/Aligner.cc

Code
Comments
Other
Rev Date Author Line
2815 28 Aug 12 peter 1 // $Id$
2815 28 Aug 12 peter 2
2815 28 Aug 12 peter 3 /*
4308 10 Feb 23 peter 4   Copyright (C) 2012, 2014, 2016, 2019, 2020, 2022, 2023 Peter Johansson
2815 28 Aug 12 peter 5
2815 28 Aug 12 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
2815 28 Aug 12 peter 7
2815 28 Aug 12 peter 8   The yat library is free software; you can redistribute it and/or
2815 28 Aug 12 peter 9   modify it under the terms of the GNU General Public License as
2815 28 Aug 12 peter 10   published by the Free Software Foundation; either version 3 of the
2815 28 Aug 12 peter 11   License, or (at your option) any later version.
2815 28 Aug 12 peter 12
2815 28 Aug 12 peter 13   The yat library is distributed in the hope that it will be useful,
2815 28 Aug 12 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
2815 28 Aug 12 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2815 28 Aug 12 peter 16   General Public License for more details.
2815 28 Aug 12 peter 17
2815 28 Aug 12 peter 18   You should have received a copy of the GNU General Public License
2815 28 Aug 12 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
2815 28 Aug 12 peter 20 */
2815 28 Aug 12 peter 21
2881 18 Nov 12 peter 22 #include <config.h>
2881 18 Nov 12 peter 23
2815 28 Aug 12 peter 24 #include "Aligner.h"
3212 05 May 14 peter 25 #include "Cigar.h"
2815 28 Aug 12 peter 26 #include "Matrix.h"
2815 28 Aug 12 peter 27
3462 27 Jan 16 peter 28 #include <algorithm>
2815 28 Aug 12 peter 29 #include <cassert>
4003 15 Oct 20 peter 30 #include <cctype>
2815 28 Aug 12 peter 31 #include <limits>
3295 25 Jul 14 peter 32 #include <ostream>
2815 28 Aug 12 peter 33
2815 28 Aug 12 peter 34 namespace theplu {
2815 28 Aug 12 peter 35 namespace yat {
2815 28 Aug 12 peter 36 namespace utility {
2815 28 Aug 12 peter 37
2815 28 Aug 12 peter 38   Aligner::Aligner(double gap, double open_gap)
2815 28 Aug 12 peter 39     : vertical_gap_(gap), open_vertical_gap_(open_gap), horizon_gap_(gap),
2815 28 Aug 12 peter 40       open_horizon_gap_(open_gap)
2815 28 Aug 12 peter 41   {
2815 28 Aug 12 peter 42   }
2815 28 Aug 12 peter 43
2815 28 Aug 12 peter 44
2815 28 Aug 12 peter 45   Aligner::Aligner(double vertical_gap, double open_vertical_gap,
2815 28 Aug 12 peter 46                    double horizon_gap, double open_horizon_gap)
2815 28 Aug 12 peter 47     : vertical_gap_(vertical_gap), open_vertical_gap_(open_vertical_gap),
2815 28 Aug 12 peter 48       horizon_gap_(horizon_gap), open_horizon_gap_(open_horizon_gap)
2815 28 Aug 12 peter 49   {
2815 28 Aug 12 peter 50   }
2815 28 Aug 12 peter 51
2815 28 Aug 12 peter 52
4125 14 Jan 22 peter 53   double Aligner::needleman_wunsch(const MatrixBase& s)
2815 28 Aug 12 peter 54   {
2853 23 Sep 12 peter 55     Matrix m(s.rows()+1,s.columns()+1, -std::numeric_limits<double>::max());
2853 23 Sep 12 peter 56     m(0, 0) = 0.0;
2815 28 Aug 12 peter 57     // Init upper and left border of matrix
2815 28 Aug 12 peter 58     for (size_t i=1; i<m.rows(); i++)
4308 10 Feb 23 peter 59       m(i,0) = - (i * vertical_gap_ + open_vertical_gap_);
2815 28 Aug 12 peter 60     for (size_t i=1; i<m.columns(); i++)
4308 10 Feb 23 peter 61       m(0,i) = - (i * horizon_gap_ + open_horizon_gap_);
2815 28 Aug 12 peter 62     (*this)(s, m);
2815 28 Aug 12 peter 63     // return lower right corner element
2815 28 Aug 12 peter 64     return m(s.rows(), s.columns());
2815 28 Aug 12 peter 65   }
2815 28 Aug 12 peter 66
2815 28 Aug 12 peter 67
4125 14 Jan 22 peter 68   double Aligner::smith_waterman(const MatrixBase& s)
2815 28 Aug 12 peter 69   {
2815 28 Aug 12 peter 70     Matrix m(s.rows()+1,s.columns()+1);
2815 28 Aug 12 peter 71     (*this)(s, m);
2815 28 Aug 12 peter 72     return max(m);
2815 28 Aug 12 peter 73   }
2815 28 Aug 12 peter 74
2815 28 Aug 12 peter 75
4125 14 Jan 22 peter 76   void Aligner::operator()(const MatrixBase& score, MatrixMutable& x)
2815 28 Aug 12 peter 77   {
2815 28 Aug 12 peter 78     assert(x.columns() == score.columns()+1);
2815 28 Aug 12 peter 79     assert(x.rows() == score.rows()+1);
2815 28 Aug 12 peter 80
2815 28 Aug 12 peter 81     columns_ = x.columns();
3342 06 Nov 14 peter 82     alignment_.assign(columns_*x.rows(), none);
2815 28 Aug 12 peter 83
2815 28 Aug 12 peter 84     for (size_t i=1; i<x.rows(); ++i)
2815 28 Aug 12 peter 85       for (size_t j=1; j<x.columns(); ++j) {
2815 28 Aug 12 peter 86         // match (or mismatch)
2815 28 Aug 12 peter 87         if (x(i-1, j-1) + score(i-1, j-1) > x(i,j)) {
2815 28 Aug 12 peter 88           x(i,j) = x(i-1, j-1) + score(i-1, j-1);
3194 22 Apr 14 peter 89           directions(i, j) = diagonal;
2815 28 Aug 12 peter 90         }
2815 28 Aug 12 peter 91
3194 22 Apr 14 peter 92         if (directions(i, j-1) == right) {
2853 23 Sep 12 peter 93           if (x(i, j-1)-horizon_gap_ > x(i,j)) {
2815 28 Aug 12 peter 94             x(i,j) = x(i,j-1) - horizon_gap_;
3194 22 Apr 14 peter 95             directions(i,j) = right;
2815 28 Aug 12 peter 96           }
2815 28 Aug 12 peter 97         }
2815 28 Aug 12 peter 98         else if(x(i,j-1)-horizon_gap_-open_horizon_gap_ > x(i,j)) {
2815 28 Aug 12 peter 99           x(i,j) = x(i,j-1) - horizon_gap_ - open_horizon_gap_;
3194 22 Apr 14 peter 100           directions(i,j) = right;
2815 28 Aug 12 peter 101         }
2815 28 Aug 12 peter 102
3194 22 Apr 14 peter 103         if (directions(i-1,j) == down) {
2815 28 Aug 12 peter 104           if (x(i-1,j) - vertical_gap_ > x(i,j)) {
2815 28 Aug 12 peter 105             x(i,j) = x(i-1,j) - vertical_gap_;
3194 22 Apr 14 peter 106             directions(i,j) = down;
2815 28 Aug 12 peter 107           }
2815 28 Aug 12 peter 108         }
2815 28 Aug 12 peter 109         else if (x(i-1,j) - vertical_gap_ - open_vertical_gap_ > x(i,j)) {
2815 28 Aug 12 peter 110           x(i,j) = x(i-1,j) - vertical_gap_ - open_vertical_gap_;
3194 22 Apr 14 peter 111           directions(i,j) = down;
2815 28 Aug 12 peter 112         }
2815 28 Aug 12 peter 113       }
2815 28 Aug 12 peter 114   }
2815 28 Aug 12 peter 115
2815 28 Aug 12 peter 116
2816 28 Aug 12 peter 117   const Aligner::direction& Aligner::alignment(size_t i, size_t j) const
2815 28 Aug 12 peter 118   {
2816 28 Aug 12 peter 119     return alignment_[i*columns_ + j];
2815 28 Aug 12 peter 120   }
2815 28 Aug 12 peter 121
2815 28 Aug 12 peter 122
3200 03 May 14 peter 123   const Aligner::Cigar Aligner::cigar(size_t row, size_t column) const
3200 03 May 14 peter 124   {
3200 03 May 14 peter 125     Cigar cig;
3200 03 May 14 peter 126     while (row || column) {
3200 03 May 14 peter 127       switch (alignment(row, column)) {
3200 03 May 14 peter 128       case(diagonal):
3200 03 May 14 peter 129         assert(row);
3200 03 May 14 peter 130         --row;
3200 03 May 14 peter 131         assert(column);
3200 03 May 14 peter 132         --column;
3200 03 May 14 peter 133         cig.push_front(BAM_CMATCH);
3200 03 May 14 peter 134         break;
3200 03 May 14 peter 135       case (right):
3200 03 May 14 peter 136         assert(column);
3200 03 May 14 peter 137         --column;
3200 03 May 14 peter 138         cig.push_front(BAM_CINS);
3200 03 May 14 peter 139         break;
3200 03 May 14 peter 140       case (down):
3200 03 May 14 peter 141         assert(row);
3200 03 May 14 peter 142         --row;
3200 03 May 14 peter 143         cig.push_front(BAM_CDEL);
3200 03 May 14 peter 144         break;
3200 03 May 14 peter 145       default:
3200 03 May 14 peter 146         return cig;
3200 03 May 14 peter 147       }
3200 03 May 14 peter 148     }
3200 03 May 14 peter 149     return cig;
3200 03 May 14 peter 150   }
3200 03 May 14 peter 151
3200 03 May 14 peter 152
3194 22 Apr 14 peter 153   Aligner::direction& Aligner::directions(size_t i, size_t j)
2815 28 Aug 12 peter 154   {
2816 28 Aug 12 peter 155     return alignment_[i*columns_ + j];
2815 28 Aug 12 peter 156   }
2815 28 Aug 12 peter 157
3200 03 May 14 peter 158
4334 25 Mar 23 peter 159   double Aligner::insertion_penalty(void) const
4334 25 Mar 23 peter 160   {
4334 25 Mar 23 peter 161     return vertical_gap_;
4334 25 Mar 23 peter 162   }
4334 25 Mar 23 peter 163
4334 25 Mar 23 peter 164
4334 25 Mar 23 peter 165   double Aligner::open_insertion_penalty(void) const
4334 25 Mar 23 peter 166   {
4334 25 Mar 23 peter 167     return open_vertical_gap_;
4334 25 Mar 23 peter 168   }
4334 25 Mar 23 peter 169
4334 25 Mar 23 peter 170
4334 25 Mar 23 peter 171   double Aligner::deletion_penalty(void) const
4334 25 Mar 23 peter 172   {
4334 25 Mar 23 peter 173     return horizon_gap_;
4334 25 Mar 23 peter 174   }
4334 25 Mar 23 peter 175
4334 25 Mar 23 peter 176
4334 25 Mar 23 peter 177   double Aligner::open_deletion_penalty(void) const
4334 25 Mar 23 peter 178   {
4334 25 Mar 23 peter 179     return open_horizon_gap_;
4334 25 Mar 23 peter 180   }
4334 25 Mar 23 peter 181
4334 25 Mar 23 peter 182
3200 03 May 14 peter 183   // ================ Cigar =======================
3200 03 May 14 peter 184
4003 15 Oct 20 peter 185   Aligner::Cigar::Cigar(const std::string& s)
4003 15 Oct 20 peter 186   {
4003 15 Oct 20 peter 187     try {
4003 15 Oct 20 peter 188       std::istringstream is(s);
4003 15 Oct 20 peter 189       while (parse_element(is));
4003 15 Oct 20 peter 190     }
4003 15 Oct 20 peter 191     catch (invalid_argument& e) {
4003 15 Oct 20 peter 192       std::ostringstream msg;
4003 15 Oct 20 peter 193       msg << "Aligner::Cigar: " << s;
4003 15 Oct 20 peter 194       std::throw_with_nested(invalid_argument(msg.str()));
4003 15 Oct 20 peter 195     }
4003 15 Oct 20 peter 196   }
4003 15 Oct 20 peter 197
4003 15 Oct 20 peter 198
3840 13 Aug 19 peter 199   Aligner::Cigar::const_iterator Aligner::Cigar::begin(void) const
3840 13 Aug 19 peter 200   {
3840 13 Aug 19 peter 201     return const_iterator(cigar_.begin());
3840 13 Aug 19 peter 202   }
3840 13 Aug 19 peter 203
3840 13 Aug 19 peter 204
3200 03 May 14 peter 205   void Aligner::Cigar::clear(void)
3200 03 May 14 peter 206   {
3200 03 May 14 peter 207     cigar_.clear();
3200 03 May 14 peter 208   }
3200 03 May 14 peter 209
3200 03 May 14 peter 210
3840 13 Aug 19 peter 211   Aligner::Cigar::const_iterator Aligner::Cigar::end(void) const
3840 13 Aug 19 peter 212   {
3840 13 Aug 19 peter 213     return const_iterator(cigar_.end());
3840 13 Aug 19 peter 214   }
3840 13 Aug 19 peter 215
3840 13 Aug 19 peter 216
3200 03 May 14 peter 217   uint8_t Aligner::Cigar::op(size_t i) const
3200 03 May 14 peter 218   {
3200 03 May 14 peter 219     assert(i<cigar_.size());
3200 03 May 14 peter 220     return bam_cigar_op(cigar_[i]);
3200 03 May 14 peter 221   }
3200 03 May 14 peter 222
3200 03 May 14 peter 223
3213 05 May 14 peter 224   uint32_t Aligner::Cigar::length(uint8_t mask) const
3213 05 May 14 peter 225   {
3213 05 May 14 peter 226     uint32_t res = 0;
3213 05 May 14 peter 227     for (size_t i=0; i<cigar_.size(); ++i)
3213 05 May 14 peter 228       if (bam_cigar_type(op(i)) & mask)
3213 05 May 14 peter 229         res += oplen(i);
3213 05 May 14 peter 230     return res;
3213 05 May 14 peter 231   }
3213 05 May 14 peter 232
3213 05 May 14 peter 233
3213 05 May 14 peter 234   uint32_t Aligner::Cigar::length(void) const
3213 05 May 14 peter 235   {
3213 05 May 14 peter 236     return length(3);
3213 05 May 14 peter 237   }
3213 05 May 14 peter 238
3213 05 May 14 peter 239
3200 03 May 14 peter 240   char Aligner::Cigar::opchr(size_t i) const
3200 03 May 14 peter 241   {
3200 03 May 14 peter 242     assert(i<cigar_.size());
3200 03 May 14 peter 243     return bam_cigar_opchr(cigar_[i]);
3200 03 May 14 peter 244   }
3200 03 May 14 peter 245
3200 03 May 14 peter 246
3200 03 May 14 peter 247   uint32_t Aligner::Cigar::oplen(size_t i) const
3200 03 May 14 peter 248   {
3200 03 May 14 peter 249     assert(i<cigar_.size());
3200 03 May 14 peter 250     return bam_cigar_oplen(cigar_[i]);
3200 03 May 14 peter 251   }
3200 03 May 14 peter 252
3200 03 May 14 peter 253
4003 15 Oct 20 peter 254   bool Aligner::Cigar::parse_element(std::istream& is)
4003 15 Oct 20 peter 255   {
4003 15 Oct 20 peter 256     char c = is.get();
4003 15 Oct 20 peter 257     if (!is.good())
4003 15 Oct 20 peter 258       return false;
4003 15 Oct 20 peter 259     if (!isdigit(c)) {
4003 15 Oct 20 peter 260       std::string msg("not a digit: '");
4003 15 Oct 20 peter 261       msg += c;
4003 15 Oct 20 peter 262       msg += "'";
4003 15 Oct 20 peter 263       throw invalid_argument(msg);
4003 15 Oct 20 peter 264     }
4003 15 Oct 20 peter 265     uint32_t len = c - '0';
4003 15 Oct 20 peter 266
4003 15 Oct 20 peter 267     while (is.get(c).good()) {
4003 15 Oct 20 peter 268       if (isdigit(c))
4003 15 Oct 20 peter 269         len = 10 * len + (c - '0');
4003 15 Oct 20 peter 270       else {
4003 15 Oct 20 peter 271         push_back(cigar_table(static_cast<uint8_t>(c)), len);
4003 15 Oct 20 peter 272         return true;
4003 15 Oct 20 peter 273       }
4003 15 Oct 20 peter 274     }
4003 15 Oct 20 peter 275     throw invalid_argument("Cigar::parse_element: missing operator");
4003 15 Oct 20 peter 276     return true;
4003 15 Oct 20 peter 277   }
4003 15 Oct 20 peter 278
4003 15 Oct 20 peter 279
3200 03 May 14 peter 280   void Aligner::Cigar::pop_back(uint32_t len)
3200 03 May 14 peter 281   {
3200 03 May 14 peter 282     while (len) {
3200 03 May 14 peter 283       assert(size());
3200 03 May 14 peter 284       size_t i = size()-1;
3200 03 May 14 peter 285       if (len < oplen(i)) {
3200 03 May 14 peter 286         cigar_[i] -= (len<<BAM_CIGAR_SHIFT);
3200 03 May 14 peter 287         return;
3200 03 May 14 peter 288       }
3200 03 May 14 peter 289       else {
3200 03 May 14 peter 290         len -= oplen(i);
3200 03 May 14 peter 291         cigar_.pop_back();
3200 03 May 14 peter 292       }
3200 03 May 14 peter 293     }
3200 03 May 14 peter 294   }
3200 03 May 14 peter 295
3200 03 May 14 peter 296
3200 03 May 14 peter 297   void Aligner::Cigar::pop_front(uint32_t len)
3200 03 May 14 peter 298   {
3200 03 May 14 peter 299     while (len) {
3200 03 May 14 peter 300       assert(size());
3200 03 May 14 peter 301       if (len < oplen(0)) {
3200 03 May 14 peter 302         cigar_[0] -= (len<<BAM_CIGAR_SHIFT);
3200 03 May 14 peter 303         return;
3200 03 May 14 peter 304       }
3200 03 May 14 peter 305       else {
3200 03 May 14 peter 306         len -= oplen(0);
3200 03 May 14 peter 307         cigar_.pop_front();
3200 03 May 14 peter 308       }
3200 03 May 14 peter 309     }
3200 03 May 14 peter 310   }
3200 03 May 14 peter 311
3200 03 May 14 peter 312
3200 03 May 14 peter 313   void Aligner::Cigar::push_back(uint8_t op, uint32_t len)
3200 03 May 14 peter 314   {
3205 04 May 14 peter 315     if (len==0)
3205 04 May 14 peter 316       return;
3200 03 May 14 peter 317     if (cigar_.empty() || this->op(size()-1)!=op)
3201 03 May 14 peter 318       cigar_.push_back(bam_cigar_gen(len, op));
3200 03 May 14 peter 319     else
3200 03 May 14 peter 320       cigar_[size()-1] += (len<<BAM_CIGAR_SHIFT);
3200 03 May 14 peter 321   }
3200 03 May 14 peter 322
3200 03 May 14 peter 323
3200 03 May 14 peter 324   void Aligner::Cigar::push_front(uint8_t op, uint32_t len)
3200 03 May 14 peter 325   {
3205 04 May 14 peter 326     if (len==0)
3205 04 May 14 peter 327       return;
3200 03 May 14 peter 328     if (cigar_.empty() || this->op(0)!=op)
3200 03 May 14 peter 329       cigar_.push_front(bam_cigar_gen(len, op));
3200 03 May 14 peter 330     else
3200 03 May 14 peter 331       cigar_[0] += (len<<BAM_CIGAR_SHIFT);
3200 03 May 14 peter 332   }
3200 03 May 14 peter 333
3200 03 May 14 peter 334
3213 05 May 14 peter 335   uint32_t Aligner::Cigar::query_length(void) const
3213 05 May 14 peter 336   {
3213 05 May 14 peter 337     return length(1);
3213 05 May 14 peter 338   }
3213 05 May 14 peter 339
3213 05 May 14 peter 340
4334 25 Mar 23 peter 341   Aligner::Cigar::const_reverse_iterator Aligner::Cigar::rbegin(void) const
4334 25 Mar 23 peter 342   {
4334 25 Mar 23 peter 343     return const_reverse_iterator(cigar_.rbegin());
4334 25 Mar 23 peter 344   }
4334 25 Mar 23 peter 345
4334 25 Mar 23 peter 346
3213 05 May 14 peter 347   uint32_t Aligner::Cigar::reference_length(void) const
3213 05 May 14 peter 348   {
3213 05 May 14 peter 349     return length(2);
3213 05 May 14 peter 350   }
3213 05 May 14 peter 351
3213 05 May 14 peter 352
4334 25 Mar 23 peter 353   Aligner::Cigar::const_reverse_iterator Aligner::Cigar::rend(void) const
4334 25 Mar 23 peter 354   {
4334 25 Mar 23 peter 355     return const_reverse_iterator(cigar_.rend());
4334 25 Mar 23 peter 356   }
4334 25 Mar 23 peter 357
4334 25 Mar 23 peter 358
3462 27 Jan 16 peter 359   void Aligner::Cigar::reverse(void)
3462 27 Jan 16 peter 360   {
3462 27 Jan 16 peter 361     std::reverse(cigar_.begin(), cigar_.end());
3462 27 Jan 16 peter 362   }
3462 27 Jan 16 peter 363
3462 27 Jan 16 peter 364
3200 03 May 14 peter 365   size_t Aligner::Cigar::size(void) const
3200 03 May 14 peter 366   {
3200 03 May 14 peter 367     return cigar_.size();
3200 03 May 14 peter 368   }
3200 03 May 14 peter 369
3200 03 May 14 peter 370
3200 03 May 14 peter 371   uint32_t Aligner::Cigar::operator[](size_t i) const
3200 03 May 14 peter 372   {
3200 03 May 14 peter 373     assert(i<cigar_.size());
3200 03 May 14 peter 374     return cigar_[i];
3200 03 May 14 peter 375   }
3200 03 May 14 peter 376
3200 03 May 14 peter 377
4334 25 Mar 23 peter 378   bool Aligner::Cigar::operator==(const Aligner::Cigar& other) const
4334 25 Mar 23 peter 379   {
4334 25 Mar 23 peter 380     return cigar_ == other.cigar_;
4334 25 Mar 23 peter 381   }
4334 25 Mar 23 peter 382
4334 25 Mar 23 peter 383
4334 25 Mar 23 peter 384   bool Aligner::Cigar::operator!=(const Aligner::Cigar& other) const
4334 25 Mar 23 peter 385   {
4334 25 Mar 23 peter 386     return cigar_ != other.cigar_;
4334 25 Mar 23 peter 387   }
4334 25 Mar 23 peter 388
4334 25 Mar 23 peter 389
3200 03 May 14 peter 390   std::ostream& operator<<(std::ostream& os, const Aligner::Cigar& cigar)
3200 03 May 14 peter 391   {
3200 03 May 14 peter 392     for (size_t i=0; i<cigar.size(); ++i)
3200 03 May 14 peter 393       os << cigar.oplen(i) << cigar.opchr(i);
3200 03 May 14 peter 394     return os;
3200 03 May 14 peter 395   }
3200 03 May 14 peter 396
2815 28 Aug 12 peter 397 }}} // of namespace utility, yat, and theplu

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