4353 |
23 Aug 23 |
peter |
// $Id$ |
4353 |
23 Aug 23 |
peter |
2 |
|
4353 |
23 Aug 23 |
peter |
3 |
/* |
4353 |
23 Aug 23 |
peter |
Copyright (C) 2023 Peter Johansson |
4353 |
23 Aug 23 |
peter |
5 |
|
4353 |
23 Aug 23 |
peter |
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 |
The yat library is free software; you can redistribute it and/or |
4353 |
23 Aug 23 |
peter |
modify it under the terms of the GNU General Public License as |
4353 |
23 Aug 23 |
peter |
published by the Free Software Foundation; either version 3 of the |
4353 |
23 Aug 23 |
peter |
License, or (at your option) any later version. |
4353 |
23 Aug 23 |
peter |
12 |
|
4353 |
23 Aug 23 |
peter |
The yat library is distributed in the hope that it will be useful, |
4353 |
23 Aug 23 |
peter |
but WITHOUT ANY WARRANTY; without even the implied warranty of |
4353 |
23 Aug 23 |
peter |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
4353 |
23 Aug 23 |
peter |
General Public License for more details. |
4353 |
23 Aug 23 |
peter |
17 |
|
4353 |
23 Aug 23 |
peter |
You should have received a copy of the GNU General Public License |
4353 |
23 Aug 23 |
peter |
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 |
// 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 |
// 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 |
// 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 |
}}} |