+ set proper logging in optimizer code
[qpalma.git] / qpalma / DatasetUtils.py
index b98de77..6c41062 100644 (file)
 import array
 import cPickle
 import os
+import os.path
 import pdb
 
-from sequence_utils import SeqSpliceInfo,DataAccessWrapper,get_flat_file_size,reverse_complement,unbracket_seq,create_bracket_seq,reconstruct_dna_seq
+from sequence_utils import SeqSpliceInfo,DataAccessWrapper,get_flat_file_size,unbracket_seq,create_bracket_seq,reconstruct_dna_seq
+
+jp = os.path.join
 
 illumina_ga_range = (-5,40)
 #roche_454_range   =
 
 
+def processQuality(raw_qualities,prb_offset,quality_interval,perform_checks):
+   """
+   In order to save some space we use a signed char to store the
+   qualities. Each quality element can range as follows: -128 <= elem <= 127
+   """
+
+   q_values = map(lambda x: ord(x)-prb_offset,raw_qualities)
+
+   if perform_checks:
+      for entry in q_values:
+         assert quality_interval[0] <= entry <= quality_interval[1], 'Error: Your read has invalid quality values: %d' % entry
+
+   return array.array('b',q_values)
+
+
+def checkExons(dna,exons,readAlignment,exampleKey):
+   """
+    
+   """
+
+   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'% exampleKey
+      return False
+
+   if not ( acceptor_elem == 'ag' ):
+      print 'invalid acceptor in example %d'% exampleKey
+      return False
+
+   read = unbracket_seq(readAlignment)
+   read = read.replace('-','')
+
+   assert len(fetched_dna_subseq) == len(read), pdb.set_trace()
+
+   return True
+
+
 def generateTrainingDataset(settings):
    """
    This function creates a training dataset.
+
+
+   Create lower case original ("bracketed") reads
+
    """
+
    dataset = []
 
-   saveData('training',dataset,settings)
+   half_window_size = settings['half_window_size']
 
+   instance_counter = 0
 
-def generatePredictionDataset(settings):
-   """
-   This function is responsible for the prediction dataset generation for the matches of
-   the second vmatch run and the heuristic-filtered matches of the first run.
+   accessWrapper = DataAccessWrapper(settings)
+   seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
 
-   An example in our dataset consists of:
+   for line in open(map_file):
+      line = line.strip()
+      if line.startswith('#') or line == '':
+         continue
 
-   - information on the original read:
-         * id
-         * nucleotide sequence
-         * quality vectors (at the moment: prb and chastity)
+      if instance_counter > 0 and instance_counter % 5000 == 0:
+         print 'processed %d examples' % instance_counter
 
-   - information on the target dna sequence:
-         * chromosome
-         * strand
-         * fragment start 
-         * fragment end positions
+      slist    = line.split()
 
-   this information is stored in tuples which are stored then in a dictionary
-   indexed by the reads unique id.
+      id       = int(slist[0])
+      chromo   = int(slist[1])
 
-   CAUTION: The dictionary stores a list of the above defined tuples under each
-   'id'. This is because one read can be found several times on the genome.
-   """
+      # for the time being we only consider chromosomes 1-5
+      if not chromo in settings['allowed_fragments']:
+         continue
+
+      # convert D/P strand info to +/-
+      strand = ['-','+'][slist[2] == 'D']
+
+      seqBeginning   = int(slist[3])
+      seqEnd         = int(slist[4])
+
+      readAlignment  = slist[5]
+      prb = processQuality(slist[6],prb_offset,quality_intervals,perform_checks)
+      
+      exons = numpy.mat([0,0,0,0]).reshape((2,2))
+
+      exons[0,0]     = int(slist[7])
+      exons[0,1]     = int(slist[8])
+      exons[1,0]     = int(slist[9])
+      exons[1,1]     = int(slist[10])
+
+      dna = seqInfo.get_seq_and_scores(chr,strand,up_cut,down_cut,only_seq=True)
+      relative_exons = exons - up_cut
+
+      assert checkExons(dna,relative_exons,readAlignment,id)
+
+      currentSeqInfo = (id,chromo)
+      
+      dataset.setdefault(id, []).append((currentSeqInfo,readAlignment,[prb],exons))
+
+   saveData('training',dataset,settings)
 
-   #map_file = settings['map_file']
-   map_file = jp(result_dir  = self.global_settings['approximation_dir'],'map.vm')
+
+def addToDataset(map_file,dataset,settings):
    assert os.path.exists(map_file), 'Error: Can not find map file'
 
-   dataset = {}
+   perform_checks = settings['perform_checks']
 
-   prb_offset = 64
+   prb_offset = settings['prb_offset']
 
    # This tuple specifies an interval for valid Illumina Genome Analyzer quality values
    if settings['platform'] == 'IGA':
