+ got rid of some legacy code
authorFabio <fabio@congo.fml.local>
Wed, 8 Oct 2008 09:29:54 +0000 (11:29 +0200)
committerFabio <fabio@congo.fml.local>
Wed, 8 Oct 2008 09:29:54 +0000 (11:29 +0200)
+ added test code
+ fixed offsets in alignment file calculation
+ added license information to files with missing info

21 files changed:
DynProg/QPalmaDP.i
README
qpalma/Configuration.py [deleted file]
qpalma/DatasetUtils.py
qpalma/Lookup.py
qpalma/OutputFormat.py
qpalma/PipelineHeuristic.py
qpalma/Run.py [deleted file]
qpalma/SIQP.py
qpalma/SIQP_CPX.py
qpalma/SIQP_CVXOPT.py
qpalma/computeSpliceAlignWithQuality.py
qpalma/computeSpliceWeights.py
qpalma/gridtools.py
qpalma/parsers.py [deleted file]
qpalma/sequence_utils.py
scripts/qpalma_pipeline.py
tests/.test_sequence_utils.py.swp [deleted file]
tests/test_data/chr1.dna.flat [deleted file]
tests/test_qpalma_prediction.py
tests/test_sequence_utils.py

index 5eedd23..47fa728 100644 (file)
 %include "penalty_info.h"
 %include "qpalma_dp.h"
 
-/*
 %array_class(int, intArray) ;
 %array_class(double, doubleArray) ;
-%array_functions(double, doubleFArray) ;
 %array_functions(struct penalty_struct, penaltyArray) ;
+
+/*
+%array_functions(double, doubleFArray) ;
 %array_class(Pre_score, Pre_scoreArray) ;
+*/
 
 %pythoncode
 %{
@@ -57,6 +59,4 @@ def createListFromDoubleArray(array, array_len):
    for i in range(array_len):
       list[i] = doubleArray_getitem(array,i)
    return list
-
 %} 
-*/
diff --git a/README b/README
index 20cbdbd..9fec71f 100644 (file)
--- a/README
+++ b/README
@@ -5,14 +5,44 @@ source, so the header files and libraries for Python must be installed,
 as well as the core binaries.
 
 
-Instructions
-------------
+Installation Instructions
+-------------------------
 
-Type
+You need the following packages
+
+- SWIG  ( >=1.3.31 )
+- numpy ( >=1.0.4 )
+
+- Pythongrid
+- Genefinding
+
+You can obtain SWIG from http://www.swig.org and numpy from http://numpy.scipy.org/
+or ask your sysadmin to install the packages.
+
+Pythongrid and Genefinding are two utility packages from our lab. These can be obtained via:
+   
+   http://
+
+Once you have the above packages installed you can start installing QPalma.
+
+Assuming you downloaded QPalma-1.0.tar.gz type:
+
+   tar -xzvf QPalma-1.0.tar.gz
+
+   cd QPalma-1.0
 
    python setup.py build
 
-to build the necessary C++ extension modules and copy the QPalma python code.
+to build the necessary C++ extension modules and the python code.
+
+
+Setting environment variables
+-----------------------------
+
+You will need the following variables set:
+
+- SGE_ROOT should point to 
+- 
 
 
 
@@ -26,9 +56,7 @@ Troubleshooting
 
 
 
-HELP!
------
-
-If however everything failes contact:
+Contact
+-------
 
