source: branches/0.12-stable/yat/utility/Aligner.cc @ 3339

Last change on this file since 3339 was 3339, checked in by Peter, 7 years ago

Fix problem that Aligner::member_ is not cleared to 'none' before
aligning, which means that the beginning of alignment could be
incorrect. Fixes #820.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.0 KB
Line 
1// $Id: Aligner.cc 3339 2014-11-05 23:43:01Z peter $
2
3/*
4  Copyright (C) 2012 Peter Johansson
5  Copyright (C) 2014 Jari Häkkinen, Peter Johansson
6
7  This file is part of the yat library, http://dev.thep.lu.se/yat
8
9  The yat library is free software; you can redistribute it and/or
10  modify it under the terms of the GNU General Public License as
11  published by the Free Software Foundation; either version 3 of the
12  License, or (at your option) any later version.
13
14  The yat library is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17  General Public License for more details.
18
19  You should have received a copy of the GNU General Public License
20  along with yat. If not, see <http://www.gnu.org/licenses/>.
21*/
22
23#include <config.h>
24
25#include "Aligner.h"
26#include "Cigar.h"
27#include "Matrix.h"
28
29#include <cassert>
30#include <limits>
31#include <ostream>
32
33namespace theplu {
34namespace yat {
35namespace utility {
36
37  Aligner::Aligner(double gap, double open_gap)
38    : vertical_gap_(gap), open_vertical_gap_(open_gap), horizon_gap_(gap),
39      open_horizon_gap_(open_gap)
40  {
41  }
42
43
44  Aligner::Aligner(double vertical_gap, double open_vertical_gap,
45                   double horizon_gap, double open_horizon_gap)
46    : vertical_gap_(vertical_gap), open_vertical_gap_(open_vertical_gap),
47      horizon_gap_(horizon_gap), open_horizon_gap_(open_horizon_gap)
48  {
49  }
50
51
52  double Aligner::needleman_wunsch(const Matrix& s)
53  {
54    Matrix m(s.rows()+1,s.columns()+1, -std::numeric_limits<double>::max());
55    m(0, 0) = 0.0;
56    // Init upper and left border of matrix
57    for (size_t i=1; i<m.rows(); i++)
58      m(i,0) = -i * vertical_gap_ - open_vertical_gap_;
59    for (size_t i=1; i<m.columns(); i++)
60      m(0,i) = -i * horizon_gap_ - open_horizon_gap_;
61    (*this)(s, m);
62    // return lower right corner element
63    return m(s.rows(), s.columns());
64  }
65
66
67  double Aligner::smith_waterman(const Matrix& s)
68  {
69    Matrix m(s.rows()+1,s.columns()+1);
70    (*this)(s, m);
71    return max(m);
72  }
73
74
75  void Aligner::operator()(const Matrix& score, Matrix& x)
76  {
77    assert(x.columns() == score.columns()+1);
78    assert(x.rows() == score.rows()+1);
79
80    columns_ = x.columns();
81    alignment_.assign(columns_*x.rows(), none);
82
83    for (size_t i=1; i<x.rows(); ++i)
84      for (size_t j=1; j<x.columns(); ++j) {
85        // match (or mismatch)
86        if (x(i-1, j-1) + score(i-1, j-1) > x(i,j)) {
87          x(i,j) = x(i-1, j-1) + score(i-1, j-1);
88          directions(i, j) = diagonal;
89        }
90
91        if (directions(i, j-1) == right) {
92          if (x(i, j-1)-horizon_gap_ > x(i,j)) {
93            x(i,j) = x(i,j-1) - horizon_gap_;
94            directions(i,j) = right;
95          }
96        }
97        else if(x(i,j-1)-horizon_gap_-open_horizon_gap_ > x(i,j)) {
98          x(i,j) = x(i,j-1) - horizon_gap_ - open_horizon_gap_;
99          directions(i,j) = right;
100        }
101
102        if (directions(i-1,j) == down) {
103          if (x(i-1,j) - vertical_gap_ > x(i,j)) {
104            x(i,j) = x(i-1,j) - vertical_gap_;
105            directions(i,j) = down;
106          }
107        }
108        else if (x(i-1,j) - vertical_gap_ - open_vertical_gap_ > x(i,j)) {
109          x(i,j) = x(i-1,j) - vertical_gap_ - open_vertical_gap_;
110          directions(i,j) = down;
111        }
112      }
113  }
114
115
116  const Aligner::direction& Aligner::alignment(size_t i, size_t j) const
117  {
118    return alignment_[i*columns_ + j];
119  }
120
121
122  const Aligner::Cigar Aligner::cigar(size_t row, size_t column) const
123  {
124    Cigar cig;
125    while (row || column) {
126      switch (alignment(row, column)) {
127      case(diagonal):
128        assert(row);
129        --row;
130        assert(column);
131        --column;
132        cig.push_front(BAM_CMATCH);
133        break;
134      case (right):
135        assert(column);
136        --column;
137        cig.push_front(BAM_CINS);
138        break;
139      case (down):
140        assert(row);
141        --row;
142        cig.push_front(BAM_CDEL);
143        break;
144      default:
145        return cig;
146      }
147    }
148    return cig;
149  }
150
151
152  Aligner::direction& Aligner::directions(size_t i, size_t j)
153  {
154    return alignment_[i*columns_ + j];
155  }
156
157
158  // ================ Cigar =======================
159
160  void Aligner::Cigar::clear(void)
161  {
162    cigar_.clear();
163  }
164
165
166  uint8_t Aligner::Cigar::op(size_t i) const
167  {
168    assert(i<cigar_.size());
169    return bam_cigar_op(cigar_[i]);
170  }
171
172
173  uint32_t Aligner::Cigar::length(uint8_t mask) const
174  {
175    uint32_t res = 0;
176    for (size_t i=0; i<cigar_.size(); ++i)
177      if (bam_cigar_type(op(i)) & mask)
178        res += oplen(i);
179    return res;
180  }
181
182
183  uint32_t Aligner::Cigar::length(void) const
184  {
185    return length(3);
186  }
187
188
189  char Aligner::Cigar::opchr(size_t i) const
190  {
191    assert(i<cigar_.size());
192    return bam_cigar_opchr(cigar_[i]);
193  }
194
195
196  uint32_t Aligner::Cigar::oplen(size_t i) const
197  {
198    assert(i<cigar_.size());
199    return bam_cigar_oplen(cigar_[i]);
200  }
201
202
203  void Aligner::Cigar::pop_back(uint32_t len)
204  {
205    while (len) {
206      assert(size());
207      size_t i = size()-1;
208      if (len < oplen(i)) {
209        cigar_[i] -= (len<<BAM_CIGAR_SHIFT);
210        return;
211      }
212      else {
213        len -= oplen(i);
214        cigar_.pop_back();
215      }
216    }
217  }
218
219
220  void Aligner::Cigar::pop_front(uint32_t len)
221  {
222    while (len) {
223      assert(size());
224      if (len < oplen(0)) {
225        cigar_[0] -= (len<<BAM_CIGAR_SHIFT);
226        return;
227      }
228      else {
229        len -= oplen(0);
230        cigar_.pop_front();
231      }
232    }
233  }
234
235
236  void Aligner::Cigar::push_back(uint8_t op, uint32_t len)
237  {
238    if (len==0)
239      return;
240    if (cigar_.empty() || this->op(size()-1)!=op)
241      cigar_.push_back(bam_cigar_gen(len, op));
242    else
243      cigar_[size()-1] += (len<<BAM_CIGAR_SHIFT);
244  }
245
246
247  void Aligner::Cigar::push_front(uint8_t op, uint32_t len)
248  {
249    if (len==0)
250      return;
251    if (cigar_.empty() || this->op(0)!=op)
252      cigar_.push_front(bam_cigar_gen(len, op));
253    else
254      cigar_[0] += (len<<BAM_CIGAR_SHIFT);
255  }
256
257
258  uint32_t Aligner::Cigar::query_length(void) const
259  {
260    return length(1);
261  }
262
263
264  uint32_t Aligner::Cigar::reference_length(void) const
265  {
266    return length(2);
267  }
268
269
270  size_t Aligner::Cigar::size(void) const
271  {
272    return cigar_.size();
273  }
274
275
276  uint32_t Aligner::Cigar::operator[](size_t i) const
277  {
278    assert(i<cigar_.size());
279    return cigar_[i];
280  }
281
282
283  std::ostream& operator<<(std::ostream& os, const Aligner::Cigar& cigar)
284  {
285    for (size_t i=0; i<cigar.size(); ++i)
286      os << cigar.oplen(i) << cigar.opchr(i);
287    return os;
288  }
289
290}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.