+ changed filter behavior (tendency towards smaller read fragments)
authorfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Sun, 10 Feb 2008 22:38:55 +0000 (22:38 +0000)
committerfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Sun, 10 Feb 2008 22:38:55 +0000 (22:38 +0000)
+ removed monotonicity assertion for objective

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

qpalma/Configuration.py
qpalma/SIQP_CPX.py
scripts/compile_dataset.py
scripts/qpalma_main.py
tools/data_tools/filterReads.c

index 2e06711..0cd82a3 100644 (file)
@@ -13,7 +13,7 @@ import cPickle
 
 fixedParamQ = cPickle.load(open('/fml/ag-raetsch/home/fabio/svn/projects/QPalma/scripts/randInitParam.pickle'))
 
-###########################################################
+###############################################################################
 #
 # The parameters for the QPalma algorithm
 #
@@ -76,7 +76,7 @@ mode = 'using_quality_scores'
 
 read_size = 36
 
-extension = (1,20)
+extension = (250,500)
 
 #numLengthSuppPoints  = 30 
 #numDonSuppPoints     = 30
@@ -89,8 +89,8 @@ extension = (1,20)
 #numQualSuppPoints    =  8
 
 numLengthSuppPoints  = 10
-numDonSuppPoints     = 10
-numAccSuppPoints     = 10
+numDonSuppPoints     = 2
+numAccSuppPoints     = 2
 numQualSuppPoints    = 10
 
 min_qual = -5
@@ -116,7 +116,6 @@ totalQualSuppPoints = numQualPlifs*numQualSuppPoints
 numFeatures = numDonSuppPoints + numAccSuppPoints\
 + numLengthSuppPoints + sizeMatchmatrix[0]*sizeMatchmatrix[1] + totalQualSuppPoints 
 
-
 ###############################################################################
 #
 # GENERAL SETTINGS CONCERNING THE SOLVER
@@ -144,12 +143,11 @@ else:
 #
 #
 ###############################################################################
-
 training_begin    =  0
 training_end      =  500
 
 prediction_begin  =  500
-prediction_end    =  2000
+prediction_end    =  5000
 
 joinp = os.path.join
 
@@ -185,3 +183,4 @@ assert os.path.exists(remapped_fn), 'EST/Reads data does not exist!'
 
 extended_alphabet = ['-','a','c','g','t','n','[',']']
 alphabet          = ['-','a','c','g','t','n']
+
index e77b92e..d6f859d 100644 (file)
@@ -233,8 +233,8 @@ class SIQPSolver(SIQP):
          print >> self.protocol, 'SIQP_CPX: Objective is: %f'%objval
 
          # Check that objective value is monotonically increasing
-         if math.fabs(objval - self.old_objective_value) >= 10e-5:
-            assert objval >= self.old_objective_value
+         #if math.fabs(objval - self.old_objective_value) >= 1e-4:
+         #   assert objval >= self.old_objective_value
 
          self.old_objective_value = objval
 
index 0b6b0bf..9464c3b 100644 (file)
@@ -21,6 +21,14 @@ from genome_utils import load_genomic
 
 import qpalma.Configuration as Conf
 
