ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/PrimerMatch/pattern_alignment.cc
Revision: 1.2
Committed: Wed May 4 18:03:45 2005 UTC (11 years, 7 months ago) by nje01
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +130 -46 lines
Log Message:
Small bug fixes, plus codon based edit distance for peptide searching.

Line File contents
1 /**************************************************************************
2 * This code is part of the supporting infrastructure for ATA Mapper.
3 * Copyright (C) 2002,2003,2004 Applera Corporation. All rights reserved.
4 * Author: Nathan Edwards
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received (LICENSE.txt) a copy of the GNU General Public
17 * License along with this program; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *************************************************************************/
20
21
22 #include "pattern_alignment.h"
23 #include "util.h"
24
25 #include <list>
26 #include <iomanip>
27
28 #if !defined(NO_STD_NAMESPACE)
29 using namespace std;
30 #endif
31
32 istream & operator>>(istream & is, pattern_alignment & pa) {
33 pa.read(is);
34 return is;
35 }
36
37 ostream & operator<<(ostream & os, pattern_alignment const & pa) {
38 pa.write(os);
39 return os;
40 }
41
42 pattern_alignment::~pattern_alignment() {}
43
44 void pattern_alignment::write(ostream & os) const {
45 os << pattern_ << " " << end_;
46 }
47
48
49 bool exact_alignment::align(CharacterProducer & /*cp*/, std::string const & pat) {
50 if (yesno_) {
51 return true;
52 }
53 matching_text_ = pat;
54 start_ = end() - matching_text_.length();
55 alignment_.resize(matching_text_.length());
56 for (int i=0;i<matching_text_.length();i++) {
57 alignment_[i] = alignment_equal;
58 }
59 stats_[alignment_equal] = matching_text_.length();
60 alignment_done_ = true;
61 return true;
62 }
63
64 bool exact_peptide_alignment::align(CharacterProducer & cp, std::string const & pat) {
65 if (yesno_) {
66 return true;
67 }
68 matching_text_ = pat;
69 start_ = end() - matching_text_.length();
70 alignment_.resize(matching_text_.length());
71 for (int i=0;i<matching_text_.length();i++) {
72 alignment_[i] = alignment_equal;
73 }
74 stats_[alignment_equal] = matching_text_.length();
75 cp.pos(start_-1);
76 lcontext_ = cp.getch();
77 cp.pos(end());
78 rcontext_ = cp.getch();
79 alignment_done_ = true;
80 return true;
81 }
82
83 bool exact_wc_alignment::align(CharacterProducer & cp, std::string const & pat) {
84 start_ = end() - pat.length();
85 alignment_.resize(pat.length());
86 char *buffer = new char[pat.length()+1];
87 cp.pos(start_);
88 for (int i=0;i<pat.length();i++) {
89 buffer[i] = cp.getch();
90 if (buffer[i]==pat[i]) {
91 alignment_[i] = alignment_equal;
92 } else if (iupac_compatible(buffer[i],pat[i]) &&
93 (textn_ || buffer[i] != 'N')) {
94 alignment_[i] = alignment_wildcard_equal;
95 } else {
96 alignment_[i] = alignment_substitution;
97 }
98 stats_[alignment_[i]]++;
99 }
100 buffer[pat.length()] = '\0';
101 matching_text_ = buffer;
102 delete [] buffer;
103 alignment_done_ = true;
104 if (editdist()<=0) return true;
105 return false;
106 }
107
108 bool mismatch_alignment::align(CharacterProducer & cp, std::string const & pat) {
109 start_ = end() - pat.length();
110 alignment_.resize(pat.length());
111 char *buffer = new char[pat.length()+1];
112 cp.pos(start_);
113 for (int i=0;i<pat.length();i++) {
114 buffer[i] = cp.getch();
115 if (buffer[i]==pat[i]) {
116 alignment_[i] = alignment_equal;
117 } else {
118 alignment_[i] = alignment_substitution;
119 }
120 stats_[alignment_[i]]++;
121 }
122 buffer[pat.length()] = '\0';
123 matching_text_ = buffer;
124 delete [] buffer;
125 alignment_done_ = true;
126 return true;
127 }
128
129 bool editdist_alignment::align(CharacterProducer & cp, std::string const & pat) {
130 // const bool problem_case(end() <= 217834338 && 217834338 <= end2_);
131 // if (problem_case) {
132 // checkpoint;
133 // cerr << "end():" << end() << endl;
134 // cerr << "end2_:" << end2_ << endl;
135 // cerr << "pat.length():" << pat.length() << endl;
136 // cerr << "k_:" << k_ << endl;
137 // cerr << "wc_:" << ((wc_)?"true":"false") << endl;
138 // cerr << "1indels_:" << ((indels_)?"true":"false") << endl;
139 // cerr << "yesno_: " << ((yesno_)?"true":"false") << endl;
140 // cerr << "eos_: " << (int)eos_ << endl;
141 // }
142 int const_viol_penalty = 5*k_+1;
143 FILE_POSITION_TYPE textstart=0;
144 if (end() > pat.length() + k_) {
145 textstart = end() - pat.length() - k_;
146 }
147 // if (problem_case) checkpoint;
148 // if (problem_case) cerr << textstart << endl;
149 FILE_SIZE_TYPE buflen = end2_ - textstart;
150 if (maxpatlen_ < pat.length()) {
151 maxpatlen_ = pat.length();
152 }
153 // if (problem_case) cerr << buflen << endl;
154 bool realloc=false;
155 if (!buffer_ || (buflen+1)>bufsize_) {
156 // if (problem_case) checkpoint;
157 if (buffer_) {
158 delete [] buffer_;
159 }
160 // size should always be larger than buflen+1;
161 bufsize_ = 2*maxpatlen_+((indels_)?(dna_mut_?1:k_):0)+1;
162 if (bufsize_<(buflen+1)) {
163 bufsize_ = (buflen+1);
164 }
165 // if (problem_case) cerr << maxpatlen_ << " " << bufsize_ << endl;
166 buffer_ = new char[bufsize_+1];
167 realloc=true;
168 // if (problem_case) checkpoint;
169 }
170 // if (problem_case) checkpoint;
171 char *buffer;
172 if (textstart < bufstart_ || textstart + buflen >= bufend_ || realloc) {
173 // if (problem_case) checkpoint;
174 cp.pos(textstart);
175 for (int i=0;i<bufsize_;i++) {
176 buffer_[i] = cp.getch();
177 }
178 buffer_[bufsize_]='\0';
179 // if (problem_case) cerr << bufsize_ << " " << buffer_ << endl;
180 bufstart_ = textstart;
181 bufend_ = textstart+bufsize_;
182 buffer = buffer_;
183 } else {
184 // if (problem_case) checkpoint;
185 buffer = buffer_+(textstart-bufstart_);
186 }
187 // if (problem_case) checkpoint;
188 // if (problem_case) cerr << pat.length() << " " << pat << endl;
189 // if (problem_case) cerr << buflen << " ";
190 // if (problem_case) for (int i=0;i<buflen;i++)
191 // if (problem_case) cerr << buffer[i];
192 // if (problem_case) cerr << endl;
193
194 // if (problem_case) checkpoint;
195 if (!dp_ ||
196 matsize_ < (buflen+1)*(pat.length()+1)) {
197 if (dp_) {
198 delete [] dp_;
199 }
200 if (best_) {
201 delete [] best_;
202 }
203 // if (problem_case) checkpoint;
204 matsize_ = (bufsize_+1)*(maxpatlen_+1);
205 if (matsize_ < (buflen+1)*(pat.length()+1)) {
206 matsize_ = (buflen+1)*(pat.length()+1);
207 }
208 // if (problem_case) cerr << matsize_ << endl;
209 dp_ = new unsigned int[matsize_];
210 // if (problem_case) checkpoint;
211 best_ = new int[matsize_];
212 // if (problem_case) checkpoint;
213 }
214 // if (problem_case) checkpoint;
215 int lbexact=0;
216 int rbexact=pat.length()+1;
217 if (lconst_ > 0) rbexact = pat.length()+1-lconst_;
218 if (rconst_ > 0) lbexact = rconst_;
219
220 // if (problem_case) checkpoint;
221 // if (problem_case) cerr << lbexact << " " << rbexact << endl;
222
223 unsigned int *dp = dp_;
224 int *best = best_;
225 #define index(i,j) ((i)*(buflen+1)+(j))
226 int p,t;
227 dp[index(0,0)] = 0;
228 best[index(0,0)] = alignment_mask_alignment_end;
229 // if (problem_case) cerr << "dp " << 0 << " " << 0 << " " << dp[index(0,0)] << endl;
230 // if (problem_case) cerr << "bt " << 0 << " " << 0 << " " << best[index(0,0)] << endl;
231 int lb,ub;
232 ub = (indels_?(dna_mut_?1:k_):0);
233 if (ub > pat.length()) {
234 ub = pat.length();
235 }
236 // if (problem_case) checkpoint;
237 // if (problem_case) cerr << "ub: " << ub << endl;
238 for (p=1;p<=ub;p++) {
239 if (!indels_ || p < lbexact || p >= rbexact) {
240 dp[index(p,0)] = const_viol_penalty;
241 best[index(p,0)] = alignment_mask_constraint_violation;
242 } else {
243 if (!dna_mut_) {
244 dp[index(p,0)] = dp[index(p-1,0)] + 1; /* cost for a deletion */
245 best[index(p,0)] = alignment_mask_deletion;
246 } else {
247 dp[index(p,0)] = dp[index(p-1,0)] + 3; /* cost for a deletion */
248 best[index(p,0)] = alignment_mask_deletion_3;
249 }
250 }
251 // if (problem_case) cerr << "dp " << p << " " << 0 << " " << dp[index(p,0)] << endl;
252 // if (problem_case) cerr << "bt " << p << " " << 0 << " " << best[index(p,0)] << endl;
253 }
254 ub = end2_-end()+(indels_?(dna_mut_?1:k_):0);
255 if (buflen < ub) {
256 ub = buflen;
257 }
258 // if (problem_case) checkpoint;
259 // if (problem_case) cerr << "ub: " << ub << endl;
260 // if (problem_case) cerr << "end2_ - end(): " << end2_ - end() << endl;
261 for (t=1;t<=ub;t++) {
262 if (t <= end2_ - end()) {
263 dp[index(0,t)] = 0;
264 best[index(0,t)] = alignment_mask_alignment_end;
265 } else if (!indels_ || lbexact > 0) {
266 dp[index(0,t)] = const_viol_penalty;
267 best[index(0,t)] = alignment_mask_constraint_violation;
268 } else {
269 if (!dna_mut_) {
270 dp[index(0,t)] = dp[index(0,t-1)] + 1; /* cost for a insertion */
271 best[index(0,t)] = alignment_mask_insertion;
272 } else {
273 dp[index(0,t)] = dp[index(0,t-1)] + 3; /* cost for a insertion */
274 best[index(0,t)] = alignment_mask_insertion_3;
275 }
276 }
277 // if (problem_case) cerr << "dp " << 0 << " " << t << " " << index(0,t) << " " << dp[index(0,t)] << endl;
278 // if (problem_case) cerr << "bt " << 0 << " " << t << " " << index(0,t) << " " << best[index(0,t)] << endl;
279 }
280 // if (problem_case) checkpoint;
281 for (p=1;p<=pat.length();p++) {
282 lb = p-(indels_?(dna_mut_?1:k_):0);
283 if (lb < 1) {
284 lb = 1;
285 }
286 ub = p+end2_-end()+(indels_?(dna_mut_?1:k_):0);
287 if (buflen < ub) {
288 ub = buflen;
289 }
290 // if (problem_case) checkpoint;
291 // cerr << 0 << " " << 27 << " " << index(0,27) << " " << dp[index(0,27)] << endl;
292 // cerr << 0 << " " << 27 << " " << index(0,27) << " " << best[index(0,27)] << endl;
293 // if (problem_case) cerr << "lb: " << lb << endl;
294 // if (problem_case) cerr << "ub: " << ub << endl;
295 int bestscorerow=const_viol_penalty;
296 for (t=lb;t<=ub;t++) {
297 unsigned int v,v1;
298 int ac=0;
299 if (buffer[buflen-t] == pat[pat.length()-p]) {
300 v = dp[index(p-1,t-1)]; /* characters match, no charge */
301 ac = alignment_mask_equal;
302 } else if (wc_ && iupac_compatible(pat[pat.length()-p],buffer[buflen-t]) && (buffer[buflen-t]!='N' || textn_)) {
303 v = dp[index(p-1,t-1)]; /* characters match, no charge */
304 ac = alignment_mask_wildcard_equal;
305 } else if (buffer[buflen-t] == eos_ ||
306 p <= lbexact || p >=rbexact) {
307 v = const_viol_penalty;
308 ac = alignment_mask_constraint_violation;
309 } else {
310 int mut=-1;
311 if (!dna_mut_) {
312 v = dp[index(p-1,t-1)]+1; /* cost for a substitution */
313 ac = alignment_mask_substitution;
314 } else if ((mut=aasubdist(buffer[buflen-t],pat[pat.length()-p]))>=0) {
315 v = dp[index(p-1,t-1)]+mut; /* cost for subst*/
316 switch(mut) {
317 case 1:
318 ac = alignment_mask_substitution_1;
319 break;
320 case 2:
321 ac = alignment_mask_substitution_2;
322 break;
323 case 3:
324 ac = alignment_mask_substitution_3;
325 break;
326 default:
327 checkpoint;
328 assert(0);
329 }
330 } else {
331 v = const_viol_penalty;
332 ac = alignment_mask_constraint_violation;
333 }
334 }
335 if (buffer[buflen-t] == eos_ || !indels_ || t <= lb ||
336 p < lbexact || p >= rbexact) {
337 v1 = const_viol_penalty; /* cost for an insertion */
338 if (v1 < v) {
339 v = v1;
340 ac = alignment_mask_constraint_violation;
341 }
342 } else {
343 if (!dna_mut_) {
344 v1 = dp[index(p,t-1)]+1; /* cost for an insertion */
345 if (v1 < v) {
346 v = v1;
347 ac = alignment_mask_insertion;
348 } else if (v1 == v) {
349 ac |= alignment_mask_insertion;
350 }
351 } else {
352 v1 = dp[index(p,t-1)]+3; /* cost for an insertion */
353 if (v1 < v) {
354 v = v1;
355 ac = alignment_mask_insertion_3;
356 } else if (v1 == v) {
357 ac |= alignment_mask_insertion_3;
358 }
359 }
360 }
361 if (!indels_ || t >= ub || p <= lbexact || p >= rbexact) {
362 v1 = const_viol_penalty; /* cost for a deletion */
363 if (v1 < v) {
364 v = v1;
365 ac = alignment_mask_constraint_violation;
366 }
367 } else {
368 if (!dna_mut_) {
369 v1 = dp[index(p-1,t)]+1; /* cost for a deletion */
370 if (v1 < v) {
371 v = v1;
372 ac = alignment_mask_deletion;
373 } else if (v1 == v) {
374 ac |= alignment_mask_deletion;
375 }
376 } else {
377 v1 = dp[index(p-1,t)]+3; /* cost for a deletion */
378 if (v1 < v) {
379 v = v1;
380 ac = alignment_mask_deletion_3;
381 } else if (v1 == v) {
382 ac |= alignment_mask_deletion_3;
383 }
384 }
385 }
386 dp[index(p,t)] = v;
387 best[index(p,t)] = ac;
388 // if (problem_case) cerr << "dp " << p << " " << t << " " << index(p,t) << " " << dp[index(p,t)] << endl;
389 // if (problem_case) cerr << "bt " << p << " " << t << " " << index(p,t) << " " << best[index(p,t)] << endl;
390 if (v < bestscorerow) {
391 bestscorerow = v;
392 }
393 }
394 if (bestscorerow > k_) {
395 // if (problem_case) checkpoint;
396 // if (problem_case) cerr << bestscorerow << " " << k_ << endl;
397 if (!yesno_) {
398 // checkpoint;
399 alignment_.push_back(alignment_constraint_violation);
400 stats_[alignment_constraint_violation]++;
401 }
402 alignment_done_ = true;
403 // if (problem_case) checkpoint;
404 return false;
405 }
406 }
407
408 // cerr << 0 << " " << 27 << " " << index(0,27) << " " << dp[index(0,27)] << endl;
409 // cerr << 0 << " " << 27 << " " << index(0,27) << " " << best[index(0,27)] << endl;
410 // if (problem_case) checkpoint;
411
412 int beststart=pat.length()-(indels_?(dna_mut_?1:k_):0);
413 if (beststart>buflen) {
414 beststart = buflen;
415 }
416 if (beststart < 0) {
417 beststart = 0;
418 }
419 // cerr << "beststart: " << beststart << endl;
420 int bestval=dp[index(pat.length(),beststart)];
421 ub = pat.length()+end2_-end()+(indels_?(dna_mut_?1:k_):0);
422 if (buflen < ub) {
423 ub = buflen;
424 }
425 // checkpoint;
426 // cerr << "ub: " << ub << endl;
427 // if (problem_case) checkpoint;
428 // if (problem_case) cerr << bestval << " " << beststart << endl;
429 // if (problem_case) cerr << dp[index(pat.length(),beststart)] << " " << beststart << endl;
430 for (int k=beststart+1;k<=ub;k++) {
431 // if (problem_case) cerr << dp[index(pat.length(),k)] << " " << k << endl;
432 if ((dp[index(pat.length(),k)]) < bestval ||
433 ((dp[index(pat.length(),k)]) <= bestval &&
434 (best[index(pat.length(),k)]&(alignment_mask_equal|
435 alignment_mask_wildcard_equal|
436 alignment_mask_substitution|
437 alignment_mask_substitution_1|
438 alignment_mask_substitution_2|
439 alignment_mask_substitution_3)
440 ))) {
441 bestval = dp[index(pat.length(),k)];
442 beststart = k;
443 }
444 }
445
446 // if (problem_case) checkpoint;
447 // if (problem_case) cerr << bestval << " " << beststart << endl;
448
449 p=pat.length(); t=beststart;
450 // if (problem_case) cerr << p << " " << t << endl;
451 if (t < p - (indels_?(dna_mut_?1:k_):0) || t > p + (indels_?(dna_mut_?1:k_):0) + (end2_ - end())) {
452 if (!yesno_) {
453 alignment_.push_back(alignment_constraint_violation);
454 stats_[alignment_constraint_violation]++;
455 }
456 // if (problem_case) checkpoint;
457 alignment_done_ = true;
458 return false;
459 }
460
461 // if (problem_case) checkpoint;
462
463 std::list<alignment_code> alignment;
464 std::list<alignment_code>::iterator alit;
465
466 // if (problem_case) checkpoint;
467
468 int ac;
469 bool match,ins,del,wc,sub;
470 // checkpoint;
471
472 // cerr << 0 << " " << 27 << " " << index(0,27) << " " << dp[index(0,27)] << endl;
473 // cerr << 0 << " " << 27 << " " << index(0,27) << " " << best[index(0,27)] << endl;
474
475 assert(dp[index(0,0)] == 0);
476 assert(best[index(0,0)] == alignment_mask_alignment_end);
477
478 alignment_code lastac=alignment_none;
479 // checkpoint;
480 // cerr << p << " " << t << endl;
481 // cerr << p << " " << t << " " << dp[index(p,t)] << endl;
482 // cerr << p << " " << t << " " << best[index(p,t)] << endl;
483 while (!(best[index(p,t)] & alignment_mask_alignment_end)) {
484 // if (problem_case) cerr << p << " " << t << " " << index(p,t) << " " << dp[index(p,t)] << endl;
485 // if (problem_case) cerr << p << " " << t << " " << index(p,t) << " " << best[index(p,t)] << endl;
486 assert(p >= 0 && t >= 0);
487 ac = best[index(p,t)];
488 match = (ac & (alignment_mask_equal|
489 alignment_mask_wildcard_equal|
490 alignment_mask_substitution|
491 alignment_mask_substitution_1|
492 alignment_mask_substitution_2|
493 alignment_mask_substitution_3)
494 );
495 wc = (ac & alignment_mask_wildcard_equal);
496 sub = (ac & (alignment_mask_substitution|
497 alignment_mask_substitution_1|
498 alignment_mask_substitution_2|
499 alignment_mask_substitution_3)
500 );
501 ins = (ac & (alignment_mask_insertion|
502 alignment_mask_insertion_3));
503 del = (ac & (alignment_mask_deletion|
504 alignment_mask_deletion_3));
505 // cerr << "eodp: " << ((ac & alignment_mask_alignment_end)?"true":"false") << " ";
506 // cerr << "match: " << (match?"true":"false") << " ";
507 // cerr << "wc: " << (wc?"true":"false") << " ";
508 // cerr << "sub: " << (sub?"true":"false") << " ";
509 // cerr << "ins: " << (ins?"true":"false") << " ";
510 // cerr << "del: " << (del?"true":"false") << " ";
511 // cerr << endl;
512 if (match && !(((lastac == alignment_insertion ||
513 lastac == alignment_insertion_3) && ins) ||
514 ((lastac == alignment_deletion ||
515 lastac == alignment_deletion_3) && del) ||
516 (lastac == alignment_wildcard_equal && !wc && (ins || del)))) {
517 p--; t--;
518 if ((ac & alignment_mask_equal) &&
519 !((lastac == alignment_wildcard_equal && wc) ||
520 (lastac == alignment_substitution && sub))) {
521 lastac = alignment_equal;
522 } else if (wc) {
523 lastac = alignment_wildcard_equal;
524 } else if (sub) {
525 if (ac & alignment_mask_substitution) {
526 lastac = alignment_substitution;
527 } else if (ac & alignment_mask_substitution_1) {
528 lastac = alignment_substitution_1;
529 } else if (ac & alignment_mask_substitution_2) {
530 lastac = alignment_substitution_2;
531 } else if (ac & alignment_mask_substitution_3) {
532 lastac = alignment_substitution_3;
533 }
534 }
535 } else if (del) {
536 p--;
537 if (ac & alignment_mask_deletion) {
538 lastac = alignment_deletion;
539 } else if (ac & alignment_mask_deletion_3) {
540 lastac = alignment_deletion_3;
541 }
542 } else if (ins) {
543 t--;
544 if (ac & alignment_mask_insertion) {
545 lastac = alignment_insertion;
546 } else if (ac & alignment_mask_insertion_3) {
547 lastac = alignment_insertion_3;
548 }
549 } else if (ac & alignment_mask_constraint_violation) {
550 p = 0; t = 0;
551 lastac = alignment_constraint_violation;
552 } else {
553 assert(0);
554 }
555 if (!yesno_) {
556 stats_[lastac]++;
557 alignment.push_back(lastac);
558 }
559 }
560 if (!yesno_) {
561 alignment_.resize(alignment.size());
562 for (int i=0;i<alignment_.size();i++) {
563 alignment_[i] = alignment.front();
564 alignment.pop_front();
565 }
566 }
567 int endt = t;
568 int endp = p;
569
570 // if (problem_case) cerr << buflen << " " << beststart << " " << endt << endl;
571
572 matching_text_ = "";
573 for (int i=0;i<(buflen-(buflen-beststart)-endt);i++) {
574 matching_text_ += *(buffer+(buflen-beststart)+i);
575 }
576 // if (problem_case) cerr << matching_text_ << endl;
577 start_ = end2_ - beststart;
578 end(start_ + matching_text_.length());
579 // if (problem_case) cerr << start_ << endl;
580 // if (problem_case) cerr << end() << endl;
581
582 // if (problem_case) checkpoint;
583
584 if (/*problem_case || */ endp < 0 || endt < 0) {
585 cerr << "Something is wrong with the dp..." << endl;
586 // Dump everything we know about the DP!
587 checkpoint;
588 cerr << "Abort DP: Dumping everything!\n";
589 cerr << "Text: ";
590 for (t=0;t<buflen;t++) {
591 cerr << buffer[t];
592 }
593 cerr << endl;
594 cerr << "Pattern: " << pat << endl;
595 cerr << "DP Matrix: " << endl;
596 for (t=0;t<=buflen;t++) {
597 if (t>0) {
598 cerr << setw(3) << (char)buffer[buflen-t] << " ";
599 } else {
600 cerr << " ";
601 }
602 }
603 cerr << endl;
604 // cerr << end2_ << " " << end() << endl;
605 for (p=0;p<=pat.length();p++) {
606 lb = p-(indels_?(dna_mut_?1:k_):0);
607 if (lb < 0) {
608 lb = 0;
609 }
610 ub = p+end2_-end()+(indels_?(dna_mut_?1:k_):0);
611 if (buflen < ub) {
612 ub = buflen;
613 }
614 // cerr << "p = " << p << ", t lb = " << lb << " ub = " << ub << endl;
615 if (p > 0) {
616 cerr << pat[pat.length()-p] << ": ";
617 } else {
618 cerr << " : ";
619 }
620 for (t=0;t<lb;t++) {
621 cerr << " ";
622 }
623 for (t=lb;t<=ub;t++) {
624 cerr << setw(3) << dp[index(p,t)] << " ";
625 }
626 cerr << endl;
627 }
628 cerr << "Best Matrix: " << endl;
629 for (t=0;t<=buflen;t++) {
630 if (t>0) {
631 cerr << setw(3) << (char)buffer[buflen-t] << " ";
632 } else {
633 cerr << " ";
634 }
635 }
636 cerr << endl;
637 for (p=0;p<=pat.length();p++) {
638 lb = p-(indels_?(dna_mut_?1:k_):0);
639 if (lb < 0) {
640 lb = 0;
641 }
642 ub = p+end2_-end()+(indels_?(dna_mut_?1:k_):0);
643 if (buflen < ub) {
644 ub = buflen;
645 }
646 // cerr << "p = " << p << ", t lb = " << lb << " ub = " << ub << endl;
647 if (p > 0) {
648 cerr << pat[pat.length()-p] << ": ";
649 } else {
650 cerr << " : ";
651 }
652 for (t=0;t<lb;t++) {
653 cerr << " ";
654 }
655 for (t=lb;t<=ub;t++) {
656 cerr << setw(3) << best[index(p,t)] << " "; }
657 cerr << endl;
658 }
659 cerr << "Bestval: " << bestval << endl;
660 cerr << "Bestpos: " << beststart << endl;
661 cerr << "Matching Text: " << matching_text_ << endl;
662 cerr << "Alignment chars: ";
663
664 for (int i=0;i<alignment_.size();i++) {
665 cerr << alignment_[i] << " ";
666 }
667 cerr << endl;
668 cerr << "Current p,t: " << endp << "," << endt << endl;
669 // abort();
670 }
671 alignment_done_ = true;
672 return (bestval<=k_);
673 }
674
675 // void alignment_write(ostream & os,
676 // FILE_POSITION_TYPE const seq_pos,
677 // std::string const & pattern,
678 // long unsigned int id, bool revcomp) {
679 // assert(alignment_done_);
680 // std::string const & mt = matching_text();
681 // int p=0;
682 // os << " ";
683 // for (int i=0; i<size(); i++) {
684 // if ((*this)[i] != alignment_deletion) {
685 // os << mt[p];
686 // p++;
687 // } else {
688 // os << "-";
689 // }
690 // }
691 // os /* << " " << start() << " " << end() */
692 // << " " << seq_pos - length() + 1 << " " << seq_pos
693 // << " " << editdist();
694 // os << endl;
695 // os << " ";
696 // for (int i=0; i<size(); i++) {
697 // if ((*this)[i] == alignment_equal) {
698 // os << "|";
699 // p++;
700 // } else if ((*this)[i] == alignment_substitution) {
701 // os << "*";
702 // } else {
703 // os << " ";
704 // }
705 // }
706 // os << endl;
707 // std::string kw = pattern;
708 // p = 0;
709 // os << " ";
710 // for (int i=0; i<size(); i++) {
711 // if ((*this)[i] != alignment_insertion) {
712 // os << kw[p];
713 // p++;
714 // } else {
715 // os << "-";
716 // }
717 // }
718 // os << " " << id;
719 // if (revcomp) {
720 // os << " REVCOMP";
721 // }
722 // os << endl;
723 // }