+ changed interface array_function to array_class
[qpalma.git] / QPalmaDP / qpalma_dp.cpp
1 #include "qpalma_dp.h"
2 #include <cstring>
3 using namespace std;
4
5 /*[splice_align, est_align, weightMatch, alignmentscores, dnaest] = ...
6 myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
7 print_matrix) */
8
9 /* nlhs - Number of expected mxArrays (Left Hand Side)
10 plhs - Array of pointers to expected outputs
11 nrhs - Number of inputs (Right Hand Side)
12 prhs - Array of pointers to input data. */
13
14 Alignment::Alignment() {
15 len = 0;
16 limits = 0;
17 penalties = 0;
18 max_len = 0;
19 min_len = 0;
20 cache = 0;
21 enum ETransformType transform = T_LINEAR;
22 id = 0;
23 name = 0;
24 use_svm = 0;
25 }
26
27 Alignment::~Alignment() {}
28
29 void Alignment::setQualityMatrix(double* qMat, int length){
30 qualityMatrix = new double[length];
31 for(int i=0; i<length; i++)
32 qualityMatrix[i] = qMat[i];
33 }
34 void Alignment::getSpliceAlign(){}
35 void Alignment::getEstAlign(){}
36 void Alignment::getWeightMatch(){}
37 void Alignment::getTotalScores(){}
38 void Alignment::getDNAEST(){}
39
40 void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
41 int est_len_p, struct penalty_struct h, double* matchmatrix, int mm_len,
42 double* donor, int d_len, double* acceptor, int a_len,
43 bool remove_duplicate_scores, bool print_matrix) {
44
45 printf("Entering myalign...\n");
46
47 mlen = 6; // score matrix: length of 6 for "- A C G T N"
48 dna_len = dna_len_p;
49 est_len = est_len_p;
50 nr_paths = nr_paths_p;
51
52 #if 0
53 /***************************************************************************/
54 // 3. h (cell array)
55 /***************************************************************************/
56 mwSize anz_func = 1;
57
58 functions = read_penalty_struct_from_cell(prhs[arg], anz_func);
59 //std::cout << " bla " << functions->max_len << " \n";
60 //std::cout << "lookup_penalty(5) = " << lookup_penalty(functions, 5, NULL, 0)<<std::endl ;
61 arg++; // arg -> 4
62
63
64 /***************************************************************************/
65 // 4. match matrix (columnwise)
66 /***************************************************************************/
67 int nr_entries = mxGetN(prhs[arg]) * mxGetM(prhs[arg]);
68 //mexPrintf("\nNumber of entries in matchmatrix: %d\n", nr_entries) ;
69 double *matchmatrix;
70 matchmatrix = mxGetPr(prhs[arg]);
71 //mexPrintf("matchmatrix\n") ;
72 //for (int i=0; i< nr_entries; i++) {
73 // mexPrintf("entry %d: <%f>\n", i, matchmatrix[i]);
74 //}
75 arg++; // arg -> 5
76 #endif
77
78 /***************************************************************************/
79 // initialize alignment matrices and call fill_matrix()
80 /***************************************************************************/
81
82 //possible donor positions
83 int nr_donor_sites = 0 ;
84 for (int ii=0; ii<d_len; ii++) {
85 if(isfinite(donor[ii])) {
86 nr_donor_sites++ ;
87 }
88 }
89
90 int* donor_sites = new int[nr_donor_sites];
91 int donor_idx = 0 ;
92 for (int ii=0; ii<d_len; ii++) {
93 if(isfinite(donor[ii])) {
94 donor_sites[donor_idx] = ii+1 ;
95 donor_idx++ ;
96 }
97 }
98
99 int* max_score_positions = new int[nr_paths*2];
100
101 Pre_score** matrices = new Pre_score*[nr_paths];
102 for (int z=0; z<nr_paths; z++) {
103 matrices[z] = new Pre_score[dna_len * est_len];
104 }
105
106 printf("calling fill_matrix...\n");
107
108 fill_matrix(nr_paths, matrices, est_len, dna_len, est, dna, &h, matchmatrix, donor, acceptor, remove_duplicate_scores, nr_donor_sites, donor_sites, max_score_positions);
109
110 printf("after call to fill_matrix...\n");
111 /***************************************************************************/
112 // return arguments etc.
113 /***************************************************************************/
114 int result_length; //Eine Variable fuer alle Matrizen
115
116 printf("before init of splice_align and rest...\n");
117
118 splice_align = new int[(dna_len-1)*nr_paths];
119 est_align = new int[(est_len-1)*nr_paths];
120 mmatrix_param = new int[(mlen*mlen)*nr_paths];
121 alignmentscores = new double[nr_paths]; //alignment score for each path/matrix
122
123 printf("Lengths are %d %d %d...\n",dna_len,est_len,mlen);
124 printf("before memset...\n");
125
126 memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
127 memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
128 memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
129 memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
130
131 printf("after memset...\n");
132
133 for (int z=0; z<nr_paths; z++) {
134 result_length = 0 ;
135
136 int* s_align = splice_align + (dna_len-1)*z; //pointer
137 int* e_align = est_align + (est_len-1)*z ; //pointer
138 int* mparam = mmatrix_param + (mlen*mlen)*z; //pointer
139
140 bool no_more_path = result_align(matrices, z, est_len, dna_len, &result_length, est, dna, s_align, e_align, mparam, alignmentscores, max_score_positions);
141
142 // dnaest
143 double *DNA = new double[result_length];
144 double *EST = new double[result_length];
145
146 //backtracking
147 int i = max_score_positions[2*z] ; //i (est)
148 int j = max_score_positions[2*z +1] ; //j (dna)
149 int prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i;
150 int prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
151 int prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
152
153 for (int ii=result_length; ii>0; ii--) {
154 if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
155 DNA[ii-1] = check_char(dna[j-1]) ;
156 EST[ii-1] = check_char(est[i-1]) ;
157 }
158 else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST
159 DNA[ii-1] = check_char(dna[j-1]) ;
160 EST[ii-1] = 0 ; //gap
161 }
162 else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA
163 DNA[ii-1] = 0 ; //gap
164 EST[ii-1] = check_char(est[i-1]) ;
165 }
166 else {// splice site
167 for (j; j > prev_j; j--) {
168 DNA[ii-1] = check_char(dna[j-1]) ;
169 EST[ii-1] = 6 ; //intron
170 ii-- ;
171 }
172 ii++ ; // last ii-- too much (because done in for-loop)
173 }
174
175 i = prev_i;
176 j = prev_j;
177 prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i;
178 prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j;
179 prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
180 }
181
182 } //end of z
183
184 printf("Leaving myalign...\n");
185 }
186
187 void Alignment::getAlignmentResults(int* s_align, int* e_align,
188 int* mmatrix_p, double* alignscores) {
189
190 printf("Entering getAlignmentResults...\n");
191 uint splice_align_size = (dna_len-1)*nr_paths;
192 uint est_align_size = (est_len-1)*nr_paths;
193 uint mmatrix_param_size = (mlen*mlen)*nr_paths;
194 uint alignmentscores_size = nr_paths; //alignment score for each path/matrix
195
196 for(int idx=0; idx<splice_align_size; idx++)
197 s_align[idx] = splice_align[idx];
198
199 for(int idx=0; idx<est_align_size; idx++)
200 e_align[idx] = est_align[idx];
201
202 for(int idx=0; idx<mmatrix_param_size; idx++)
203 mmatrix_p[idx] = mmatrix_param[idx];
204
205 for(int idx=0; idx<alignmentscores_size; idx++)
206 alignscores[idx] = alignmentscores[idx];
207
208 printf("Leaving getAlignmentResults...\n");
209 }
210