QPalmaDP/qpalma_dp.cpp

   1 #include "qpalma_dp.h"
   2
   3 /*[splice_align, est_align, weightMatch, alignmentscores, dnaest] = ...
   4   myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
   5   print_matrix) */
   6
   7 /* nlhs - Number of expected mxArrays (Left Hand Side)
   8    plhs - Array of pointers to expected outputs
   9    nrhs - Number of inputs (Right Hand Side)
  10    prhs - Array of pointers to input data. */
  11
  12   myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
  13   print_matrix) */
  14
  15 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
  16 {
  17
  18
  19
  20   /***************************************************************************/
  21   // right number of arguments?
  22   /***************************************************************************/
  23   if(nrhs != 9)
  24   {
  25     mexPrintf("10 input arguments please: \n [splice_align, est_align, weightMatch, alignmentscores, dnaest] = myalign(nr_paths, dna, est, cellarray, matchmatrix, donor, acceptor, remove_duplicate_scores, print_matrix)\n");
  26     return ;
  27   }
  28
  29   if(nlhs != 5)
  30   {
  31     mexPrintf("5 output arguments please\n [splice_align, est_align, weightMatch, alignmentscores, dnaest] = myalign(nr_paths, dna, est, cellarray, matchmatrix, donor, acceptor, remove_duplicate_scores, print_matrix)\n");
  32     return ;
  33   }
  34
  35
  36   /***************************************************************************/
  37   // variables
  38   /***************************************************************************/
  39   const int mlen = 6; // score matrix: length of 6 for "- A C G T N"
  40   int arg = 0 ; // counter variable
  41   int failure ; //
  42
  43
  44   /***************************************************************************/
  45   // 0. nr_paths
  46   /***************************************************************************/
  47   int nr_paths = (int)*mxGetPr(prhs[arg]); // 0. input (double)
  48   arg++; // arg -> 1
  49
  50
  51   /***************************************************************************/
  52   // 1. dna
  53   /***************************************************************************/
  54   int dna_len = mxGetN(prhs[arg]) + 1; // mxGetN: number of columns (int)
  55   char* dna = (char*)mxCalloc(dna_len, sizeof(char));
  56
  57   failure = mxGetString(prhs[arg], dna, dna_len); // NULL at end of string
  58   if (failure) { // mxGetString returns 0 if success, 1 on failure
  59     mexPrintf("couldn't open dna");
  60     return;
  61   }
  62   arg++; // arg -> 2
  63
  64
  65   /***************************************************************************/
  66   // 2. est
  67   /***************************************************************************/
  68   int est_len = mxGetN(prhs[arg]) + 1;
  69   char* est = (char*)mxCalloc(est_len, sizeof(char));
  70
  71   failure = mxGetString(prhs[arg], est, est_len);
  72   if (failure) {
  73     mexPrintf("couldn't open est");
  74     return;
  75   }
  76   arg++; // arg -> 3
  77
  78
  79   /***************************************************************************/
  80   // 3. h (cell array)
  81   /***************************************************************************/
  82   penalty_struct* functions;
  83   mwSize anz_func = 1;
  84
  85   functions = read_penalty_struct_from_cell(prhs[arg], anz_func);
  86   //std::cout << " bla " << functions->max_len << " \n";
  87   //std::cout << "lookup_penalty(5) = " << lookup_penalty(functions, 5, NULL, 0)<<std::endl ;
  88   arg++; // arg -> 4
  89
  90
  91   /***************************************************************************/
  92   // 4. match matrix (columnwise)
  93   /***************************************************************************/
  94   int nr_entries = mxGetN(prhs[arg]) * mxGetM(prhs[arg]);
  95   //mexPrintf("\nNumber of entries in matchmatrix: %d\n", nr_entries) ;
  96   double *matchmatrix;
  97   matchmatrix = mxGetPr(prhs[arg]);
  98   //mexPrintf("matchmatrix\n") ;
  99   //for (int i=0; i< nr_entries; i++) {
 100   //  mexPrintf("entry %d: <%f>\n", i, matchmatrix[i]);
 101   //}
 102   arg++; // arg -> 5
 103
 104
 105   /***************************************************************************/
 106   // 5. donor
 107   /***************************************************************************/
 108   int donor_length = mxGetN(prhs[arg]);
 109   double *donor;
 110   donor = mxGetPr(prhs[arg]);
 111   //mexPrintf("Donor eingelesen: <%f>\n", donor[0]);
 112   arg++; // arg -> 6
 113
 114
 115   /***************************************************************************/
 116   // 6. acceptor
 117   /***************************************************************************/
 118   int acceptor_length = mxGetN(prhs[arg]);
 119   double *acceptor;
 120   acceptor = mxGetPr(prhs[arg]);
 121   //mexPrintf("Acceptor eingelesen: <%f>\n", acceptor[0]);
 122   arg++; // arg -> 7
 123
 124
 125   /***************************************************************************/
 126   // 7.remove_duplicate_scores
 127   /***************************************************************************/
 128   bool remove_duplicate_scores = (bool)*mxGetPr(prhs[arg]);
 129   arg++; // arg -> 8
 130
 131
 132   /***************************************************************************/
 133   // 8. print_matrix
 134   /***************************************************************************/
 135   int print_matrix = (int)*mxGetPr(prhs[arg]);
 136   arg++; // arg -> 9
 137
 138
 139   /***************************************************************************/
 140   // initialize alignment matrices and call fill_matrix()
 141   /***************************************************************************/
 142
 143   //possible donor positions
 144   int nr_donor_sites = 0 ;
 145   for (int ii=0; ii<donor_length; ii++) {
 146     if(isfinite(donor[ii])) {
 147       nr_donor_sites++ ;
 148     }
 149   }
 150   //mwSize donor_sites[nr_donor_sites] ;
 151   mwSize* donor_sites = (mwSize*)mxCalloc(nr_donor_sites,sizeof(mwSize));
 152   int donor_idx = 0 ;
 153   for (int ii=0; ii<donor_length; ii++) {
 154     if(isfinite(donor[ii])) {
 155       donor_sites[donor_idx] = ii+1 ;
 156       donor_idx++ ;
 157     }
 158   }
 159
 160
 161   //int max_score_positions[nr_paths*2]; //local alignment: where is the highest alignment score [i1,j1,i2,j2...]
 162   mwSize* max_score_positions = (mwSize*)mxCalloc(nr_paths*2,sizeof(mwSize));
 163
 164   //Pre_score* matrices[nr_paths];
 165   Pre_score** matrices = (Pre_score**)mxCalloc(nr_paths,sizeof(Pre_score*));
 166   for (int z=0; z<nr_paths; z++) {
 167       //matrices[z] = new Pre_score[dna_len * est_len];
 168       //mexPrintf("size of one matrice: %d\n", est_len*dna_len*sizeof(Pre_score)) ;
 169       matrices[z] = (Pre_score*)mxCalloc(dna_len*est_len,sizeof(Pre_score));
 170   }
 171
 172   fill_matrix(nr_paths, matrices, est_len, dna_len, est, dna, functions, matchmatrix, donor, acceptor, remove_duplicate_scores, nr_donor_sites, donor_sites, max_score_positions);
 173
 174   //just some info
 175   //for (int z=0; z<nr_paths; z++) {
 176   //  mexPrintf("(%d, %d)\n", max_score_positions[2*z], max_score_positions[2*z +1]) ;
 177   //}
 178
 179   /***************************************************************************/
 180   // return arguments etc.
 181   /***************************************************************************/
 182   int result_length; //Eine Variable fuer alle Matrizen
 183
 184   //int splice_align[(dna_len-1)*nr_paths]; //SpliceAlign for DNA - 1 line per matrix
 185   //int est_align[(est_len-1)*nr_paths]; //SpliceAlign for EST - 1 line per matrix
 186   //int mmatrix_param[(mlen*mlen)*nr_paths];//Matrix fuer Scorematrixparameter - 1 Zeile pro Matrix
 187   int* splice_align = (int*)mxCalloc((dna_len-1)*nr_paths,sizeof(int));
 188   int* est_align = (int*)mxCalloc((est_len-1)*nr_paths,sizeof(int));
 189   int* mmatrix_param = (int*)mxCalloc((mlen*mlen)*nr_paths,sizeof(int));
 190   double alignmentscores[nr_paths]; //alignment score for each path/matrix
 191
 192
 193   memset((char*)splice_align, -1, (dna_len-1)*nr_paths*sizeof(int)); // fills splice_align with zeros
 194   memset((char*)est_align, -1, (est_len-1)*nr_paths*sizeof(int)); // fills est_align with zeros
 195   memset((char*)mmatrix_param, 0, mlen*mlen*nr_paths*sizeof(int)); //fills mmatrix_param with zeros
 196   memset(alignmentscores, -1, nr_paths*sizeof(double)); //fills alignmentscores with zeros
 197
 198   // dnaest
 199   plhs[nlhs-1] = mxCreateCellMatrix(2, nr_paths);
 200
 201
 202
 203   //mexPrintf("DEBUG: decoding %d paths\n",nr_paths);
 204   for (int z=0; z<nr_paths; z++) {
 205     result_length = 0 ;
 206
 207     int* s_align = splice_align + (dna_len-1)*z;  //pointer
 208     int* e_align = est_align + (est_len-1)*z ;    //pointer
 209     int* mparam  = mmatrix_param + (mlen*mlen)*z; //pointer
 210
 211     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);
 212
 213     // dnaest
 214     mxArray* mxdna = mxCreateDoubleMatrix(result_length, 1, mxREAL);
 215     mxArray* mxest = mxCreateDoubleMatrix(result_length, 1, mxREAL);
 216     double *DNA = mxGetPr(mxdna) ;
 217     double *EST = mxGetPr(mxest) ;
 218     mxSetCell(plhs[nlhs-1], z*2, mxdna) ;
 219     mxSetCell(plhs[nlhs-1], z*2+1, mxest) ;
 220
 221     //backtracking
 222     mwSize i = max_score_positions[2*z] ; //i (est)
 223     mwSize j = max_score_positions[2*z +1] ; //j (dna)
 224     mwSize prev_i = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_i;
 225     mwSize prev_j = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_j;
 226     mwSize prev_m_no = ((Pre_score*)matrices[z] + i*dna_len +j)->prev_matrix_no;
 227
 228     for (int ii=result_length; ii>0; ii--) {
 229       if ((prev_i == (i-1)) && (prev_j == (j-1))) { // match or mismatch
 230         DNA[ii-1] = check_char(dna[j-1]) ;
 231         EST[ii-1] = check_char(est[i-1]) ;
 232       }
 233       else if ((prev_i == (i)) && (prev_j == (j-1))) {// gap on EST
 234         DNA[ii-1] = check_char(dna[j-1]) ;
 235         EST[ii-1] = 0 ; //gap
 236       }
 237       else if ((prev_i == (i-1)) && (prev_j == (j))) { // gap on DNA
 238         DNA[ii-1] = 0 ; //gap
 239         EST[ii-1] = check_char(est[i-1]) ;
 240       }
 241       else {// splice site
 242         for (j; j > prev_j; j--) {
 243           DNA[ii-1] = check_char(dna[j-1]) ;
 244           EST[ii-1] = 6 ; //intron
 245           ii-- ;
 246         }
 247         ii++ ; // last ii-- too much (because done in for-loop)
 248       }
 249
 250       i = prev_i;
 251       j = prev_j;
 252       prev_i = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_i;
 253       prev_j = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_j;
 254       prev_m_no = ((Pre_score*)matrices[prev_m_no] + i*dna_len + j)->prev_matrix_no;
 255     }
 256
 257   } //end of z
 258
 259
 260
 261   //mexPrintf("DEBUG:Returning arguments\n");
 262   int index = 0 ;
 263
 264   //splice_align (mxArray):
 265   const mwSize dims0[] = {(dna_len-1), nr_paths};
 266   //const int dims0[] = {(dna_len-1), nr_paths};
 267   plhs[index] = mxCreateNumericArray(2, dims0, mxINT32_CLASS, mxREAL);
 268   mwSize* start0 = (mwSize*)(mxGetData(plhs[index]));
 269   memcpy(start0, splice_align, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
 270   index++ ; // ->1
 271
 272   //est_align (mxArray):
 273   const mwSize dims1[] = {(est_len-1), nr_paths};
 274   //const int dims1[] = {(est_len-1), nr_paths};
 275   plhs[index] = mxCreateNumericArray(2, dims1, mxINT32_CLASS, mxREAL);
 276   mwSize* start1 = (mwSize*)(mxGetData(plhs[index]));
 277   memcpy(start1, est_align, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
 278   index++ ; // ->2
 279
 280   //weightmatch (mxArray)
 281   const mwSize dims2[] = {mlen*mlen, nr_paths};
 282   //const int dims2[] = {mlen*mlen, nr_paths};
 283   plhs[index] = mxCreateNumericArray(2, dims2, mxINT32_CLASS, mxREAL);
 284   mwSize* start2= (mwSize*)(mxGetData(plhs[index]));
 285   memcpy(start2, mmatrix_param, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
 286   index++ ; // ->3
 287
 288   //alignmentscore
 289   plhs[index] = mxCreateDoubleMatrix(1,nr_paths,mxREAL);
 290   if (plhs[index] == NULL)
 291   {
 292     printf("%s : Out of memory on line %d\n", __FILE__, __LINE__);
 293     printf("Unable to create mxArray.\n");
 294     exit(1);
 295   }
 296   double* start3 = (double*)(mxGetData(plhs[index]));
 297   memcpy(start3, alignmentscores, mxGetM(plhs[index])*mxGetN(plhs[index])*mxGetElementSize(plhs[index]));
 298   index++ ; // ->4
 299
 300   assert(nlhs==index+1) ;
 301
 302   //Free Memory
 303   //mexPrintf("DEBUG:Free Memory\n");
 304   mxFree(donor_sites);
 305   mxFree(max_score_positions);
 306   for (int i=nr_paths-1;i>=0; i--)
 307   {
 308       //delete[] matrices[i];
 309       mxFree(matrices[i]);
 310   }
 311   mxFree(matrices);
 312   mxFree(splice_align);
 313   mxFree(est_align);
 314   mxFree(mmatrix_param);
 315   mxFree(dna);
 316   mxFree(est);
 317
 318   return;
 319
 320 }
 321