1 #include <iostream>
2 #include "chaining.h"
3 #include "../output/output_format.h"
4 #include "../dp/dp.h"
5 #include "../dp/ungapped.h"
6
7 using std::cout;
8 using std::endl;
9
print_diag(int i0,int j0,int l,int score,const Diag_graph & diags,const Sequence & query,const Sequence & subject)10 void print_diag(int i0, int j0, int l, int score, const Diag_graph &diags, const Sequence &query, const Sequence &subject)
11 {
12 Diagonal_segment ds(i0, j0, l, 0);
13 unsigned n = 0;
14 int path_max, path_min;
15 for (vector<Diagonal_node>::const_iterator d = diags.nodes.begin(); d != diags.nodes.end(); ++d) {
16 if (d->intersect(ds).len > 0) {
17 if (d->score == 0)
18 continue;
19 const int diff = score_range(query, subject, d->query_end(), d->subject_end(), j0 + l);
20 if (n > 0)
21 cout << "(";
22 cout << "Diag n=" << d - diags.nodes.begin() << " i=" << i0 << " j=" << j0 << " len=" << l
23 << " prefix_score=" << score + score_range(query, subject, i0 + l, j0 + l, d->subject_end()) - std::min(diff, 0)
24 << " prefix_score2=" << diags.prefix_score((unsigned)(d - diags.nodes.begin()), j0 + l, path_max, path_min);
25 if (n > 0)
26 cout << ")";
27 cout << endl;
28 ++n;
29 }
30 }
31 if (n == 0)
32 cout << "Diag n=x i=" << i0 << " j=" << j0 << " len=" << l << " prefix_score=" << score << endl;
33 }
34
smith_waterman(Sequence q,Sequence s,const Diag_graph & diags)35 void smith_waterman(Sequence q, Sequence s, const Diag_graph &diags)
36 {
37 Hsp hsp(true);
38 //smith_waterman(q, s, hsp);
39 Hsp::Iterator i = hsp.begin();
40 int i0 = -1, j0 = -1, l = 0, score = 0;
41 for (; i.good(); ++i) {
42 switch (i.op()) {
43 case op_match:
44 case op_substitution:
45 if (i0 < 0) {
46 i0 = i.query_pos.translated;
47 j0 = i.subject_pos;
48 l = 0;
49 }
50 score += score_matrix(q[i.query_pos.translated], s[i.subject_pos]);
51 ++l;
52 break;
53 case op_deletion:
54 case op_insertion:
55 if (i0 >= 0) {
56 print_diag(i0, j0, l, score, diags, q, s);
57 score -= score_matrix.gap_open() + score_matrix.gap_extend();
58 i0 = -1;
59 j0 = -1;
60 }
61 else
62 score -= score_matrix.gap_extend();
63 break;
64 case op_frameshift_forward:
65 case op_frameshift_reverse:
66 ;
67 }
68 }
69 print_diag(i0, j0, l, score, diags, q, s);
70 print_hsp(hsp, TranslatedSequence(q));
71 }