+ added c code from the palma project
[qpalma.git] / QPalmaDP / result_align.cpp
1 #include "assert.h"
2 #include "fill_matrix.h"
3
4 using namespace std;
5
6 bool result_align(Pre_score* matrices[], int z, int est_len, int dna_len, int* result_length_ptr, char* est, char* dna, int* s_align, int* e_align, int* mparam, double* alignmentscores, int* max_score_positions)
7 //Gibt 1 zurueck, wenn kein weiterer Pfad existiert, sonst 0
8
9 {
10 const int mlen=6; // length of matchmatrix
11 int startend[4] ; // [i_start, j_start, i_end, j_end] ;
12 startend[3] = max_score_positions[2*z] ; //i (est)
13 startend[4] = max_score_positions[2*z +1] ; //j (dna)
14 //printf("in result_align: z=%d:(%d, %d)\n", z, startend[3], startend[4]) ;
15 int dnanum ; //using check_char: letter -> number
16 int estnum ; //using check_char: letter -> number
17
18 int dna_pos ;
19 int est_pos ;
20 int splice_state ;
21 int est_state ;
22 for (dna_pos=dna_len-1; dna_pos>startend[4]; dna_pos--) {
23 s_align[dna_pos-1] = 4 ; // not aligned (-1 because array starts with index 0)
24 }
25 for (est_pos=est_len-1; est_pos>startend[3]; est_pos--) {
26 e_align[est_pos-1] = 4 ;
27 }
28 //splice_align: 0: exon, 1: donor, 2: acceptor, 3: intro,. 4: dangling end
29
30
31 /* starting point for backtracking */
32 int i = startend[3] ;
33 int j = startend[4] ;
34 int act_m_no = z ;
35 alignmentscores[z] = ((Pre_score*)matrices[z] + i*dna_len +j)->value;
36 //printf("alignment_score for z=%d: %3.8f\n", z, alignmentscores[z]) ;
37
38 int prev_i = ((Pre_score*)matrices[act_m_no] + i*dna_len +j)->prev_i;
39 int prev_j = ((Pre_score*)matrices[act_m_no] + i*dna_len +j)->prev_j;
40 int prev_m_no = ((Pre_score*)matrices[act_m_no] + i*dna_len +j)->prev_matrix_no;
41
42 if(!(isfinite(((Pre_score*)matrices[act_m_no] + i*dna_len +j)->value)))
43 {
44 cout << "Kein weiterer Pfad vorhanden\n";
45 return 1;
46 }
47
48 if(((Pre_score*)matrices[act_m_no] + i*dna_len +j)->value==0) {
49 cout << "Kein weiterer Pfad vorhanden\n";
50 return 1;
51 }
52
53 assert(dna_pos == j) ;
54 assert(est_pos == i) ;
55 splice_state = 0 ; //exon
56 est_state = 1 ; //[112222111 ... for exons of length 2,4,3 ...
57
58 //compute length of alignment (between max(m,n) and m+n)
59 //local_alignment: backtracking to first Zero: look at actual value
60 while(!(((Pre_score*)matrices[act_m_no] + i*dna_len +j)->value)==0.0) {
61
62 //just some info
63 //if ((z>0) && (!(z==prev_m_no))) {
64 // printf("(%d,%d): next step in next matrix\n", i,j) ;
65 //}
66
67 if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
68 (*result_length_ptr)= (*result_length_ptr) + 1;
69
70 s_align[dna_pos-1] = splice_state;
71 dna_pos-- ;
72 e_align[est_pos-1] = est_state ; //1 or 2, depended
73 est_pos-- ;
74
75 dnanum = check_char(dna[j-1]) ; //-1 because index starts with 0
76 estnum = check_char(est[i-1]) ; //-1 because index starts with 0
77 mparam[mlen*dnanum +estnum] ++ ;
78 }
79 else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST
80 (*result_length_ptr)= (*result_length_ptr) + 1;
81
82 s_align[dna_pos-1] = splice_state;
83 dna_pos-- ;
84
85 dnanum = check_char(dna[j-1]) ; //-1 because index starts with 0
86 estnum = 0 ; //gap
87 mparam[mlen*dnanum +estnum] ++ ;
88 }
89 else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA
90 (*result_length_ptr)= (*result_length_ptr) + 1;
91
92 e_align[est_pos-1] = est_state ; //1 or 2, depended
93 est_pos-- ;
94
95 dnanum = 0 ; //gap
96 estnum = check_char(est[i-1]) ; //-1 because index starts with 0
97 mparam[mlen*dnanum +estnum] ++ ;
98 }
99 else {// splice site
100 (*result_length_ptr) = (*result_length_ptr) + (j - prev_j);
101 if (((Pre_score*)matrices[act_m_no] + i*dna_len + j)->issplice==0) {
102 printf("matrices[act_m_no](i,j): %3.15f\n",((Pre_score*)matrices[act_m_no] + i*dna_len +j)->value) ;
103 printf("matrices[prev_m_no](i,j): %3.15f\n",((Pre_score*)matrices[prev_m_no] + i*dna_len +j)->value) ;
104 printf("issplice(act): %d\n",((Pre_score*)matrices[act_m_no] + i*dna_len + j)->issplice) ;
105 printf("issplice(prev): %d\n",((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->issplice) ;
106 printf("act_m_no: %d, prev_m_no: %d\n", act_m_no, prev_m_no) ;
107 printf("j: %d -> prev_j:%d\n", j, prev_j) ;
108 printf("i: %d -> prev_i:%d\n", i, prev_i) ;
109 }
110 if (est_state==1) { //was exon labeled "1"
111 est_state = 2 ; //new exon is labeled "2"
112 }
113 else {
114 est_state = 1 ; //last exon was labeled "2", new exon is labeled "1"
115 }
116 for (j; j > prev_j; j--) {
117 if (splice_state == 0) { //coming from exon
118 splice_state = 2 ; //first intron_state: acceptor
119 } else {
120 splice_state = 3 ; //intron
121 }
122 if (j-1 == prev_j) { //last intron_state: donor
123 splice_state = 1 ;//donor
124 }
125 s_align[dna_pos-1] = splice_state;
126 dna_pos-- ;
127 }
128 splice_state = 0 ; //exon again
129
130 }
131
132 startend[1] = i;
133 startend[2] = j;
134 i = prev_i;
135 j = prev_j;
136 act_m_no = prev_m_no ;
137
138 prev_i = ((Pre_score*)matrices[act_m_no] + i*dna_len + j)->prev_i;
139 prev_j = ((Pre_score*)matrices[act_m_no] + i*dna_len + j)->prev_j;
140 prev_m_no = ((Pre_score*)matrices[act_m_no] + i*dna_len +j)->prev_matrix_no;
141 }
142 for (dna_pos; dna_pos>0; dna_pos--) {
143 s_align[dna_pos-1] = 4 ; // not aligned (-1 because array starts with index 0)
144 }
145 for (est_pos; est_pos>0; est_pos--) {
146 e_align[est_pos-1] = 4 ; // not aligned (-1 because array starts with index 0)
147 }
148
149 //printf("(%d, %d) --> (%d, %d)\n", startend[1], startend[2], startend[3], startend[4]) ;
150 //printf("result_length: %d\n", *result_length_ptr) ;
151
152
153 return 0;
154 }
155
156