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