+warning = """
+   #########################################################
+   WARNING ! Disabled the use of non-bracket reads WARNING !
+   #########################################################
+   
+   """
+
+
 help = """
 
    Usage of this script is:
@@ -141,7 +149,7 @@ def compile_d(gff_file,dna_flat_files,filtered_reads,remapped_reads,tmp_dir,data
       if instance_counter % 100 == 0:
          print 'processed %d examples' % instance_counter
 
-      if instance_counter == 2000:
+      if instance_counter == 5000:
          break
 
    print 'Full dataset has size %d' % len(Sequences)
@@ -166,13 +174,16 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
 
    chr            = fRead['chr']
    strand         = fRead['strand']
-   quality        = fRead['prb']
+   #quality        = fRead['prb']
+   quality        = fRead['cal_prb']
+   #quality        = fRead['chastity']
    spos           = fRead['splitpos']
    currentReadSeq = fRead['seq']
    currentReadSeq = currentReadSeq.lower()
 
-   if currentReadSeq.find('[') == -1:
-      return nil
+   #if currentReadSeq.find('[') == -1:
+   #   print warning
+   #   return nil
 
    originalReadSeq = fRead['seq']
 
@@ -187,12 +198,17 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
 
    seq_has_indels = False
    for elem in currentReadSeq:
-      assert elem in extended_alphabet, 'Invalid char (extended_alphabet) in read: \"%s\"' % elem
+      #assert elem in extended_alphabet, 'Invalid char (extended_alphabet) in read: \"%s\"' % elem
+      if not elem in extended_alphabet:
+         return nil
+
       if elem ==  ']':
          seq_has_indels = True
 
    for elem in currentDNASeq:
-      assert elem in alphabet, 'Invalid char in dna: \"%s\"' % elem
+      #assert elem in alphabet, 'Invalid char in dna: \"%s\"' % elem
+      if not elem in alphabet:
+         return nil
 
    if strand == '-':
       try:
@@ -305,14 +321,18 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
 
       if not dna_fragment_1.replace('-','') == dna_annot_1:
          print 'genomic seq mismatch'
+         print 'orignal read:\t %s ' % originalReadSeq
+         print '1st/2nd read:\t %s %s' % (firstReadSeq,secondReadSeq)
+         print 'dna_fragment_1:\t %s' % dna_fragment_1
+         print 'dna_annot 1+2:\t %s %s' % (dna_annot_1,dna_annot_2)
          return nil
 
       if not dna_fragment_2.replace('-','') == dna_annot_2:
          print 'genomic seq mismatch'
-         print originalReadSeq
-         print firstReadSeq
-         print secondReadSeq
-         print dna_annot_1 + dna_annot_2
+         print 'orignal read:\t %s ' % originalReadSeq
+         print '1st/2nd read:\t %s %s' % (firstReadSeq,secondReadSeq)
+         print 'dna_fragment_2:\t %s' % dna_fragment_2
+         print 'dna_annot 1+2:\t %s %s' % (dna_annot_1,dna_annot_2)
          return nil
 
       #print 'successful'
@@ -320,7 +340,7 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
       for elem in currentReadSeq:
          assert elem in alphabet, 'Invalid char (alphabet) in read: \"%s\"' % elem
 
-   assert p_start < exon_stop < exon_start < p_stop, 'Invalid Indices'
+   #assert p_start < exon_stop < exon_start < p_stop, pdb.set_trace()
    assert p_stop - p_start >= Conf.read_size
 
    currentExons = zeros((2,2),dtype=numpy.int)
@@ -373,7 +393,8 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
          return nil
 
    except:
-      pdb.set_trace()
+      return nil
+      #pdb.set_trace()
 
    # now we want to use interval_query to get the predicted splice scores
    # trained on the TAIR sequence and annotation
@@ -436,19 +457,30 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
             continue
          don[position-1] = scores[i]
    except:
-      pdb.set_trace()
+      return nil
+      #pdb.set_trace()
 
    don = [-inf] + don[:-1]
 
-   assert ag_tuple_pos == [p for p,e in enumerate(acc) if e != -inf], pdb.set_trace()
-   assert gt_tuple_pos == [p for p,e in enumerate(don) if e != -inf], pdb.set_trace()
+   #assert ag_tuple_pos == [p for p,e in enumerate(acc) if e != -inf], pdb.set_trace()
+   #assert gt_tuple_pos == [p for p,e in enumerate(don) if e != -inf], pdb.set_trace()
+   #assert (len(currentDNASeq)) == len(acc), pdb.set_trace()
+   #assert (len(currentDNASeq)) == len(don), pdb.set_trace()
 
-   assert (len(currentDNASeq)) == len(acc), pdb.set_trace()
-   assert (len(currentDNASeq)) == len(don), pdb.set_trace()
 
-   acc = [-inf] + acc[:-1]
+   if not ag_tuple_pos == [p for p,e in enumerate(acc) if e != -inf]:
+      return nil
 
+   if not gt_tuple_pos == [p for p,e in enumerate(don) if e != -inf]:
+      return nil
 
+   if not (len(currentDNASeq)) == len(acc):
+      return nil
+
+   if not (len(currentDNASeq)) == len(don):
+      return nil
+
+   acc = [-inf] + acc[:-1]
 
    return currentDNASeq, acc, don, currentExons, currentReadSeq, originalReadSeq, quality, spos
 
index f0b7893..13b69fd 100644 (file)
@@ -59,8 +59,6 @@ class QPalma:
       #compile_d(Conf.gff_fn,Conf.dna_flat_fn,Conf.filtered_fn,Conf.remapped_fn,Conf.tmp_dir,Conf.dataset_fn,False)
       #compile_d(Conf.gff_fn,Conf.dna_flat_fn,Conf.filtered_fn,Conf.remapped_fn,Conf.tmp_dir,Conf.dataset_fn,True)
 
-      pdb.set_trace()
-
       data_filename = Conf.dataset_fn
       Sequences, Acceptors, Donors, Exons, Ests, OriginalEsts, Qualities, SplitPos = paths_load_data(data_filename,'training',None,self.ARGS)
 
@@ -82,14 +80,15 @@ class QPalma:
       self.Donors      = Donors
       self.Acceptors   = Acceptors
 
-      calc_info(self.Acceptors,self.Donors,self.Exons)
+      calc_info(self.Acceptors,self.Donors,self.Exons,self.Qualities)
+      pdb.set_trace()
       print 'leaving constructor...'
 
    def plog(self,string):
       self.logfh.write(string)
       self.logfh.flush()
 
-   def train(self,value):
+   def train(self):
       self.logfh = open('_qpalma_train.log','w+')
 
       beg = Conf.training_begin
@@ -129,13 +128,13 @@ class QPalma:
       numq        = Conf.numQualSuppPoints
       totalQualSP = Conf.totalQualSuppPoints
 
-      i1 = range(0,lengthSP)
-      i2 = range(lengthSP,lengthSP+donSP)
-      i3 = range(lengthSP+donSP,lengthSP+donSP+accSP)
-      i4 = range(lengthSP+donSP+accSP,lengthSP+donSP+accSP+mmatrixSP)
-      i5 = range(lengthSP+donSP+accSP+mmatrixSP,lengthSP+donSP+accSP+mmatrixSP+totalQualSP)
+      #i1 = range(0,lengthSP)
+      #i2 = range(lengthSP,lengthSP+donSP)
+      #i3 = range(lengthSP+donSP,lengthSP+donSP+accSP)
+      #i4 = range(lengthSP+donSP+accSP,lengthSP+donSP+accSP+mmatrixSP)
+      #i5 = range(lengthSP+donSP+accSP+mmatrixSP,lengthSP+donSP+accSP+mmatrixSP+totalQualSP)
       
-      offset =lengthSP+donSP+accSP+mmatrixSP
+      #offset =lengthSP+donSP+accSP+mmatrixSP
       #intervals = [[offset+2,offset+3]] #i5,\ #i2,#i3,#i4,#i5]
       #intervals = []
 
@@ -202,20 +201,24 @@ class QPalma:
             if Conf.mode == 'using_quality_scores':
                quality = Qualities[exampleIdx]
 
-            #quality = [(int(math.fabs(e))) for e in quality]
+            # if Conf.disable_quality_scores
             #quality = [40]*len(est)
 
             exons = Exons[exampleIdx] 
             don_supp = Donors[exampleIdx] 
             acc_supp = Acceptors[exampleIdx] 
 
-            #for idx,elem in enumerate(don_supp):
-            #   if elem != -inf:
-            #      don_supp[idx] = 0.0
+            # if Conf.disable_splice_signals
+            #don_supp = [0.0]*len(don_supp)
+            #acc_supp = [0.0]*len(acc_supp)
 
-            #for idx,elem in enumerate(acc_supp):
-            #   if elem != -inf:
-            #      acc_supp[idx] = 0.0
+            for idx,elem in enumerate(don_supp):
+               if elem != -inf:
+                  don_supp[idx] = 0.0
+
+            for idx,elem in enumerate(acc_supp):
+               if elem != -inf:
+                  acc_supp[idx] = 0.0
 
             # Berechne die Parameter des wirklichen Alignments (but with untrained d,a,h ...)    
             if Conf.mode == 'using_quality_scores':
@@ -262,12 +265,12 @@ class QPalma:
             # check that splice site scores are at dna positions as expected by
             # the dynamic programming component
 
-            for d_pos in [pos for pos,elem in enumerate(donor) if elem != -inf]:
-               assert dna[d_pos] == 'g' and (dna[d_pos+1] == 'c'\
-               or dna[d_pos+1] == 't'), pdb.set_trace()
-                
-            for a_pos in [pos for pos,elem in enumerate(acceptor) if elem != -inf]:
-               assert dna[a_pos-1] == 'a' and dna[a_pos] == 'g', pdb.set_trace()
+            #for d_pos in [pos for pos,elem in enumerate(donor) if elem != -inf]:
+            #   assert dna[d_pos] == 'g' and (dna[d_pos+1] == 'c'\
+            #   or dna[d_pos+1] == 't'), pdb.set_trace()
+            #    
+            #for a_pos in [pos for pos,elem in enumerate(acceptor) if elem != -inf]:
+            #   assert dna[a_pos-1] == 'a' and dna[a_pos] == 'g', pdb.set_trace()
 
             ps = h.convert2SWIG()
 
@@ -515,6 +518,7 @@ class QPalma:
    def evaluate(self,param_filename):
       beg = Conf.prediction_begin
       end = Conf.prediction_end
+      self.logfh = open('_qpalma_predict.log','w+')
 
       # predict on training set
       print '##### Prediction on the training set #####'
@@ -524,10 +528,10 @@ class QPalma:
       print '##### Prediction on the test set #####'
       self.predict(param_filename,beg,end)
    
+      self.logfh.close()
       #pdb.set_trace()
 
    def predict(self,param_filename,beg,end):
-      self.logfh = open('_qpalma_predict.log','w+')
 
       Sequences   = self.Sequences[beg:end]
       Exons       = self.Exons[beg:end]
@@ -603,6 +607,9 @@ class QPalma:
          don_supp = Donors[exampleIdx] 
          acc_supp = Acceptors[exampleIdx] 
 
+         #don_supp = [0.0]*len(don_supp)
+         #acc_supp = [0.0]*len(acc_supp)
+
          # for idx,elem in enumerate(don_supp):
          #    if elem != -inf:
          #       don_supp[idx] = 0.0
@@ -611,6 +618,14 @@ class QPalma:
          #    if elem != -inf:
          #       acc_supp[idx] = 0.0
 
+         for idx,elem in enumerate(don_supp):
+            if elem != -inf:
+               don_supp[idx] = 0.0
+
+         for idx,elem in enumerate(acc_supp):
+            if elem != -inf:
+               acc_supp[idx] = 0.0
+
          ## Berechne die Parameter des wirklichen Alignments (but with untrained d,a,h ...)    
          #if Conf.mode == 'using_quality_scores':
          #   trueSpliceAlign, trueWeightMatch, trueWeightQuality =\
@@ -801,7 +816,6 @@ class QPalma:
       print 'Mean of pos. offset:\t\t%.2f\t%.2f\t%.2f\t%.2f' % (mean_e1Begin_neg,mean_e1End_neg,mean_e2Begin_neg,mean_e2End_neg)
 
       print 'Prediction completed'
-      self.logfh.close()
 
    def evaluatePositions(self,eBegin,eEnd):
 
@@ -968,7 +982,7 @@ def pprint_alignment(_newSpliceAlign,_newEstAlign, dna_array, est_array):
 
    return line1,line2,line3
 
-def calc_info(acc,don,exons):
+def calc_info(acc,don,exons,qualities):
    min_score = 100
    max_score = -100
 
@@ -1004,6 +1018,18 @@ def calc_info(acc,don,exons):
    mean_intron_len /= 1.0*len(exons)
    print 'Intron length max/min: %d / %d mean length %f' % (max_score,min_score,mean_intron_len)
 
+   min_score = 10000
+   max_score = 0
+   mean_quality = 0
+
+   for qVec in qualities:
+      for elem in qVec:
+         min_score = min(min_score,elem)
+         max_score = max(max_score,elem)
+
+   print 'Quality values max/min: %d / %d mean' % (max_score,min_score)
+
+
 def zero_out(vec,intervals):
    for interval in intervals:
       for pos in interval:
@@ -1014,8 +1040,8 @@ if __name__ == '__main__':
 
    mode = sys.argv[1]
 
-   if len(sys.argv) == 3 and mode == 'train':
-      qpalma.train(int(sys.argv[2]))
+   if len(sys.argv) == 2 and mode == 'train':
+      qpalma.train()
 
    elif len(sys.argv) == 3 and mode == 'predict':
       param_filename = sys.argv[2]
index 0a6c30f..e12f432 100644 (file)
@@ -52,8 +52,8 @@ const char* line_format = "%d\t%d\t%s\t%lu\t%c\t%d\t%d\t%d\t%d\t%s\t%s\t%s\n";
 
 const int read_size = 36;
 
-const int min_overlap = 6;
-const int max_overlap = 30;
+const int min_overlap = 1;
+const int max_overlap = 35;
 
 int read_nr = 1;
 
@@ -491,14 +491,24 @@ int join_reads(int exon_stop, int exon_start, Read* up_read, Read* down_read, FI
 
    } else {
 
+      //if(up_range < down_range) {
+      //   u_size = up_range;
+      //   d_size = read_size - u_size;
+      //}
+
+      //if(up_range > down_range) {
+      //   d_size = down_range;
+      //   u_size = read_size - d_size;
+      //}
+      
       if(up_range < down_range) {
-         u_size = up_range;
-         d_size = read_size - u_size;
+         d_size = down_range;
+         u_size = read_size - d_size;
       }
 
       if(up_range > down_range) {
-         d_size = down_range;
-         u_size = read_size - d_size;
+         u_size = up_range;
+         d_size = read_size - u_size;
       }
    }
 
@@ -581,7 +591,7 @@ int join_reads(int exon_stop, int exon_start, Read* up_read, Read* down_read, FI
    }
 #endif // _FDEBUG_
 
-   int status = fitting(up_read->prb,down_read->prb,u_off,d_off,u_size,d_size);
+   int status = 1; //fitting(up_read->prb,down_read->prb,u_off,d_off,u_size,d_size);
 
    if(status != 1) {
       retval = 0;