@@ -73,6 +145,10 @@ def generatePredictionDataset(settings):
    seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
 
    for line in open(map_file):
+      line = line.strip()
+      if line.startswith('#') or line == '':
+         continue
+
       if instance_counter > 0 and instance_counter % 5000 == 0:
          print 'processed %d examples' % instance_counter
 
@@ -97,55 +173,87 @@ def generatePredictionDataset(settings):
       assert not '[' in read_seq and not ']' in read_seq
 
       # we use an area of +/-  `self.half_window_size` nucleotides around the seed position
-      if pos > self.half_window_size+1:
-         us_offset = self.half_window_size
+      if pos > half_window_size+1:
+         us_offset = half_window_size
       else:
          us_offset = pos - 1 
 
-      if pos+self.half_window_size < seqInfo.chromo_sizes[chromo]:
-         ds_offset = self.half_window_size
+      if pos+half_window_size < seqInfo.getFragmentSize(chromo):
+         ds_offset = half_window_size
       else:
-         ds_offset = seqInfo.chromo_sizes[chromo]-pos-1
+         ds_offset = seqInfo.getFragmentSize(chromo)-pos-1
          
       genomicSeq_start = pos - us_offset
       genomicSeq_stop  = pos + ds_offset
 
       dna_seq = seqInfo.get_seq_and_scores(chromo,strand,genomicSeq_start,genomicSeq_stop,True)
       
-      # In order to save some space we use a signed char to store the
-      # qualities. Each quality element can range as follows: -128 <= elem <= 127
-      
-      q_values = map(lambda x: ord(x)-self.prb_offset,slist[5])
-
-      if settings['perform_checks']:
-         for entry in q_values:
-            assert quality_interval[0] <= entry <= quality_interval[1], 'Error: Your read has invalid quality values: %d' % entry
-
-      prb = array.array('b',q_values)
-
+      prb = processQuality(slist[5],prb_offset,quality_interval,perform_checks)
+   
       currentSeqInfo = (id,chromo,strand,genomicSeq_start,genomicSeq_stop)
-      currentQualities = [prb]
 
       # As one single read can have several vmatch matches we store all these
       # matches under the unique id of the read
-      dataset.setdefault(id, []).append((currentSeqInfo,read_seq,currentQualities))
+      dataset.setdefault(id, []).append((currentSeqInfo,read_seq,[prb]))
 
       instance_counter += 1
 
    print 'Total examples processed: %d' % instance_counter
 
+   return dataset
+
+
+def generatePredictionDataset(settings):
+   """
+   This function is responsible for the prediction dataset generation for the matches of
+   the second vmatch run and the heuristic-filtered matches of the first run.
+
+   An example in our dataset consists of:
+
+   - information on the original read:
+         * id
+         * nucleotide sequence
+         * quality vectors (at the moment: prb and chastity)
+
+   - information on the target dna sequence:
+         * chromosome
+         * strand
+         * fragment start 
+         * fragment end positions
+
+   this information is stored in tuples which are stored then in a dictionary
+   indexed by the reads unique id.
+
+   CAUTION: The dictionary stores a list of the above defined tuples under each
+   'id'. This is because one read can be found several times on the genome.
+   """
+
+   dataset = {}
+
+   map_file = settings['spliced_reads_fn']
+   dataset = addToDataset(map_file,dataset,settings)
+
+   map_file = jp(settings['approximation_dir'],'map.vm.spliced')
+   dataset = addToDataset(map_file,dataset,settings)
+
    saveData('prediction',dataset,settings)
 
 
 def saveData(prefix,dataset,settings):
    """
    """
-   ddir = settings['dataset_dir']
-   dataset_fn        = jp(ddir,'%s_data.pickle'%prefix)
-   dataset_keys_fn   = jp(ddir,'%s_data.keys.pickle'%prefix)
 
-   assert not os.path.exists(dataset_fn), 'The data_file already exists!'
-   assert not os.path.exists(dataset_keys_fn), 'The data_keys file already exists!'
+   if prefix == 'prediction':
+      dataset_fn        = settings['prediction_dataset_fn']
+      dataset_keys_fn   = settings['prediction_dataset_keys_fn']
+   elif prefix == 'training':
+      dataset_fn        = settings['training_dataset_fn']
+      dataset_keys_fn   = settings['training_dataset_keys_fn']
+   else:
+      assert False
+
+   #assert not os.path.exists(dataset_fn), 'The data file already exists!'
+   #assert not os.path.exists(dataset_keys_fn), 'The data keys file already exists!'
 
    # saving new dataset and single keys as well
    cPickle.dump(dataset,open(dataset_fn,'w+'),protocol=2)