-fabio@tuebingen.mpg.de
+If however everything failes contact: fabio@tuebingen.mpg.de
diff --git a/qpalma/Configuration.py b/qpalma/Configuration.py
deleted file mode 100644 (file)
index e1cb4c7..0000000
+++ /dev/null
@@ -1,192 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import os.path
-import cPickle
-
-#
-# choose a path where all results of the QPalma pipeline will be stored
-#
-
-result_dir = '/fml/ag-raetsch/home/fabio/tmp/transcriptome_data/prediction'
-
-
-###############################################################################
-# Load a random but fixed initial parameter vector this makes debugging easier
-###############################################################################
-
-#fixedParamQ = cPickle.load(open('/fml/ag-raetsch/home/fabio/svn/projects/QPalma/scripts/randInitParam.pickle'))
-
-###############################################################################
-#
-# The parameters for the QPalma algorithm
-#
-#min_intron_len = 20
-#max_intron_len = 2000
-#
-#min_svm_score = 0.0 
-#max_svm_score = 1.0
-#
-#min_qual = -5
-#max_qual = 40
-
-numConstraintsPerRound = 50
-
-###############################################################################
-# 
-# CHOOSING THE MODE 
-# 
-# 'normal' means work like Palma 'using_quality_scores' means work like Palma
-# plus using sequencing quality scores
-#
-###############################################################################
-
-#mode = 'normal'
-mode = 'using_quality_scores'
-
-###############################################################################
-# 
-# When using quality scores our scoring function is defined as
-#
-# f: S_e x R x S -> R, where S_e is {A,C,G,T,N} and S = {A,C,G,T,N,-}
-# 
-# as opposed to a usage without quality scores when we only have
-# 
-# f: S x S -> R 
-#
-# The matrix of plifs is defined as follows:
-#
-#  elem |   -  a  c  g  t  n
-#  -------------------------  
-#  idx  |   0  1  2  3  4  5
-#   
-#              dna
-#
-#           -  a  c  g  t  n
-#        a  
-# est    c
-#        g
-#        t
-#        n
-#  
-# so the index can be calculated as (estnum-1)*6 + dnanum.
-# Ests do not have gaps with quality scores so we look up the matchmatrix
-# instead.
-###############################################################################
-
-read_size = 38
-#read_size = 36
-extension = (250,500)
-
-numLengthSuppPoints  = 10
-numDonSuppPoints     = 10
-numAccSuppPoints     = 10
-numQualSuppPoints    = 10
-
-if mode == 'normal':
-   sizeMatchmatrix   = (6,6)
-   estPlifs = 0
-   dnaPlifs = 0
-   numQualPlifs = estPlifs*dnaPlifs
-elif mode == 'using_quality_scores':
-   sizeMatchmatrix   = (6,1)
-   estPlifs = 5
-   dnaPlifs = 6
-   numQualPlifs = estPlifs*dnaPlifs
-else:
-   assert False, 'Wrong operation mode specified'
-
-###############################################################################
-#
-# GENERAL SETTINGS CONCERNING THE SOLVER
-#
-###############################################################################
-
-iter_steps = 40
-remove_duplicate_scores = False
-print_matrix            = False
-anzpath                 = 2
-
-if mode == 'normal':
-   #fixedParam = fixedParamQ
-   fixedParam = None
-elif mode == 'using_quality_scores':
-   fixedParam = None
-else:
-   assert False, 'Wrong operation mode specified'
-
-###############################################################################
-#
-# DATA SETTINGS CONCERNING THE SPLITS AND FILE LOCATIONS
-#
-###############################################################################
-training_begin    =  0
-training_end      =  10000
-
-prediction_begin  =  10000
-prediction_end    =  40000
-
-joinp = os.path.join
-
-tmp_dir              = '/fml/ag-raetsch/home/fabio/tmp/solexa_tmp'
-data_path            = '/fml/ag-raetsch/share/projects/qpalma/solexa'  
-
-original_path  = joinp(data_path,'original_solexa_data')
-annot_path     = joinp(data_path,'annotation_data')
-remapped_path  = joinp(data_path,'remapped_solexa_data')
-
-dna_flat_fn    =  '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
-
-gff_fn         = joinp(annot_path,'TAIR7_GFF3_genes_Chr%s.gff_v1')
-
-###############################################################################
-#
-# SANITY CHECKS
-#
-###############################################################################
-#assert numQualPlifs       >= 0
-#assert numDonSuppPoints    > 1
-#assert numAccSuppPoints    > 1
-#assert numLengthSuppPoints > 1 
-#assert numQualSuppPoints   > 1
-assert os.path.exists(dna_flat_fn), 'DNA data does not exist!'
-
-extended_alphabet = ['-','a','c','g','t','n','[',']']
-alphabet          = ['-','a','c','g','t','n']
-
-###############################################################################
-#
-# Settings for the VMatch pipeline steps
-#
-###############################################################################
-
-reads_location = ''
-
-#
-# First VMatch step 
-#
-
-mismatches_1      = 2
-end_gap_1         = 0 
-repeat_mapping_1  = 1
-seedlength_1      = 9
-suffixtree_1      = '/media/oka_raid/nobackup/data/Vmatch/ATH/ATH1.v5.seed9.fa'
-
-#
-# Second VMatch step 
-#
-
-mismatches_2      = 1
-sub_mismatches_2  = 1
-min_short_end_2   = 14
-repeat_mapping_2  = 1
-seedlength_2      = 9
-suffixtree_2      = '/media/oka_raid/nobackup/data/Vmatch/ATH/ATH1.v5.seed9.fa'
-
-#
-#
-# do not modify anything below this line
-#
-#
-
-conf_object_path = os.path.join(result_dir,'config_object.pickle')
index 0da50fa..981455c 100644 (file)
@@ -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
@@ -140,6 +112,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)
 
 
index f22ecd7..da9365c 100644 (file)
@@ -1,5 +1,10 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# Written (W) 2008 Fabio De Bona
+# Copyright (C) 2008 Max-Planck-Society
 
 import sys
 import os
