source: trunk/yat/utility/Aligner.cc @ 3342

Last change on this file since 3342 was 3342, checked in by Peter, 8 years ago

merge 0.12.2 into trunk

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