+ removed duplicated code
authorfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Mon, 17 Mar 2008 13:08:41 +0000 (13:08 +0000)
committerfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Mon, 17 Mar 2008 13:08:41 +0000 (13:08 +0000)
+ found bug in the evaluation -> wrong position calculation for chr 6 and 7

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

scripts/qpalma_main.py

index 4180006..a7a2e1f 100644 (file)
@@ -29,8 +29,8 @@ from numpy.linalg import norm
 import QPalmaDP
 import qpalma
 
-#from qpalma.SIQP_CPX import SIQPSolver
-from qpalma.SIQP_CVXOPT import SIQPSolver
+from qpalma.SIQP_CPX import SIQPSolver
+#from qpalma.SIQP_CVXOPT import SIQPSolver
 
 from qpalma.DataProc import *
 from qpalma.computeSpliceWeights import *
@@ -69,6 +69,69 @@ def unbracket_est(est):
    return "".join(new_est).lower()
 
 
+def getData(SeqInfo,OriginalEsts,Exons,exampleIdx,run):
+   currentSeqInfo = SeqInfo[exampleIdx]
+   chr,strand,up_cut,down_cut = currentSeqInfo 
+
+   est = OriginalEsts[exampleIdx] 
+   est = "".join(est)
+   est = est.lower()
+   est = unbracket_est(est)
+   est = est.replace('-','')
+
+   assert len(est) == run['read_size'], pdb.set_trace()
+   est_len = len(est)
+
+   original_est = OriginalEsts[exampleIdx] 
+   original_est = "".join(original_est)
+   original_est = original_est.lower()
+
+   dna_flat_files =  '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
+   dna, acc_supp, don_supp = get_seq_and_scores(chr,strand,up_cut,down_cut,dna_flat_files)
+
+   # splice score is located at g of ag
+   ag_tuple_pos = [p for p,e in enumerate(dna) if p>1 and dna[p-1]=='a' and dna[p]=='g' ]
+   assert ag_tuple_pos == [p for p,e in enumerate(acc_supp) if e != -inf and p > 1], pdb.set_trace()
+
+   gt_tuple_pos = [p for p,e in enumerate(dna) if p>0 and p<len(dna)-1 and e=='g' and (dna[p+1]=='t' or dna[p+1]=='c')]
+   assert gt_tuple_pos == [p for p,e in enumerate(don_supp) if e != -inf and p > 0], pdb.set_trace()
+
+   # check whether the whole business is really related to invalid
+   # splice scores
+   #acc_supp = [0.0]*len(acc_supp)
+   #don_supp = [0.0]*len(don_supp)
+
+   original_exons = Exons[exampleIdx]
+
+   exons = original_exons - (up_cut-1)
+   exons[0,0] -= 1
+   exons[1,0] -= 1
+
+   fetched_dna_subseq = dna[exons[0,0]:exons[0,1]] + dna[exons[1,0]:exons[1,1]]
+   
+   donor_elem = dna[exons[0,1]:exons[0,1]+2]
+   acceptor_elem = dna[exons[1,0]-2:exons[1,0]]
+
+   if not ( donor_elem == 'gt' or donor_elem == 'gc' ):
+      print 'invalid donor in example %d'% exampleIdx
+
+   if not ( acceptor_elem == 'ag' ):
+      print 'invalid acceptor in example %d'% exampleIdx
+
+   assert len(fetched_dna_subseq) == len(est), pdb.set_trace()
+
+   #new_string = ''
+   #for idx in range(len(est)):
+   #   dna_char = fetched_dna_subseq[idx] 
+   #   est_char = est[idx]
+   #   if dna_char == est_char:
+   #      new_string += est_char
+   #   else:
+   #      new_string += '[%s%s]'%(dna_char,est_char)
+
+   return dna,est,acc_supp,don_supp,exons,original_est
+
+
 class QPalma:
    """
    A training method for the QPalma project
@@ -309,32 +372,14 @@ class QPalma:
             if (exampleIdx%100) == 0:
                print 'Current example nr %d' % exampleIdx
 
-            currentSeqInfo = SeqInfo[exampleIdx]
-            chr,strand,up_cut,down_cut = currentSeqInfo 
-
-            est = OriginalEsts[exampleIdx] 
-            est = "".join(est)
-            est = est.lower()
-            est = unbracket_est(est)
-            est = est.replace('-','')
-
-            assert len(est) == run['read_size'], pdb.set_trace()
-            est_len = len(est)
+            dna,est,acc_supp,don_supp,exons,original_est =\
+            getData(SeqInfo,OriginalEsts,Exons,exampleIdx,run)
 
-            original_est = OriginalEsts[exampleIdx] 
-            original_est = "".join(original_est)
-            original_est = original_est.lower()
-
-            dna_flat_files =  '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
-            dna, acc_supp, don_supp = get_seq_and_scores(chr,strand,up_cut,down_cut,dna_flat_files)
             dna_len = len(dna)
 
-            # splice score is located at g of ag
-            ag_tuple_pos = [p for p,e in enumerate(dna) if p>1 and dna[p-1]=='a' and dna[p]=='g' ]
-            assert ag_tuple_pos == [p for p,e in enumerate(acc_supp) if e != -inf and p > 1], pdb.set_trace()
-
-            gt_tuple_pos = [p for p,e in enumerate(dna) if p>0 and p<len(dna)-1 and e=='g' and (dna[p+1]=='t' or dna[p+1]=='c')]
-            assert gt_tuple_pos == [p for p,e in enumerate(don_supp) if e != -inf and p > 0], pdb.set_trace()
+            #if new_string != original_est:
+            #   print new_string,original_est
+            #   continue
 
             if run['mode'] == 'normal':
                quality = [40]*len(est)
@@ -345,42 +390,6 @@ class QPalma:
             if not run['enable_quality_scores']:
                quality = [40]*len(est)
 
-            # check whether the whole business is really related to invalid
-            # splice scores
-            #acc_supp = [0.0]*len(acc_supp)
-            #don_supp = [0.0]*len(don_supp)
-
-            original_exons = Exons[exampleIdx]
-            exons = original_exons - (up_cut-1)
-            exons[0,0] -= 1
-            exons[1,0] -= 1
-
-            fetched_dna_subseq = dna[exons[0,0]:exons[0,1]] + dna[exons[1,0]:exons[1,1]]
-            
-            donor_elem = dna[exons[0,1]:exons[0,1]+2]
-            acceptor_elem = dna[exons[1,0]-2:exons[1,0]]
-
-            if not ( donor_elem == 'gt' or donor_elem == 'gc' ):
-               print 'invalid donor in example %d'% exampleIdx
-
-            if not ( acceptor_elem == 'ag' ):
-               print 'invalid acceptor in example %d'% exampleIdx
-   
-            assert len(fetched_dna_subseq) == len(est), pdb.set_trace()
-
-            new_string = ''
-            for idx in range(len(est)):
-               dna_char = fetched_dna_subseq[idx] 
-               est_char = est[idx]
-               if dna_char == est_char:
-                  new_string += est_char
-               else:
-                  new_string += '[%s%s]'%(dna_char,est_char)
-
-            if new_string != original_est:
-               print new_string,original_est
-               continue
-
             if not run['enable_splice_signals']:
                for idx,elem in enumerate(don_supp):
                   if elem != -inf:
@@ -732,29 +741,16 @@ class QPalma:
          currentSeqInfo = SeqInfo[exampleIdx]
          chr,strand,up_cut,down_cut = currentSeqInfo 
 
-         est = OriginalEsts[exampleIdx] 
-         est = "".join(est)
-         est = est.lower()
-         est = unbracket_est(est)
-         est = est.replace('-','')
-
-         assert len(est) == run['read_size'], pdb.set_trace()
-         est_len = len(est)
-
-         original_est = OriginalEsts[exampleIdx] 
-         original_est = "".join(original_est)
-         original_est = original_est.lower()
+         dna,est,acc_supp,don_supp,exons,original_est =\
+         getData(SeqInfo,OriginalEsts,Exons,exampleIdx,run)
 
-         dna_flat_files =  '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
-         dna, acc_supp, don_supp = get_seq_and_scores(chr,strand,up_cut,down_cut,dna_flat_files)
          dna_len = len(dna)
 
-         # splice score is located at g of ag
-         ag_tuple_pos = [p for p,e in enumerate(dna) if p>1 and dna[p-1]=='a' and dna[p]=='g' ]
-         assert ag_tuple_pos == [p for p,e in enumerate(acc_supp) if e != -inf and p > 1], pdb.set_trace()
-
-         gt_tuple_pos = [p for p,e in enumerate(dna) if p>0 and p<len(dna)-1 and e=='g' and (dna[p+1]=='t' or dna[p+1]=='c')]
-         assert gt_tuple_pos == [p for p,e in enumerate(don_supp) if e != -inf and p > 0], pdb.set_trace()
+         #if new_string != original_est:
+         #   print 'seq. inconsistency'
+         #   print new_string,original_est
+         #   print exampleIdx
+         #   continue
 
          if run['mode'] == 'normal':
             quality = [40]*len(est)
@@ -765,43 +761,6 @@ class QPalma:
          if not run['enable_quality_scores']:
             quality = [40]*len(est)
 
-         # check whether the whole business is really related to invalid
-         # splice scores
-         #acc_supp = [0.0]*len(acc_supp)
-         #don_supp = [0.0]*len(don_supp)
-
-         original_exons = Exons[exampleIdx]
-         exons = original_exons - (up_cut-1)
-         exons[0,0] -= 1
-         exons[1,0] -= 1
-
-         fetched_dna_subseq = dna[exons[0,0]:exons[0,1]] + dna[exons[1,0]:exons[1,1]]
-         
-         donor_elem = dna[exons[0,1]:exons[0,1]+2]
-         acceptor_elem = dna[exons[1,0]-2:exons[1,0]]
-
-         if not ( donor_elem == 'gt' or donor_elem == 'gc' ):
-            print 'invalid donor in example %d'% exampleIdx
-
-         if not ( acceptor_elem == 'ag' ):
-            print 'invalid acceptor in example %d'% exampleIdx
-
-         assert len(fetched_dna_subseq) == len(est), pdb.set_trace()
-
-         new_string = ''
-         for idx in range(len(est)):
-            dna_char = fetched_dna_subseq[idx] 
-            est_char = est[idx]
-            if dna_char == est_char:
-               new_string += est_char
-            else:
-               new_string += '[%s%s]'%(dna_char,est_char)
-
-         if new_string != original_est:
-            print new_string,original_est
-            print exampleIdx
-            continue
-
          if not run['enable_splice_signals']:
             for idx,elem in enumerate(don_supp):
                if elem != -inf:
@@ -826,6 +785,10 @@ class QPalma:
          # the vmatch results
          for alternative_alignment in currentAlternatives:
             chr, strand, genomicSeq_start, genomicSeq_stop, currentLabel = alternative_alignment
+
+            if not chr in range(1,6):
+               continue
+
             currentDNASeq, currentAcc, currentDon = get_seq_and_scores(chr,strand,genomicSeq_start,genomicSeq_stop,run['dna_flat_files'])
 
             if not run['enable_splice_signals']: