+ fixed alignment index calculation for negative strand
[qpalma.git] / qpalma / DatasetUtils.py
index 0da50fa..b56df4a 100644 (file)
@@ -12,7 +12,7 @@ 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
 
@@ -29,37 +29,9 @@ def generateTrainingDataset(settings):
    saveData('training',dataset,settings)
 
 
-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.
-   """
-
-   #map_file = settings['map_file']
-   map_file = jp(settings['approximation_dir'],'map.vm.spliced')
+def addToDataset(map_file,dataset,settings):
    assert os.path.exists(map_file), 'Error: Can not find map file'
 
-   dataset = {}
-
    prb_offset = settings['prb_offset']
 
    # This tuple specifies an interval for valid Illumina Genome Analyzer quality values
@@ -76,7 +48,8 @@ def generatePredictionDataset(settings):
    seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
 
    for line in open(map_file):
-      if line.startswith('#'):
+      line = line.strip()
+      if line.startswith('#') or line == '':
          continue
 
       if instance_counter > 0 and instance_counter % 5000 == 0:
@@ -108,10 +81,10 @@ def generatePredictionDataset(settings):
       else:
          us_offset = pos - 1 
 
-      if pos+half_window_size < seqInfo.chromo_sizes[chromo]:
+      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
@@ -140,6 +113,42 @@ def generatePredictionDataset(settings):
 
    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)