index af7178c..d9ba9b0 100644 (file)
@@ -109,8 +109,12 @@ def recalculatePositions(chromo,start_pos,strand,starts,ends,seqInfo,window_size
       total_size = seqInfo.chromo_sizes[chromo]
       start_pos = total_size - start_pos
 
-   starts   = [int(elem) + start_pos for elem in starts]
-   ends     = [int(elem) + start_pos for elem in ends]
+   if strand == '+':
+      starts   = [int(elem) + start_pos-2 for elem in starts]
+      ends     = [int(elem) + start_pos-2 for elem in ends]
+   else:
+      starts   = [int(elem) + start_pos for elem in starts]
+      ends     = [int(elem) + start_pos for elem in ends]
 
    if strand == '-':
       starts   = [e-window_size for e in starts]
@@ -167,7 +171,30 @@ def createAlignmentOutput(settings,chunk_fn,result_fn):
 def createBlatOutput(current_prediction,settings):
    """
    This function takes one prediction as input and returns the corresponding
-   alignment in BLAT format.
+   alignment in BLAT format:
+
+    1. matches     - Number of bases that match that aren't repeats
+    2. misMatches  - Number of bases that don't match
+    3. repMatches  - Number of bases that match but are part of repeats
+    4. nCount      - Number of 'N' bases
+    5. qNumInsert  - Number of inserts in query
+    6. qBaseInsert - Number of bases inserted in query
+    7. tNumInsert  - Number of inserts in target
+    8. tBaseInsert - Number of bases inserted in target
+    9. strand      - '+' or '-' for query strand. In mouse, second '+'or '-' is for genomic strand
+   10. qName       - Query sequence name
+   11. qSize       - Query sequence size
+   12. qStart      - Alignment start position in query
+   13. qEnd        - Alignment end position in query
+   14. tName       - Target sequence name
+   15. tSize       - Target sequence size
+   16. tStart      - Alignment start position in target
+   17. tEnd        - Alignment end position in target
+   18. blockCount  - Number of blocks in the alignment
+   19. blockSizes  - Comma-separated list of sizes of each block
+   20. qStarts     - Comma-separated list of starting positions of each block in query
+   21. tStarts     - Comma-separated list of starting positions of each block in target 
+
    """
 
    # fetch the data needed
index 3fd1418..847652a 100644 (file)
@@ -110,9 +110,7 @@ class PipelineHeuristic:
       self.mmatrix = mmatrix
       self.qualityPlifs = qualityPlifs
 
-      self.read_size    = 38
-      self.prb_offset   = 64
-      #self.prb_offset   = 50
+      self.prb_offset   = settings['prb_offset']
 
       # parameters of the heuristics to decide whether the read is spliced
       #self.splice_thresh      = 0.005
diff --git a/qpalma/Run.py b/qpalma/Run.py
deleted file mode 100644 (file)
index a9755bd..0000000
+++ /dev/null
@@ -1,30 +0,0 @@
-#!/usr/bin/env python
-# -*- coding:latin-1 -*-
-
-
-class Run(dict):
-   """
-   This class represents the needed parameters for one single run of an
-   algorithm. This means that if we for example perform a model selection over
-   C then for each value of C we have to create one Run object storing the
-   respective parameters.
-   """
-
-   def __init__(self):
-      pass
-
-   def __repr__(self):
-      
-      result = ""
-      for key,val in self.iteritems():
-         result += "%s:\t\t %s\n" % (key,str(val))
-
-      return result
-
-if __name__ == '__main__':
-   r1 = Run()
-   r1['attrib_1'] = 12
-   r1['attrib_2'] = 22
-
-   print "%s" % r1
-
index fecfb1e..c994785 100644 (file)
@@ -1,5 +1,10 @@
-#!/usr/bin/env python
-# -*- coding:latin-1 -*-
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# Written (W) 2008 Fabio De Bona
+# Copyright (C) 2008 Max-Planck-Society
 
 import sys
 import os
index 3ca654a..3bfbf1a 100644 (file)
@@ -1,5 +1,10 @@
-#!/usr/bin/env python
-# -*- coding:latin-1 -*-
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# Written (W) 2008 Fabio De Bona
+# Copyright (C) 2008 Max-Planck-Society
 
 import sys
 import os
index c09340b..d3be2e5 100644 (file)
@@ -1,5 +1,10 @@
-#!/usr/bin/env python
-# -*- coding:latin-1 -*-
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# Written (W) 2008 Fabio De Bona
+# Copyright (C) 2008 Max-Planck-Society
 
 import sys
 import os
index 243bafa..e498098 100644 (file)
@@ -1,5 +1,10 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# Written (W) 2008 Fabio De Bona
+# Copyright (C) 2008 Max-Planck-Society
 
 import numpy
 from numpy.matlib import mat,zeros,ones,inf
index 60c65d6..c4a18ce 100644 (file)
@@ -1,5 +1,10 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# Written (W) 2008 Fabio De Bona
+# Copyright (C) 2008 Max-Planck-Society
 
 from numpy import inf
 from numpy.matlib import zeros
index 83a5da8..5d29743 100644 (file)
@@ -137,7 +137,7 @@ class ApproximationTask(ClusterTask):
 
       for idx in range(0,num_splits):
          data_fname     = jp(result_dir,'map.part_%d'%idx)
-         result_fname   = jp(result_dir,'map.vm.part_%d.heuristic.spliced'%idx)
+         result_fname   = jp(result_dir,'map.vm.part_%d'%idx)
          self.result_files.append(result_fname)
 
          current_job = KybJob(gridtools.ApproximationTaskStarter,[run_fname,data_fname,param_fname,result_fname,self.settings])
@@ -152,6 +152,7 @@ class ApproximationTask(ClusterTask):
    def collectResults(self):
       result_dir  = self.settings['approximation_dir']
       combined_fn = jp(result_dir,'map.vm.spliced')
+      self.result_files = map(lambda x:x+'.spliced',self.result_files)
       combine_files(self.result_files,combined_fn)
       combine_files([combined_fn,self.settings['spliced_reads_fn']],'map.vm')
 
diff --git a/qpalma/parsers.py b/qpalma/parsers.py
deleted file mode 100644 (file)
index d7d1576..0000000
+++ /dev/null
@@ -1,422 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import pdb
-import os
-import sys
-import mmap
-import resource
-
-class ReadParser:
-   """
-   A base class for the Solexa reads parsers.
-   """
-
-   def __init__(self,filename):
-      self.fh = open(filename)
-
-   def __iter__(self):
-      return self
-
-   def next(self):
-      pass
-
-   def parse(self):
-      lines = []
-
-      for elem in self:
-         lines.append(elem)
-         
-      return lines
-
-
-def parse(filename):
-   entries = []
-   
-   fh = open(filename)
-
-   all_lines = fh.read()
-   all_lines = all_lines.split('\n')
-
-   for parsed_line in all_lines:
-      if parsed_line == '':
-         continue
-
-      id,chr,pos,strand,mismatches,length,offset,seq,prb,cal_prb,chastity =\
-      parsed_line.split()
-      chr   = int(chr)
-      pos   = int(pos)
-
-      mismatches  = int(mismatches)
-      length      = int(length)
-      offset      = int(offset)
-
-      line_d = {'id':id, 'chr':chr, 'pos':pos, 'strand':strand,\
-      'mismatches':mismatches, 'length':length, 'offset':offset,\
-      'seq':seq,'prb':prb,'cal_prb':cal_prb,'chastity':chastity}
-
-      entries.append((line_d,parsed_line))
-
-   return entries
-
-
-   def parse(self):
-      entries = {}
-
-      for elem in self.fh:
-         line_d = self.parseLine(elem)
-         id = line_d['id']
-         assert id not in entries, pdb.set_trace()
-         entries[id] = line_d
-         
-      return entries
-
-class RemappedReadParser(ReadParser):
-   """
-   This class offers a parser for the reads that are remapped by the vmatch
-   utility. 
-
-   According to the docu the entries are:
-   
-   ID, Chromosome, Position, Orientation (D or P), Mismatches, Alignment length, Offset, Alignment
-
-   """
-
-   def __init__(self,filename):
-      ReadParser.__init__(self,filename)
-
-   def parseLine(self,line):
-      """
-      We assume that a line has the following entires:
-
-      """
-      id,chr1,pos1,seq1,chr2,pos2,seq2,quality = line.split()
-
-      chr1 = int(chr1)
-      chr2 = int(chr2)
-
-      seq1=seq1.lower()
-      seq2=seq2.lower()
-
-      pos1 = int(pos1)
-      pos2 = int(pos2)
-
-      line_d = {'id':id, 'chr1':chr1, 'pos1':pos1, 'seq1':seq1, 'chr2':chr2,\
-      'pos2':pos2, 'seq2':seq2, 'quality':'quality'}
-
-      return line_d
-
-   def next(self):
-      for line in self.fh:
-         line = line.strip()
-         yield self.parseLine(line)
-
-      raise StopIteration
-
-   def parse(self):
-      entries = {}
-
-      for elem in self.fh:
-         line_d = self.parseLine(elem)
-         id = line_d['id']
-         try:
-            entries[id] = [line_d]
-         except:
-            old_entry = entries[id]
-            old_entry.append(line_d)
-            entries[id] = old_entry
-
-      return entries
-
-
-class MapParser(ReadParser):
-   """
-   This class offers a parser for the reads that are remapped by the vmatch
-   utility. 
-
-   According to the docu the entries are:
-   
-   ID, Chromosome, Position, Orientation (D or P), Mismatches, Alignment length, Offset, Alignment
-
-   """
-
-   def __init__(self,filename):
-      ReadParser.__init__(self,filename)
-
-   def parseLine(self,line):
-      """
-      We assume that a line has the following entries:
-
-      """
-
-      id,chr,pos,strand,mismatches,length,offset,seq = line.split()
-
-      chr   = int(chr)
-      pos   = int(pos)
-
-      if strand == 'D':
-         strand = '+'
-
-      if strand == 'P':
-         strand = '-'
-
-      mismatches  = int(mismatches)
-      length      = int(length)
-      offset      = int(offset)
-
-      seq=seq.lower()
-
-      line_d = {'id':id, 'chr':chr, 'pos':pos, 'strand':strand,\
-      'mismatches':mismatches, 'length':length, 'offset':offset,\
-      'seq':seq}
-
-      return line_d
-
-   def next(self):
-      for line in self.fh:
-         line = line.strip()
-         yield self.parseLine(line)
-
-      raise StopIteration
-
-   def parse(self):
-      entries = {}
-
-      for elem in self.fh:
-         line_d = self.parseLine(elem)
-         id = line_d['id']
-         try:
-            entries[id] = [line_d]
-         except:
-            old_entry = entries[id]
-            old_entry.append(line_d)
-            entries[id] = old_entry
-
-      return entries
-
-
-
-def parse_map_vm_heuristic(filename):
-   """
-   This function parse the map.vm file for the PipelineHeuristics
-   """
-   entries = []
-   
-   fh = open(filename)
-
-   all_lines = fh.read()
-   all_lines = all_lines.split('\n')
-
-   for parsed_line in all_lines:
-      if parsed_line == '':
-         continue
-
-      id,chr,pos,strand,mismatches,length,offset,seq,prb,cal_prb,chastity =\
-      parsed_line.split()
-      chr   = int(chr)
-      pos   = int(pos)
-
-      mismatches  = int(mismatches)
-      length      = int(length)
-      offset      = int(offset)
-
-      line_d = {'id':id, 'chr':chr, 'pos':pos, 'strand':strand,\
-      'mismatches':mismatches, 'length':length, 'offset':offset,\
-      'seq':seq,'prb':prb,'cal_prb':cal_prb,'chastity':chastity}
-
-      entries.append((line_d,parsed_line))
-
-   return entries
-
-def cpu():
-   return (resource.getrusage(resource.RUSAGE_SELF).ru_utime+\
-   resource.getrusage(resource.RUSAGE_SELF).ru_stime) 
-
-class Line:
-   pass
-
-def parse_filtered_reads(filename):
-   """
-   
-   """
-   entries = {}
-   
-   #fh = open(filename)
-   #all_lines = fh.read()
-   #all_lines = all_lines.split('\n')
-
-   #for parsed_line in all_lines:
-   #   if parsed_line == '':
-   #      continue
-
-   print 'asking for map'
-   data = map_file(filename)
-   print 'obtained map'
-
-   strand_map = ['-','+']
-
-   full_line_split_time = 0  
-   full_array_time = 0
-   full_dict_time = 0
-
-   parsing_start = cpu()
-
-   print 'start parsing'
-   
-   while True:
-
-      start = cpu()
-      parsed_line = data.readline().strip()
-
-      if parsed_line == '':
-         break
-
-      id,chr,strand,seq,splitpos,read_size,_prb,_cal_prb,_chastity,gene_id,p_start,exon_stop,exon_start,p_stop,true_cut = parsed_line.split()
-      stop = cpu()
-
-      full_line_split_time += stop-start
-   
-      id = int(id)
-      splitpos = int(splitpos)
-      read_size = int(read_size)
-
-      seq=seq.lower()
-
-      #assert strand in ['D','P']
-
-      strand = strand_map[strand == 'D']
-
-      chr = int(chr)
-   
-      start_array_time = cpu()
-
-      prb      = [0]*read_size
-      cal_prb  = [0]*read_size
-      chastity = [0]*read_size
-
-      for idx in range(read_size):
-         prb[idx]       = ord(_prb[idx])-50
-         cal_prb[idx]   = ord(_cal_prb[idx])-64
-         chastity[idx]  = ord(_chastity[idx])+10
-
-      #import array
-
-      #prb      = array.array('b',_prb)
-      #cal_prb  = array.array('b',_cal_prb)
-      #chastity = array.array('b',_chastity)
-
-      #for idx in range(read_size):
-      #   prb[idx]       -= 50
-      #   cal_prb[idx]   -= 64
-      #   chastity[idx]  += 10
-
-      #assert prb == __prb
-      #assert cal_prb == __cal_prb
-      #assert chastity == __chastity
-
-      stop_array_time = cpu()
-
-      full_array_time += stop_array_time - start_array_time 
-
-
-      start_dict_time  = cpu()
-
-      p_start = int(p_start)
-      exon_stop = int(exon_stop)
-      exon_start = int(exon_start)
-      p_stop = int(p_stop)
-      true_cut = int(true_cut)
-
-      line_d = {'id':id, 'chr':chr, 'strand':strand, 'seq':seq, 'splitpos':splitpos,\
-      'read_size':read_size, 'prb':prb, 'cal_prb':cal_prb, 'chastity':chastity, 'gene_id':gene_id,\
-      'p_start':p_start, 'exon_stop':exon_stop, 'exon_start':exon_start,\
-      'p_stop':p_stop,'true_cut':true_cut}
-
-      stop_dict_time  = cpu()
-
-      full_dict_time += stop_dict_time - start_dict_time 
-
-      entries[id] = line_d
-
-   parsing_stop = cpu()
-
-   print 'parsing took %f secs' % (stop-start)
-   print 'line split time %f ' % full_line_split_time 
-   print 'array time %f ' % full_array_time 
-   print 'dict time %f ' % full_dict_time 
-
-   return entries
-
-
-def parse_map_vm(filename):
-   """
-   Parses the map.vm files from Vmatch
-
-
-   For an original read we can have several vmatch locations
-   """
-   entries = {}
-   
-   #fh = open(filename)
-   #all_lines = fh.read()
-   #all_lines = all_lines.split('\n')
-
-   #for parsed_line in all_lines:
-   #   if parsed_line == '':
-   #      continue
-
-   data = map_file(filename)
-   
-   while True:
-
-      parsed_line = data.readline()
-
-      if parsed_line == '':
-         break
-
-      id,chr,pos,strand,mismatches,length,offset,seq,q1,q2,q3 = parsed_line.split()
-
-      id = int(id)
-
-      chr   = int(chr)
-      pos   = int(pos)
-
-      if strand == 'D':
-         strand = '+'
-
-      if strand == 'P':
-         strand = '-'
-
-      mismatches  = int(mismatches)
-      length      = int(length)
-      offset      = int(offset)
-
-      seq=seq.lower()
-
-      line_d = {'id':id, 'chr':chr, 'pos':pos, 'strand':strand,\
-      'mismatches':mismatches, 'length':length, 'offset':offset}
-   
-      if entries.has_key(id):
-         #pdb.set_trace()
-         old_entry = entries[id]
-         old_entry.append(line_d)
-         entries[id] = old_entry
-      else:
-         entries[id] = [line_d]
-
-   return entries
-
-
-def map_file(filename):
-   """
-   As the mmapped file object behaves just like plain file objects
-   you do not have to unmap but just close the file.
-   """
-
-   file = open(filename, "r+")
-   size = os.path.getsize(filename)
-
-   data = mmap.mmap(file.fileno(), size, mmap.ACCESS_READ)
-   assert size == len(data), 'Could not map whole file at once!'
-
-   return data
index a64226f..e111f3c 100644 (file)
@@ -25,7 +25,9 @@ import numpy
 
 from numpy.matlib import inf
 
-from Genefinding import doIntervalQuery
+from Genefinding import doIntervalQuery,new_intp,intp_assign
+
+from QPalmaDP import createIntArrayFromList
 
 
 def load_genomic(chromosome, strand, start, stop, genome, one_based=True):
index cffccc4..ac7194e 100644 (file)
@@ -99,10 +99,10 @@ class System:
       # Before creating a candidate spliced read dataset we have to first filter
       # the matches from the first seed finding run.
 
-      approx_task = ApproximationTask(self.settings)
-      approx_task.CreateJobs()
-      approx_task.Submit()
-      approx_task.CheckIfTaskFinished()
+      #approx_task = ApproximationTask(self.settings)
+      #approx_task.CreateJobs()
+      #approx_task.Submit()
+      #approx_task.CheckIfTaskFinished()
       
       # After filtering combine the filtered matches from the first run and the
       # found matches from the second run to a full dataset
diff --git a/tests/.test_sequence_utils.py.swp b/tests/.test_sequence_utils.py.swp
deleted file mode 100644 (file)
index e69de29..0000000
diff --git a/tests/test_data/chr1.dna.flat b/tests/test_data/chr1.dna.flat
deleted file mode 100644 (file)
index 869713a..0000000
+++ /dev/null
@@ -1 +0,0 @@
-acgtacgtagcgatctgctcgtacgtgctcgatcaggcagtczzuuctgt
index 3aca1f9..0e41bb2 100644 (file)
@@ -16,10 +16,12 @@ import os.path
 import pdb
 import unittest
 
-from qpalma_main import QPalma
+from qpalma.qpalma_main import QPalma
 from qpalma.utils import print_prediction
 
-from createAlignmentFileFromPrediction import alignment_reconstruct
+from qpalma.Run import Run
+
+from qpalma.OutputFormat import alignment_reconstruct
 
 from qpalma.sequence_utils import DataAccessWrapper,SeqSpliceInfo
 
index 00e8870..4c6483b 100644 (file)
@@ -185,19 +185,29 @@ class TestLookupTable(unittest.TestCase):
 
 
    def testTableThalianaData(self):
-      g_dir    = '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
+      """
+      """
+
+      g_dir    = '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome'
       acc_dir  = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/acc'
       don_dir  = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/don'
+      g_fmt    = 'chr%d.dna.flat'
+      s_fmt    = 'contig_%d%s'
 
-      g_fmt = 'chr%d.dna.flat'
-      s_fmt = 'contig_%d%s'
+      settings = {}
+      settings['genome_dir']           = g_dir
+      settings['acceptor_scores_loc']  = acc_dir
+      settings['donor_scores_loc']     = don_dir 
 
-      num_chromo = 2
+      settings['genome_file_fmt']      = g_fmt
+      settings['splice_score_file_fmt']= s_fmt
 
-      lt1 = LookupTable(g_dir,acc_dir,don_dir,g_fmt,s_fmt,range(1,num_chromo))
+      settings['allowed_fragments']    = [1]
 
-      accessWrapper = DataAccessWrapper(g_dir,acc_dir,don_dir,g_fmt,s_fmt)
-      seqInfo = SeqSpliceInfo(accessWrapper,range(1,num_chromo))
+      lt1 = LookupTable(settings)
+
+      accessWrapper = DataAccessWrapper(settings)
+      seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
 
       seq = 'CTGGCCAAAAGCTCAGGGAAGACGCAGCCTAGGGCTCC'.lower()
       seq = reverse_complement(seq)
@@ -205,10 +215,12 @@ class TestLookupTable(unittest.TestCase):
       strand = '-'
       pos = 10475515
       dna = seqInfo.get_seq_and_scores(chromo,strand,pos,pos+38,True)
+
       self.assertEqual(seq,dna)
 
-      dna = lt1.get_seq_and_scores(chromo,strand,pos,pos+38,True)
-      pdb.set_trace()
+      full_info = lt1.get_seq_and_scores(chromo,strand,pos,pos+38)
+      dna = full_info[0]
+
       self.assertEqual(seq,dna)
 
       dna,acc,don = seqInfo.get_seq_and_scores(1,'+',1,1369)