+ added some assertions
[qpalma.git] / QPalmaDP / qpalma_dp.cpp
1 #include "qpalma_dp.h"
2 #include "debug_tools.h"
3 #include <cstring>
4 using namespace std;
5
6 /*[splice_align, est_align, weightMatch, alignmentscores, dnaest] = ...
7 myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
8 print_matrix) */
9
10 Alignment::Alignment(int numQPlifs, int numq) {
11 len = 0;
12 limits = 0;
13 penalties = 0;
14 max_len = 0;
15 min_len = 0;
16 cache = 0;
17 enum ETransformType transform = T_LINEAR;
18 id = 0;
19 name = 0;
20 use_svm = 0;
21
22 // set ptrs to zero first
23 splice_align = 0;
24 est_align = 0;
25 mmatrix_param = 0;
26 alignmentscores = 0;
27 qualityScoresAllPaths = 0;
28
29 numQualSuppPoints = numq;
30 totalNumPlifs = numQPlifs;
31
32 FA( numQualSuppPoints >= 0 );
33 FA( totalNumPlifs >= 0 );
34 }
35
36 void Alignment::getDNAEST(){}
37
38 void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
39 int est_len_p, double* prb, double* chastity, struct penalty_struct h, double* matchmatrix, int mm_len,
40 double* donor, int d_len, double* acceptor, int a_len, struct penalty_struct* qualityScores,
41 bool remove_duplicate_scores, bool print_matrix, bool use_quality_scores) {
42
43 // printf("Entering myalign...\n");
44 mlen = 6; // score matrix: length of 6 for "- A C G T N"
45 dna_len = dna_len_p + 1 ;
46 est_len = est_len_p + 1 ;
47 nr_paths = nr_paths_p;
48
49 /***************************************************************************/
50 // initialize alignment matrices and call fill_matrix()
51 /***************************************************************************/
52
53 mode currentMode;
54 if (use_quality_scores)
55 currentMode = USE_QUALITY_SCORES;
56 else
57 currentMode = NORMAL;
58
59 // dnaest
60 double *DNA_ARRAY = 0;
61 double *EST_ARRAY = 0;
62
63 //possible donor positions
64 int nr_donor_sites = 0 ;
65 for (int ii=0; ii<d_len; ii++) {
66 if(isfinite(donor[ii])) {
67 nr_donor_sites++ ;
68 }
69 }
70
71 int* donor_sites = new int[nr_donor_sites];
72 int donor_idx = 0 ;
73 for (int ii=0; ii<d_len; ii++) {
74 if(isfinite(donor[ii])) {
75 donor_sites[donor_idx] = ii+1 ;
76 donor_idx++ ;
77 }
78 }
79
80
81 int* max_score_positions = new int[nr_paths*2];
82
83 Pre_score** matrices = new Pre_score*[nr_paths];
84 for (int z=0; z<nr_paths; z++) {
85 matrices[z] = new Pre_score[dna_len * est_len];
86 }
87
88 // printf("calling fill_matrix...\n");
89 //
90
91 fill_matrix(nr_paths, matrices, est_len, dna_len, est, dna, prb, &h, matchmatrix, qualityScores, donor, acceptor, remove_duplicate_scores, nr_donor_sites, donor_sites, max_score_positions,currentMode);
92
93 //printf("after call to fill_matrix...\n");
94 /***************************************************************************/
95 // return arguments etc.
96 /***************************************************************************/
97 int result_length; //Eine Variable fuer alle Matrizen
98
99 // printf("before init of splice_align and rest...\n");
100 //
101 splice_align_size = (dna_len-1)*nr_paths;
102 est_align_size = (est_len-1)*nr_paths;
103 mmatrix_param_size = (mlen*mlen)*nr_paths;
104 alignmentscores_size = nr_paths; //alignment score for each path/matrix
105 qScores_size = totalNumPlifs*nr_paths; //alignment score for each path/matrix
106
107 //printf("dna_len is %d est_len is %d mmatrix_len is %d",splice_align_size, est_align_size, mmatrix_param_size);
108
109 splice_align = new int[splice_align_size];
110 est_align = new int[est_align_size];
111 mmatrix_param = new int[mmatrix_param_size];
112 alignmentscores = new double[nr_paths]; //alignment score for each path/matrix
113 qualityScoresAllPaths = new penalty_struct[qScores_size];
114 int qidx, pidx;
115
116 for(qidx=0;qidx<totalNumPlifs*nr_paths;qidx++) {
117 penalty_struct p;
118 init_penalty_struct(p);
119 p.len = numQualSuppPoints;
120 p.limits = (REAL*) calloc(p.len,sizeof(REAL));
121 for(pidx=0;pidx<p.len;pidx++)
122 p.limits[pidx] = qualityScores[qidx%totalNumPlifs].limits[pidx];
123
124 p.penalties = (REAL*) calloc(p.len,sizeof(REAL));
125
126 //for(pidx=0;pidx<p.len;pidx++)
127 // printf("%f ",p.limits[pidx]);
128 //printf("\n");
129
130 //for(pidx=0;pidx<p.len;pidx++)
131 // printf("%f ",p.penalties[pidx]);
132 //printf("\n");
133
134 qualityScoresAllPaths[qidx] = p;
135 }
136
137 // printf("before memset...\n");
138 memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
139 memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
140 memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
141 memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
142 //printf("after memset...\n");
143
144
145 for (int z=0; z<nr_paths; z++) {
146 result_length = 0 ;
147
148 int* s_align = splice_align + (dna_len-1)*z; //pointer
149 int* e_align = est_align + (est_len-1)*z ; //pointer
150 int* mparam = mmatrix_param + (mlen*mlen)*z; //pointer
151 penalty_struct* qparam = qualityScoresAllPaths + (totalNumPlifs*z);
152
153 //printf("before call to result_align...\n");
154 bool no_more_path = result_align(matrices, z, est_len, dna_len, &result_length, est, dna, prb, chastity, s_align, e_align, mparam, alignmentscores, max_score_positions, qparam, currentMode);
155 //printf("after call to result_align...\n");
156
157 //for(qidx=0;qidx<totalNumPlifs*nr_paths;qidx++) {
158 // penalty_struct p;
159 // p = qualityScoresAllPaths[qidx];
160 // for(pidx=0;pidx<p.len;pidx++)
161 // printf("%f ",p.penalties[pidx]);
162 // printf("\n");
163 //}
164
165 //if(DNA_ARRAY != 0) {
166 // delete[] DNA_ARRAY;
167 // delete[] EST_ARRAY;
168 // }
169
170 //DNA_ARRAY = new double[result_length];
171 //EST_ARRAY = new double[result_length];
172
173 // //backtracking
174 // int i = max_score_positions[2*z] ; //i (est)
175 // int j = max_score_positions[2*z +1] ; //j (dna)
176 // int prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i;
177 // int prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
178 // int prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
179 //
180 // for (int ii=result_length; ii>0; ii--) {
181 // if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
182 //DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
183 //EST_ARRAY[ii-1] = check_char(est[i-1]) ;
184 // }
185 // else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST_ARRAY
186 //DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
187 //EST_ARRAY[ii-1] = 0 ; //gap
188 // }
189 // else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA_ARRAY
190 //DNA_ARRAY[ii-1] = 0 ; //gap
191 //EST_ARRAY[ii-1] = check_char(est[i-1]) ;
192 // }
193 // else {// splice site
194 //for (j; j > prev_j; j--) {
195 // DNA_ARRAY[ii-1] = check_char(dna[j-1]) ;
196 // EST_ARRAY[ii-1] = 6 ; //intron
197 // ii-- ;
198 //}
199 //ii++ ; // last ii-- too much (because done in for-loop)
200 // }
201 //
202 // i = prev_i;
203 // j = prev_j;
204 // prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i;
205 // prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j;
206 // prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
207 // }
208 //
209 } //end of z
210
211 //if(DNA_ARRAY != 0) {
212 // delete[] DNA_ARRAY;
213 // delete[] EST_ARRAY;
214 // }
215
216 for (int i=nr_paths-1;i>=0; i--)
217 delete[] matrices[i];
218
219 //printf("Leaving myalign...\n");
220 }
221
222 void Alignment::getAlignmentResults(int* s_align, int* e_align,
223 int* mmatrix_p, double* alignscores, double* qScores) {
224
225 //printf("Entering getAlignmentResults...\n");
226
227 int idx;
228 for(idx=0; idx<splice_align_size; idx++)
229 s_align[idx] = splice_align[idx];
230
231 for(idx=0; idx<est_align_size; idx++)
232 e_align[idx] = est_align[idx];
233
234 for(idx=0; idx<mmatrix_param_size; idx++)
235 mmatrix_p[idx] = mmatrix_param[idx];
236
237 for(idx=0; idx<alignmentscores_size; idx++)
238 alignscores[idx] = alignmentscores[idx];
239
240 penalty_struct currentPlif;
241 int pidx;
242 int ctr=0;
243 for(idx=0; idx<qScores_size; idx++) {
244 currentPlif = qualityScoresAllPaths[idx];
245 //printf("Size is %dx%d\n",qScores_size,currentPlif.len);
246 for(pidx=0; pidx<currentPlif.len; pidx++) {
247 qScores[ctr] = currentPlif.penalties[pidx];
248 ctr++;
249 }
250 }
251
252 //printf("ctr is %d\n",ctr);
253 //printf("Leaving getAlignmentResults...\n");
254 }