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