+ minor changes
authorfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Thu, 17 Jan 2008 14:17:01 +0000 (14:17 +0000)
committerfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Thu, 17 Jan 2008 14:17:01 +0000 (14:17 +0000)
git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@7496 e1793c9e-67f9-0310-80fc-b846ff1f7b36

QPalmaDP/fill_matrix.cpp
python/qpalma.py
python/set_param_palma.py

index ca6c851..e226266 100644 (file)
@@ -72,7 +72,7 @@ int check_char(char base)
  * quality score.
  */
 
-double getScore(struct penalty_struct* qualityScores, double* matchmatrix, int mlen, int dnaChar, int estChar, double baseScore) {
+double getScore(struct penalty_struct* qualityScores, int mlen, int dnaChar, int estChar, double baseScore) {
    double score = 0.0;
    FA(estChar > 0);
    FA(dnaChar > -1);
@@ -83,7 +83,23 @@ double getScore(struct penalty_struct* qualityScores, double* matchmatrix, int m
    int currentPos = (estChar-1)*6+dnaChar;
    struct penalty_struct currentPen = qualityScores[currentPos];
 
-   score = lookup_penalty(&currentPen,baseScore,NULL,false);
+   //currentPen.transform = T_LINEAR;
+   //score = lookup_penalty(&currentPen,(int)baseScore,NULL,false);
+
+       INT idx = 0 ;
+       for (INT i=0; i<currentPen.len; i++)
+               if (currentPen.limits[i]<=baseScore)
+                       idx++ ;
+       
+       if (idx==0)
+               score=currentPen.penalties[0] ;
+       else if (idx==currentPen.len)
+               score=currentPen.penalties[currentPen.len-1] ;
+       else
+       {
+               score = (currentPen.penalties[idx]*(baseScore-currentPen.limits[idx-1]) + currentPen.penalties[idx-1]*
+                          (currentPen.limits[idx]-baseScore)) / (currentPen.limits[idx]-currentPen.limits[idx-1]) ;  
+   }
 
    //if(baseScore <= currentPen.limits[0]) {
    //   score = currentPen.penalties[0];
@@ -295,10 +311,10 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
        if (isfinite(prevValue))
        {
 
-      if( currentMode == NORMAL ) {
-             tempValue = prevValue +(matchmatrix[mlen* dnaChar]);   /* score(gap,DNA) */
-      } else {
+      if( currentMode == USE_QUALITY_SCORES ) {
          tempValue = prevValue + matchmatrix[dnaChar];
+      } else {
+             tempValue = prevValue + matchmatrix[mlen*dnaChar];   /* score(gap,DNA) */
       }
 
 
@@ -321,10 +337,10 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
        prevValue = ((Pre_score*)actMatrix +(i-1)*dna_len + j-1)->value ; /* (j-1,i-1) -> (j,i) */
        if (isfinite(prevValue))
        {
-      if( currentMode == NORMAL ) {
-        tempValue = prevValue +(matchmatrix[mlen* dnaChar +estChar]);
+      if( currentMode == USE_QUALITY_SCORES ) {
+        tempValue = prevValue + getScore(qualityScores,mlen,dnaChar,estChar,baseScore);
       } else {
-         tempValue = prevValue + getScore(qualityScores,matchmatrix,mlen,dnaChar,estChar,baseScore);
+        tempValue = prevValue +(matchmatrix[mlen* dnaChar +estChar]);
      }
 
          if (isfinite(tempValue))
@@ -347,10 +363,10 @@ void fill_matrix(int nr_paths, Pre_score* matrices[], int est_len /*counter var:
        if (isfinite(prevValue))
        {
 
-      if( currentMode == NORMAL ) {
-             tempValue = prevValue +(matchmatrix[estChar]); /* score(EST,gap) */
+      if( currentMode == USE_QUALITY_SCORES ) {
+        tempValue = prevValue + getScore(qualityScores,mlen,0,estChar,baseScore);
       } else {
-         tempValue = prevValue + getScore(qualityScores,matchmatrix,mlen,0,estChar,baseScore);
+            tempValue = prevValue +(matchmatrix[estChar]); /* score(EST,gap) */
      }
 
          if (isfinite(tempValue))
index 8dcc48f..1be8ea1 100644 (file)
@@ -74,7 +74,6 @@ class QPalma:
          #Qualities = []
          #for i in range(len(Ests)):
          #   Qualities.append([40]*len(Ests[i]))
-
          use_quality_scores = True
       else:
          assert(False)
@@ -205,8 +204,8 @@ class QPalma:
             dna_len = len(dna)
             est_len = len(est)
 
-            #dna ='n'*dna_len
-            #est ='g'*est_len
+            dna ='t'*dna_len
+            est ='a'*est_len
             #quality = [40]*len(quality)
 
             #mmatrix = zeros((6,1))
@@ -305,6 +304,7 @@ class QPalma:
                for qidx in range(Configuration.totalQualSuppPoints):
                   decodedQualityFeatures[qidx] = newQualityPlifsFeatures[(pathNr*Configuration.totalQualSuppPoints)+qidx]
 
+               path_loss[pathNr] = 0
                # sum up positionwise loss between alignments
                for alignPosIdx in range(len(newSpliceAlign[pathNr,:])):
                   if newSpliceAlign[pathNr,alignPosIdx] != trueSpliceAlign[alignPosIdx]:
index c6f20ba..f760031 100644 (file)
@@ -124,7 +124,6 @@ def set_param_palma(param, train_with_intronlengthinformation,\
       currentPlif.limits    = linspace(Configuration.min_qual,Configuration.max_qual,Configuration.numQualSuppPoints) 
       begin                 = lengthSP+donSP+accSP+mmatrixSP+(idx*Configuration.numQualSuppPoints)
       end                   = lengthSP+donSP+accSP+mmatrixSP+((idx+1)*Configuration.numQualSuppPoints)
-      #print begin,end
       currentPlif.penalties = param[begin:end].flatten().tolist()[0]
       currentPlif.transform = '' 
       currentPlif.name      = 'q'