+ fixed some index bugs
[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 Alignment::Alignment() {
10 len = 0;
11 limits = 0;
12 penalties = 0;
13 max_len = 0;
14 min_len = 0;
15 cache = 0;
16 enum ETransformType transform = T_LINEAR;
17 id = 0;
18 name = 0;
19 use_svm = 0;
20 }
21
22 Alignment::~Alignment() {}
23
24 void Alignment::setQualityMatrix(double* qMat, int length){
25 qualityMatrix = new double[length];
26 for(int i=0; i<length; i++)
27 qualityMatrix[i] = qMat[i];
28 }
29 void Alignment::getSpliceAlign(){}
30 void Alignment::getEstAlign(){}
31 void Alignment::getWeightMatch(){}
32 void Alignment::getTotalScores(){}
33 void Alignment::getDNAEST(){}
34
35 void Alignment::myalign(int nr_paths_p, char* dna, int dna_len_p, char* est,
36 int est_len_p, struct penalty_struct h, double* matchmatrix, int mm_len,
37 double* donor, int d_len, double* acceptor, int a_len,
38 bool remove_duplicate_scores, bool print_matrix) {
39
40 printf("Entering myalign...\n");
41
42 mlen = 6; // score matrix: length of 6 for "- A C G T N"
43 dna_len = dna_len_p + 1 ;
44 est_len = est_len_p + 1 ;
45 nr_paths = nr_paths_p;
46
47 /***************************************************************************/
48 // initialize alignment matrices and call fill_matrix()
49 /***************************************************************************/
50
51 //possible donor positions
52 int nr_donor_sites = 0 ;
53 for (int ii=0; ii<d_len; ii++) {
54 if(isfinite(donor[ii])) {
55 nr_donor_sites++ ;
56 }
57 }
58
59 int* donor_sites = new int[nr_donor_sites];
60 int donor_idx = 0 ;
61 for (int ii=0; ii<d_len; ii++) {
62 if(isfinite(donor[ii])) {
63 donor_sites[donor_idx] = ii+1 ;
64 donor_idx++ ;
65 }
66 }
67
68 int* max_score_positions = new int[nr_paths*2];
69
70 Pre_score** matrices = new Pre_score*[nr_paths];
71 for (int z=0; z<nr_paths; z++) {
72 matrices[z] = new Pre_score[dna_len * est_len];
73 }
74
75 printf("calling fill_matrix...\n");
76
77 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);
78
79 printf("after call to fill_matrix...\n");
80 /***************************************************************************/
81 // return arguments etc.
82 /***************************************************************************/
83 int result_length; //Eine Variable fuer alle Matrizen
84
85 printf("before init of splice_align and rest...\n");
86
87 splice_align = new int[(dna_len-1)*nr_paths];
88 est_align = new int[(est_len-1)*nr_paths];
89 mmatrix_param = new int[(mlen*mlen)*nr_paths];
90 alignmentscores = new double[nr_paths]; //alignment score for each path/matrix
91
92 printf("before memset...\n");
93
94 memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
95 memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
96 memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
97 memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
98
99 printf("after memset...\n");
100
101 for (int z=0; z<nr_paths; z++) {
102 result_length = 0 ;
103
104 int* s_align = splice_align + (dna_len-1)*z; //pointer
105 int* e_align = est_align + (est_len-1)*z ; //pointer
106 int* mparam = mmatrix_param + (mlen*mlen)*z; //pointer
107
108 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);
109
110 // dnaest
111 double *DNA = new double[result_length];
112 double *EST = new double[result_length];
113
114 //backtracking
115 int i = max_score_positions[2*z] ; //i (est)
116 int j = max_score_positions[2*z +1] ; //j (dna)
117 int prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i;
118 int prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
119 int prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
120
121 for (int ii=result_length; ii>0; ii--) {
122 if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
123 DNA[ii-1] = check_char(dna[j-1]) ;
124 EST[ii-1] = check_char(est[i-1]) ;
125 }
126 else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST
127 DNA[ii-1] = check_char(dna[j-1]) ;
128 EST[ii-1] = 0 ; //gap
129 }
130 else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA
131 DNA[ii-1] = 0 ; //gap
132 EST[ii-1] = check_char(est[i-1]) ;
133 }
134 else {// splice site
135 for (j; j > prev_j; j--) {
136 DNA[ii-1] = check_char(dna[j-1]) ;
137 EST[ii-1] = 6 ; //intron
138 ii-- ;
139 }
140 ii++ ; // last ii-- too much (because done in for-loop)
141 }
142
143 i = prev_i;
144 j = prev_j;
145 prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i;
146 prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j;
147 prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
148 }
149
150 } //end of z
151
152 printf("Leaving myalign...\n");
153 }
154
155 void Alignment::getAlignmentResults(int* s_align, int* e_align,
156 int* mmatrix_p, double* alignscores) {
157
158 printf("Entering getAlignmentResults...\n");
159 uint splice_align_size = (dna_len-1)*nr_paths;
160 uint est_align_size = (est_len-1)*nr_paths;
161 uint mmatrix_param_size = (mlen*mlen)*nr_paths;
162 uint alignmentscores_size = nr_paths; //alignment score for each path/matrix
163
164 for(int idx=0; idx<splice_align_size; idx++)
165 s_align[idx] = splice_align[idx];
166
167 for(int idx=0; idx<est_align_size; idx++)
168 e_align[idx] = est_align[idx];
169
170 for(int idx=0; idx<mmatrix_param_size; idx++)
171 mmatrix_p[idx] = mmatrix_param[idx];
172
173 for(int idx=0; idx<alignmentscores_size; idx++)
174 alignscores[idx] = alignmentscores[idx];
175
176 printf("Leaving getAlignmentResults...\n");
177 }
178