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