git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@8608 e1793c9e...
authorraetsch <raetsch@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Thu, 17 Apr 2008 20:03:11 +0000 (20:03 +0000)
committerraetsch <raetsch@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Thu, 17 Apr 2008 20:03:11 +0000 (20:03 +0000)
qpalma/computeSpliceAlignWithQuality.py
scripts/PipelineHeuristic.py

index 4734ac4..3bbfe35 100644 (file)
@@ -175,7 +175,7 @@ def  computeSpliceAlignWithQuality(dna, exons, est, original_est, quality, quali
    return SpliceAlign, trueWeightMatch, trueWeightQuality, dna_calc
 
 
-def  computeSpliceAlignScoreWithQuality(dna, exons, est, original_est, quality, qualityPlifs, run, weightvector):
+def  computeSpliceAlignScoreWithQuality(original_est, quality, qualityPlifs, run, weightvector):
    """
    """
 
@@ -192,19 +192,9 @@ def  computeSpliceAlignScoreWithQuality(dna, exons, est, original_est, quality,
    numqP = weightvector[accSP+donSP+lengthSP+mmatrixSP:accSP+donSP+lengthSP+mmatrixSP+numq]
    numQualSuppPoints=run['numQualSuppPoints']
    
-   map = {'-':0, 'a':1, 'c':2, 'g':3, 't':4, 'n':5, '[': 10, ']':11}
+   map = {'-':0, 'a':1, 'c':2, 'g':3, 't':4, 'n':5, '[': 10, ']':11 }
    original_est_mapped = [ map[e] for e in original_est ]
 
-   if exons.shape == (2,2):
-      numberOfExons = 2
-      exonSizes = [-1]*numberOfExons
-
-      exonSizes[0] = exons[0,1] - exons[0,0]
-      exonSizes[1] = exons[1,1] - exons[1,0]
-
-      score=0.0 
-
-   elif exons.shape == (2,1):
        score=0.0
        est_ptr=0
        dna_ptr=0
@@ -219,9 +209,9 @@ def  computeSpliceAlignScoreWithQuality(dna, exons, est, original_est, quality,
                dnaletter=original_est_mapped[ptr]
                estletter=dnaletter
                ptr+=1
+
            if estletter==0:# '-':
                score+= mmatrixP[dna_letter]
-               dna_ptr+=1 
            else:
                value = quality[est_ptr]
                cur_plf = qualityPlifs[(estletter-1)*6+dnaletter]
@@ -257,8 +247,5 @@ def  computeSpliceAlignScoreWithQuality(dna, exons, est, original_est, quality,
                
            assert(dna_ptr<=len(dna))
            assert(est_ptr<=len(est))
-           
-   else:
-       assert(False)
 
    return score
index d982f79..438e525 100644 (file)
@@ -370,9 +370,6 @@ class PipelineHeuristic:
          computeSpliceAlignWithQuality(dna, exons, est, original_est,\
          quality, qualityPlifs, run)
 
-         score = computeSpliceAlignScoreWithQuality(dna, exons, est, original_est,\
-         quality, qualityPlifs, run, self.currentPhi)
-
          stop = cpu()
          self.computeSpliceAlignWithQualityTime += stop-start
          start = cpu()
@@ -391,6 +388,43 @@ class PipelineHeuristic:
          # Calculate w'phi(x,y) the total score of the alignment
          alternativeScores.append((trueWeight.T * self.currentPhi)[0,0])
 
+         original_est_cut=''
+
+         est_ptr=0
+         dna_ptr=0
+         ptr=0 
+         while ptr<len(original_est):
+             
+            if original_est[ptr]=='[':
+                dnaletter=original_est[ptr+1]
+                estletter=original_est[ptr+2]
+                if est_ptr<=exons[0,1]:
+                    original_est_cut.join(original_est[ptr:ptr+4]) 
+                else:
+                    original_est_cut.join(estletter) # EST letter
+                ptr+=4 
+            else:
+                dnaletter=original_est[ptr]
+                estletter=dnaletter
+                
+                original_est_cut.join(estletter) # EST letter
+                ptr+=1
+
+            if estletter=='-':
+                dna_ptr+=1 
+            elif dnaletter=='-':
+                est_ptr+=1
+            else:
+                dna_ptr+=1 
+                est_ptr+=1
+                         
+         assert(dna_ptr<=len(dna))
+         assert(est_ptr<=len(est))
+
+         
+         # new score
+         score = computeSpliceAlignScoreWithQuality(original_est_cut, quality, qualityPlifs, run, self.currentPhi)
+
          if score!=0.0:
              print 'diff %f' % ((trueWeight.T * self.currentPhi)[0,0] - score)
 
@@ -416,8 +450,7 @@ class PipelineHeuristic:
       # Lets start calculation
       dna, exons, est, original_est, quality, acc_supp, don_supp = alignment
 
-      score = computeSpliceAlignScoreWithQuality(dna, exons, est, original_est,\
-         quality, self.qualityPlifs, run, self.currentPhi)
+      score = computeSpliceAlignScoreWithQuality(original_est, quality, self.qualityPlifs, run, self.currentPhi)
 
       stop = cpu()
       self.calcAlignmentScoreTime += stop-start