+ added basic scripts to generate runs for an experiment
[qpalma.git] / scripts / compile_dataset.py
index 0b6b0bf..cf79622 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 == 20000:
          break
 
    print 'Full dataset has size %d' % len(Sequences)
@@ -167,12 +175,15 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
    chr            = fRead['chr']
    strand         = fRead['strand']
    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:
@@ -304,15 +320,19 @@ def process_read(reReads,fRead,currentGene,dna_flat_files,test):
       #assert dna_fragment_2.replace('-','') == dna_annot_2, pdb.set_trace()
 
       if not dna_fragment_1.replace('-','') == dna_annot_1:
-         print 'genomic seq mismatch'
+         #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 'genomic seq mismatch'
+         #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