+ started to add quality scoring function to C code

TODO
- think about apropriate replacement for global vars
- think about how dnaest should be visible to the python layer

git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@7100 e1793c9e-67f9-0310-80fc-b846ff1f7b36

diff --git a/QPalmaDP/fill_matrix.cpp b/QPalmaDP/fill_matrix.cpp
index c7a32fe..9dfebc1 100644 (file)
--- a/QPalmaDP/fill_matrix.cpp
+++ b/QPalmaDP/fill_matrix.cpp
@@ -67,23 +67,19 @@ int check_char(char base)
  * system of Palma. Otherwiese scoring_functions should point to an array of
  * plifs scoring the respective pairs of characters together with the EST
  * quality score.
- *
- *
  */
 
-double getScore(double *matchmatrix,int mlen, int dnaChar, int estChar, struct penalty_struct *scoring_functions) {
-   double score = 0.0;
+double scoreMatch(int dnaChar, int estChar, double baseScore){
+   for(int idx=0;idx<24;idx++)
+      struct penalty_struct currentPen = ::qualityScores[idx];
 
-   if( scoring_functions == 0 ) {
-      if( estChar == -1 )
-         score = matchmatrix[mlen*dnaChar];
+   return 1.0;
+}
 
-      if( dnaChar == -1 )
-         score = matchmatrix[estChar];
+double getScore(double *matchmatrix,int mlen, int dnaChar, int estChar, double baseScore, mode scoringScheme) {
+   double score = 0.0;
 
-      if( estChar != -1 && dnaChar != -1 )
-         score = matchmatrix[mlen* dnaChar +estChar];
-   } else {
+   if( scoringScheme == NORMAL ) {
       if( estChar == -1 )
          score = matchmatrix[mlen*dnaChar];
 
@@ -91,8 +87,9 @@ double getScore(double *matchmatrix,int mlen, int dnaChar, int estChar, struct p
          score = matchmatrix[estChar];
 
       if( estChar != -1 && dnaChar != -1 )
-         //scoring_functions
          score = matchmatrix[mlen* dnaChar +estChar];
+   } else {
+         score = scoreMatch(dnaChar,estChar,baseScore);
    }
 
    return score;
@@ -108,6 +105,16 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
   /*********************************************************************************************/
   const int mlen = 6; // length of matchmatrix
   int estChar, dnaChar; // 0:'-', 1:'a', 2:'c', 3:'g', 4:'t', 5:'n'
+
+
+  double baseScore;
+  double *est_scores = new double[est_len];
+  for(int i=0;i<est_len;i++)
+   est_scores[i] = 0.0;
+
+   mode currentMode = NORMAL;
+
+
   double prevValue ;
   double tempValue;
   int max_len = (int)functions->limits[functions->len-1] ; //last (len) supporting point of h
@@ -216,6 +223,7 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
     {
       dnaChar = check_char(dna[j-1]) ; // -1 because of of the gaps in the 0.line/column of the al-matrix
       estChar = check_char(est[i-1]) ;
+      baseScore = est_scores[i-1];
 
       
       //if introns of length greater than max_len possible
@@ -295,7 +303,7 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
        if (isfinite(prevValue))
        {
          // tempValue = prevValue +(matchmatrix[mlen* dnaChar]);   /* score(gap,DNA) */
-         tempValue = prevValue + getScore(matchmatrix,mlen,dnaChar,-1,scoring_functions);
+         tempValue = prevValue + getScore(matchmatrix,mlen,dnaChar,-1,baseScore,currentMode);
          if (isfinite(tempValue))
          {
            assert(num_elem<max_num_elem) ;
@@ -316,7 +324,7 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
        if (isfinite(prevValue))
        {
          //tempValue = prevValue +(matchmatrix[mlen* dnaChar +estChar]);//DNA mit EST
-         tempValue = prevValue + getScore(matchmatrix,mlen,dnaChar,estChar,scoring_functions);
+         tempValue = prevValue + getScore(matchmatrix,mlen,dnaChar,estChar,baseScore,currentMode);
          if (isfinite(tempValue))
          {
            assert(num_elem<max_num_elem) ;
@@ -337,7 +345,7 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
        if (isfinite(prevValue))
        {
          // tempValue = prevValue +(matchmatrix[estChar]); /* score(EST,gap) */
-         tempValue = prevValue + getScore(matchmatrix,mlen,-1,estChar,scoring_functions);
+         tempValue = prevValue + getScore(matchmatrix,mlen,-1,estChar,baseScore,currentMode);
          if (isfinite(tempValue))
          {
            assert(num_elem<max_num_elem) ;

diff --git a/QPalmaDP/fill_matrix.h b/QPalmaDP/fill_matrix.h
index d10c233..4c3a10e 100644 (file)
--- a/QPalmaDP/fill_matrix.h
+++ b/QPalmaDP/fill_matrix.h
@@ -50,4 +50,9 @@ extern void print_align(Pre_score* matrix, int length_est,  int length_dna, Alig
 
 double getScore(double *matchmatrix,int mlen, int dnaChar, int estChar);
 
+enum mode { NORMAL, USE_QUALITY_SCORES};
+
+double scoreMatch(int dnaChar, int estChar, double baseScore);
+
+
 #endif // _FILLMATRIX_H__

diff --git a/QPalmaDP/penalty_info.cpp b/QPalmaDP/penalty_info.cpp
index d3f11a2..d5d739f 100644 (file)
--- a/QPalmaDP/penalty_info.cpp
+++ b/QPalmaDP/penalty_info.cpp
@@ -58,6 +58,20 @@ void delete_penalty_struct_array(struct penalty_struct *PEN, INT len)
        delete[] PEN ;
 }
 
+void copy_penalty_struct(struct penalty_struct *old, struct penalty_struct *newp) {
+   newp->len       = old->len;
+   newp->limits    = old->limits;
+   newp->penalties = old->penalties;
+   newp->max_len   = old->max_len;
+   newp->min_len   = old->min_len;
+   newp->cache     = old->cache;
+   newp->transform = old->transform;
+   newp->id        = old->id;
+   newp->next_pen  = old->next_pen;
+   newp->name      = old->name;
+   newp->use_svm   = old->use_svm;
+}
+
 
 //struct penalty_struct * read_penalty_struct_from_cell(const mxArray * mx_penalty_info, int &P)
 //{

diff --git a/QPalmaDP/penalty_info.h b/QPalmaDP/penalty_info.h
index 56465ca..44ce621 100644 (file)
--- a/QPalmaDP/penalty_info.h
+++ b/QPalmaDP/penalty_info.h
@@ -31,6 +31,7 @@ struct penalty_struct
 void init_penalty_struct(struct penalty_struct &PEN) ;
 void delete_penalty_struct_palma(struct penalty_struct &PEN) ;
 //void delete_penalty_struct_array(struct penalty_struct *PEN, INT len) ;
+void copy_penalty_struct(struct penalty_struct *old, struct penalty_struct *newp);
 
 //#ifdef HAVE_MATLAB
 //struct penalty_struct * read_penalty_struct_from_cell(const mxArray * mx_penalty_info, mwSize &P) ;
@@ -39,4 +40,6 @@ void delete_penalty_struct_palma(struct penalty_struct &PEN) ;
 REAL lookup_penalty(const struct penalty_struct *PEN, INT p_value, 
                                        REAL* svm_values, bool follow_next=true) ;
 
+extern struct penalty_struct *qualityScores;
+
 #endif

diff --git a/QPalmaDP/qpalma_dp.cpp b/QPalmaDP/qpalma_dp.cpp
index 5197a25..4bd4f44 100644 (file)
--- a/QPalmaDP/qpalma_dp.cpp
+++ b/QPalmaDP/qpalma_dp.cpp
@@ -6,6 +6,8 @@ using namespace std;
   myalign([nr_paths], dna, est, {h}, matchmatrix, donor, acceptor, remove_duplicate_scores, ...
   print_matrix) */
 
+struct penalty_struct *qualityScores = 0;
+
 Alignment::Alignment() {
       len = 0;
       limits = 0;
@@ -23,7 +25,8 @@ Alignment::Alignment() {
       mmatrix_param = 0;
       alignmentscores = 0;
       qualityMatrix = 0;
-      matching_plifs = 0;
+
+      ::qualityScores = new struct penalty_struct[24];
 }
 
 void Alignment::setQualityMatrix(double* qMat, int length){
@@ -34,11 +37,11 @@ void Alignment::setQualityMatrix(double* qMat, int length){
       qualityMatrix[i] = qMat[i];
 }
 
-void Alignment::setMatchPlifs(struct penalty_struct* match_plifs) {
-   matching_plifs = new struct penalty_struct[10];
 
-   for(int i=0; i<24; i++) 
-      matching_plifs[i] = match_plifs[i];
+void Alignment::setMatchPlifs(struct penalty_struct match_plifs, int idx) {
+   struct penalty_struct *currentPlif = new struct penalty_struct();
+   copy_penalty_struct(&match_plifs,currentPlif);
+   ::qualityScores[idx] = (*currentPlif);
 }
 
 void Alignment::getDNAEST(){}

diff --git a/QPalmaDP/qpalma_dp.h b/QPalmaDP/qpalma_dp.h
index e38469f..7f9817e 100644 (file)
--- a/QPalmaDP/qpalma_dp.h
+++ b/QPalmaDP/qpalma_dp.h
@@ -29,7 +29,6 @@ typedef struct Align_pair { //24B
   char dna_char;
 };
 
-
 void fill_matrix (int nr_paths, Pre_score* matrices[], int len_est, int len_dna, char* est, char* dna, penalty_struct* functions, double* matchmatrix, double* donor, double* acceptor, bool remove_duplicate_scores, int nr_donor_sites, int* donor_sites, int* max_score_positions);
 
 int check_char(char base);
@@ -107,7 +106,7 @@ class Alignment {
       double* donor, int d_len, double* acceptor, int a_len,
       bool remove_duplicate_scores, bool print_matrix);
 
-      void setMatchPlifs(struct penalty_struct* match_plifs);
+      void setMatchPlifs(struct penalty_struct match_plifs,int idx);
 
       void penSetLength(int l) { len = l; }
       void penSetLimits(REAL* lts) { 

diff --git a/python/qpalma.py b/python/qpalma.py
index 43d4e24..57eccef 100644 (file)
--- a/python/qpalma.py
+++ b/python/qpalma.py
@@ -34,9 +34,12 @@ import Configuration
 
 
 
-def initializeQualityScoringFunctions(numPlifs,numSuppPoints):
+def initializeQualityScoringFunctions(alignment,numPlifs,numSuppPoints):
 
-   allPlifs = [0]*numPlifs
+   min_intron_len=20
+   max_intron_len=1000
+   min_svm_score=-5
+   max_svm_score=5
 
    for idx in range(numPlifs):
 
@@ -51,9 +54,8 @@ def initializeQualityScoringFunctions(numPlifs,numSuppPoints):
       curPlif.use_svm   = 0 
       curPlif.next_id   = 0 
 
-      allPlifs[idx] = d
-
-   return allPlifs
+      curPlif = curPlif.convert2SWIG()
+      alignment.setMatchPlifs(curPlif,idx)
 
 
 class QPalma:
@@ -172,8 +174,6 @@ class QPalma:
 
       numPlifs = 24
       numSuppPoints = 30
-      allMatchingPlifs = initializeQualityScoringFunctions(numPlifs,numSuppPoints)
-      # self.numQualSuppPoints
 
       #############################################################################################
       # Training
@@ -252,7 +252,7 @@ class QPalma:
             qualityMat = QPalmaDP.createDoubleArrayFromList(qualityMatrix)
             currentAlignment.setQualityMatrix(qualityMat,self.numQualSuppPoints)
 
-            currentAlignment.setMatchPlifs()
+            initializeQualityScoringFunctions(currentAlignment,numPlifs,numSuppPoints)
 
             #print 'PYTHON: Calling myalign...'
             # calculates SpliceAlign, EstAlign, weightMatch, Gesamtscores, dnaest