+ update makefiles to fetch automatically valid Python includes and libs
[qpalma.git] / qpalma / qpalma_main.py
index b3349be..daf1975 100644 (file)
@@ -12,14 +12,12 @@ import os.path
 import pdb
 import sys
 
-from qpalma.sequence_utils import SeqSpliceInfo,DataAccessWrapper,unbracket_seq
-
 import numpy
 from numpy.matlib import mat,zeros,ones,inf
 from numpy.linalg import norm
 
 #from qpalma.SIQP_CPX import SIQPSolver
-#from qpalma.SIQP_CVXOPT import SIQPSolver
+from qpalma.SIQP_CVXOPT import SIQPSolver
 
 import QPalmaDP
 import qpalma
@@ -29,43 +27,35 @@ from qpalma.computeSpliceAlignWithQuality import *
 from qpalma.TrainingParam import Param
 from qpalma.Plif import Plf,compute_donacc
 
+from qpalma.sequence_utils import SeqSpliceInfo,DataAccessWrapper,unbracket_seq
 from qpalma.utils import pprint_alignment, get_alignment
 
+jp = os.path.join
+dotProd = lambda x,y: (x.T * y)[0,0]
+
+
 class SpliceSiteException:
    pass
 
 
-def getData(training_set,exampleKey,run):
-   """ This function...  """
+def preprocessExample(training_set,exampleKey,seqInfo,settings):
+   """ 
+   This function...  
+   """
 
-   currentSeqInfo,original_est,currentQualities,currentExons = training_set[exampleKey]
+   currentSeqInfo,originalRead,currentQualities,currentExons = training_set[exampleKey]
    id,chr,strand,up_cut,down_cut = currentSeqInfo
 
-   est = original_est
-   est = "".join(est)
-   est = est.lower()
-   est = unbracket_est(est)
-   est = est.replace('-','')
-
-   assert len(est) == run['read_size'], pdb.set_trace()
-   est_len = len(est)
-
-   #original_est = OriginalEsts[exampleIdx]
-   original_est = "".join(original_est)
-   original_est = original_est.lower()
-
-   dna_flat_files =  '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
-   dna, acc_supp, don_supp = get_seq_and_scores(chr,strand,up_cut,down_cut,dna_flat_files)
+   read = unbracket_seq(originalRead)
+   read = read.replace('-','')
 
-   #currentDNASeq, currentAcc, currentDon = seqInfo.get_seq_and_scores(chromo,strand,genomicSeq_start,genomicSeq_stop)
+   dna, acc_supp, don_supp = seqInfo.get_seq_and_scores(chr,strand,up_cut,down_cut)
 
-   original_exons = currentExons
-   exons = original_exons - (up_cut-1)
-   exons[0,0] -= 1
-   exons[1,0] -= 1
+   exons = currentExons - up_cut
 
    if exons.shape == (2,2):
       fetched_dna_subseq = dna[exons[0,0]:exons[0,1]] + dna[exons[1,0]:exons[1,1]]
+      print fetched_dna_subseq
      
       donor_elem = dna[exons[0,1]:exons[0,1]+2]
       acceptor_elem = dna[exons[1,0]-2:exons[1,0]]
@@ -78,9 +68,53 @@ def getData(training_set,exampleKey,run):
          print 'invalid acceptor in example %d'% exampleKey
          raise SpliceSiteException
 
-      assert len(fetched_dna_subseq) == len(est), pdb.set_trace()
+      assert len(fetched_dna_subseq) == len(read), pdb.set_trace()
+
+   return dna,read,acc_supp,don_supp,exons,originalRead,currentQualities
+
+
+def performAlignment(dna,read,quality,mmatrix,donor,acceptor,ps,qualityPlifs,current_num_path,prediction_mode,settings):
+   """
+   Given the needed input this method calls the QPalma C module which
+   calculates a dynamic programming in order to obtain an alignment
+   """
+
+   prb = QPalmaDP.createDoubleArrayFromList(quality)
+   chastity = QPalmaDP.createDoubleArrayFromList([.0]*len(read))
+
+   matchmatrix = QPalmaDP.createDoubleArrayFromList(mmatrix.flatten().tolist()[0])
+   mm_len = settings['matchmatrixRows']*settings['matchmatrixCols']
+
+   don_len = len(donor)
+   donor = QPalmaDP.createDoubleArrayFromList(donor)
+   acc_len = len(acceptor)
+   acceptor = QPalmaDP.createDoubleArrayFromList(acceptor)
+
+   # Create the alignment object representing the interface to the C/C++ code.
+   currentAlignment = QPalmaDP.Alignment(settings['numQualPlifs'],settings['numQualSuppPoints'], settings['enable_quality_scores'])
+
+   c_qualityPlifs = QPalmaDP.createPenaltyArrayFromList([elem.convert2SWIG() for elem in qualityPlifs])
+
+   # calculates SpliceAlign, EstAlign, weightMatch, Gesamtscores, dnaest
+   currentAlignment.myalign( current_num_path, dna, len(dna),\
+    read, len(read), prb, chastity, ps, matchmatrix, mm_len, donor, don_len,\
+    acceptor, acc_len, c_qualityPlifs, settings['remove_duplicate_scores'],\
+    settings['print_matrix'] )
+
+   if prediction_mode:
+      # part that is only needed for prediction
+      result_len = currentAlignment.getResultLength()
+      dna_array,read_array = currentAlignment.getAlignmentArraysNew()
+   else:
+      dna_array = None
+      read_array = None
+
+   newSpliceAlign, newEstAlign, newWeightMatch, newDPScores, newQualityPlifsFeatures =\
+   currentAlignment.getAlignmentResultsNew()
+
+   return newSpliceAlign, newEstAlign, newWeightMatch, newDPScores,\
+   newQualityPlifsFeatures, dna_array, read_array
 
-   return dna,est,acc_supp,don_supp,exons,original_est,currentQualities
 
 
 class QPalma:
@@ -89,64 +123,23 @@ class QPalma:
    the alignment.
    """
    
-   def __init__(self,run,seqInfo,dmode=False):
+   def __init__(self,seqInfo,dmode=False):
       self.ARGS = Param()
       self.qpalma_debug_mode = dmode
-      self.run = run
       self.seqInfo = seqInfo
+      self.logfh = None
 
 
    def plog(self,string):
-      self.logfh.write(string)
-      self.logfh.flush()
-
-
-   def do_alignment(self,dna,est,quality,mmatrix,donor,acceptor,ps,qualityPlifs,current_num_path,prediction_mode):
-      """
-      Given the needed input this method calls the QPalma C module which
-      calculates a dynamic programming in order to obtain an alignment
-      """
-
-      dna_len = len(dna)
-      est_len = len(est)
-
-      prb = QPalmaDP.createDoubleArrayFromList(quality)
-      chastity = QPalmaDP.createDoubleArrayFromList([.0]*est_len)
-
-      matchmatrix = QPalmaDP.createDoubleArrayFromList(mmatrix.flatten().tolist()[0])
-      mm_len = self.run['matchmatrixRows']*self.run['matchmatrixCols']
-
-      d_len = len(donor)
-      donor = QPalmaDP.createDoubleArrayFromList(donor)
-      a_len = len(acceptor)
-      acceptor = QPalmaDP.createDoubleArrayFromList(acceptor)
-
-      # Create the alignment object representing the interface to the C/C++ code.
-      currentAlignment = QPalmaDP.Alignment(self.run['numQualPlifs'],self.run['numQualSuppPoints'], self.use_quality_scores)
-      c_qualityPlifs = QPalmaDP.createPenaltyArrayFromList([elem.convert2SWIG() for elem in qualityPlifs])
-      # calculates SpliceAlign, EstAlign, weightMatch, Gesamtscores, dnaest
-      currentAlignment.myalign( current_num_path, dna, dna_len,\
-       est, est_len, prb, chastity, ps, matchmatrix, mm_len, donor, d_len,\
-       acceptor, a_len, c_qualityPlifs, self.run['remove_duplicate_scores'],\
-       self.run['print_matrix'] )
-
-      if prediction_mode:
-         # part that is only needed for prediction
-         result_len = currentAlignment.getResultLength()
-         dna_array,est_array = currentAlignment.getAlignmentArraysNew()
+      if self.logfh != None:
+         self.logfh.write(string)
+         self.logfh.flush()
       else:
-         dna_array = None
-         est_array = None
-
-      newSpliceAlign, newEstAlign, newWeightMatch, newDPScores, newQualityPlifsFeatures =\
-      currentAlignment.getAlignmentResultsNew()
-
-      return newSpliceAlign, newEstAlign, newWeightMatch, newDPScores,\
-      newQualityPlifsFeatures, dna_array, est_array
+         print string
 
 
-   def init_train(self,training_set):
-      full_working_path = os.path.join(self.run['alignment_dir'],self.run['name'])
+   def init_training(self,dataset_fn,settings,set_name):
+      full_working_path = jp(settings['training_dir'],set_name)
 
       #assert not os.path.exists(full_working_path)
       if not os.path.exists(full_working_path):
@@ -158,37 +151,13 @@ class QPalma:
       os.chdir(full_working_path)
 
       self.logfh = open('_qpalma_train.log','w+')
-      cPickle.dump(self.run,open('run_obj.pickle','w+'))
-
-      self.plog("Settings are:\n")
-      self.plog("%s\n"%str(self.run))
-
-      if self.run['mode'] == 'normal':
-         self.use_quality_scores = False
-
-      elif self.run['mode'] == 'using_quality_scores':
-         self.use_quality_scores = True
-      else:
-         assert(False)
-
-
-   def setUpSolver(self):
-      # Initialize solver 
-      self.plog('Initializing problem...\n')
-      
-      try:
-         solver = SIQPSolver(run['numFeatures'],numExamples,run['C'],self.logfh,run)
-      except:
-         self.plog('Got no license. Telling queue to reschedule job...\n')
-         sys.exit(99)
 
-      solver.enforceMonotonicity(lengthSP,lengthSP+donSP)
-      solver.enforceMonotonicity(lengthSP+donSP,lengthSP+donSP+accSP)
+      training_set = cPickle.load(open(dataset_fn))
 
-      return solver
+      self.train(training_set,settings,set_name)
 
 
-   def train(self,training_set):
+   def train(self,training_set,settings,set_name):
       """
       The mainloop for training.
       """
@@ -199,40 +168,52 @@ class QPalma:
       self.noImprovementCtr = 0
       self.oldObjValue = 1e8
 
-      iteration_steps         = self.run['iter_steps']
-      remove_duplicate_scores = self.run['remove_duplicate_scores']
-      print_matrix            = self.run['print_matrix']
-      anzpath                 = self.run['anzpath']
+      remove_duplicate_scores = settings['remove_duplicate_scores']
+      print_matrix            = settings['print_matrix']
+
+      lengthSP    = settings['numLengthSuppPoints']
+      donSP       = settings['numDonSuppPoints']
+      accSP       = settings['numAccSuppPoints']
+      mmatrixSP   = settings['matchmatrixRows']*settings['matchmatrixCols']
+      numq        = settings['numQualSuppPoints']
+      totalQualSP = settings['totalQualSuppPoints']
+
+      # calculate the total number of features
+      numFeatures = settings['numFeatures']
 
       # Initialize parameter vector
-      param = numpy.matlib.rand(run['numFeatures'],1)
-   
-      lengthSP    = self.run['numLengthSuppPoints']
-      donSP       = self.run['numDonSuppPoints']
-      accSP       = self.run['numAccSuppPoints']
-      mmatrixSP   = self.run['matchmatrixRows']*run['matchmatrixCols']
-      numq        = self.run['numQualSuppPoints']
-      totalQualSP = self.run['totalQualSuppPoints']
+      param = numpy.matlib.rand(numFeatures,1)
+
+      # we take the first quality vector of the tuple of quality vectors
+      quality_index = 0
 
       # no intron length model
-      if not self.run['enable_intron_length']:
+      if not settings['enable_intron_length']:
          param[:lengthSP] *= 0.0
 
       # Set the parameters such as limits penalties for the Plifs
-      [h,d,a,mmatrix,qualityPlifs] = set_param_palma(param,self.ARGS.train_with_intronlengthinformation,self.run)
+      [h,d,a,mmatrix,qualityPlifs] = set_param_palma(param,self.ARGS.train_with_intronlengthinformation,settings)
+
+      # Initialize solver 
+      self.plog('Initializing problem...\n')
+      
+      try:
+         solver = SIQPSolver(numFeatures,numExamples,settings['C'],self.logfh,settings)
+      except:
+         self.plog('Got no license. Telling queue to reschedule job...\n')
+         self.plog(sys.exc_info())
+         sys.exit(99)
 
-      solver = self.setUpSolver()
+      #solver.enforceMonotonicity(lengthSP,lengthSP+donSP)
+      #solver.enforceMonotonicity(lengthSP+donSP,lengthSP+donSP+accSP)
 
       # stores the number of alignments done for each example (best path, second-best path etc.)
-      num_path = [anzpath]*numExamples
+      num_path = [settings['anzpath']]*numExamples
 
-      # stores the gap for each example
-      gap      = [0.0]*numExamples
-
-      currentPhi = zeros((run['numFeatures'],1))
+      currentPhi = zeros((numFeatures,1))
       totalQualityPenalties = zeros((totalQualSP,1))
 
-      numConstPerRound = run['numConstraintsPerRound']
+      numConstPerRound = settings['numConstraintsPerRound']
       solver_call_ctr = 0
 
       suboptimal_example = 0
@@ -240,30 +221,30 @@ class QPalma:
       param_idx = 0
       const_added_ctr = 0
 
-      featureVectors = zeros((run['numFeatures'],numExamples))
+      featureVectors = zeros((numFeatures,numExamples))
 
       self.plog('Starting training...\n')
       # the main training loop
       while True:
-         if iteration_nr == iteration_steps:
+         if iteration_nr == settings['iter_steps']:
             break
 
+         #print 'Got %d examples' % len(training_set.keys())
          for exampleIdx,example_key in enumerate(training_set.keys()):
             print 'Current example %d' % example_key
+
             try:
-               dna,est,acc_supp,don_supp,exons,original_est,currentQualities =\
-               getData(training_set,example_key,run)
+               dna,read,acc_supp,don_supp,exons,original_read,currentQualities =\
+               preprocessExample(training_set,example_key,self.seqInfo,settings)
             except SpliceSiteException:
                continue
 
-            dna_len = len(dna)
-
-            if run['enable_quality_scores']:
+            if settings['enable_quality_scores']:
                quality = currentQualities[quality_index]
             else:
                quality = [40]*len(read)
 
-            if not run['enable_splice_signals']:
+            if not settings['enable_splice_scores']:
                for idx,elem in enumerate(don_supp):
                   if elem != -inf:
                      don_supp[idx] = 0.0
@@ -273,15 +254,16 @@ class QPalma:
                      acc_supp[idx] = 0.0
 
             # Berechne die Parameter des wirklichen Alignments (but with untrained d,a,h ...)    
-            if run['mode'] == 'using_quality_scores':
+            if settings['enable_quality_scores']:
                trueSpliceAlign, trueWeightMatch, trueWeightQuality ,dna_calc =\
-               computeSpliceAlignWithQuality(dna, exons, est, original_est,\
-               quality, qualityPlifs,run)
+               computeSpliceAlignWithQuality(dna, exons, read, original_read,\
+               quality, qualityPlifs,settings)
             else:
                trueSpliceAlign, trueWeightMatch, trueWeightQuality = computeSpliceAlignWithQuality(dna, exons)
 
             dna_calc = dna_calc.replace('-','')
 
+
             #print 'right before computeSpliceWeights exampleIdx %d' % exampleIdx
             # Calculate the weights
             trueWeightDon, trueWeightAcc, trueWeightIntron =\
@@ -293,7 +275,7 @@ class QPalma:
             currentPhi[lengthSP+donSP:lengthSP+donSP+accSP]                   = mat(a.penalties[:]).reshape(accSP,1)
             currentPhi[lengthSP+donSP+accSP:lengthSP+donSP+accSP+mmatrixSP]   = mmatrix[:]
 
-            if run['mode'] == 'using_quality_scores':
+            if settings['enable_quality_scores']:
                totalQualityPenalties = param[-totalQualSP:]
                currentPhi[lengthSP+donSP+accSP+mmatrixSP:]                    = totalQualityPenalties[:]
 
@@ -303,7 +285,7 @@ class QPalma:
             # The allWeights vector is supposed to store the weight parameter
             # of the true alignment as well as the weight parameters of the
             # num_path[exampleIdx] other alignments
-            allWeights = zeros((run['numFeatures'],num_path[exampleIdx]+1))
+            allWeights = zeros((numFeatures,num_path[exampleIdx]+1))
             allWeights[:,0] = trueWeight[:,0]
 
             AlignmentScores = [0.0]*(num_path[exampleIdx]+1)
@@ -316,15 +298,16 @@ class QPalma:
 
             ps = h.convert2SWIG()
 
-            newSpliceAlign, newEstAlign, newWeightMatch, newDPScores,\
-            newQualityPlifsFeatures, unneeded1, unneeded2 =\
-            self.do_alignment(dna,est,quality,mmatrix,donor,acceptor,ps,qualityPlifs,num_path[exampleIdx],False)
-            mm_len = run['matchmatrixRows']*run['matchmatrixCols']
+            newSpliceAlign, newEstAlign, newWeightMatch, newDPScores, newQualityPlifsFeatures, unused1, unused2 =\
+            performAlignment(dna,read,quality,mmatrix,donor,acceptor,ps,qualityPlifs,num_path[exampleIdx],False,settings)
 
-            newSpliceAlign = newSpliceAlign.reshape(num_path[exampleIdx],dna_len)
-            newWeightMatch = newWeightMatch.reshape(num_path[exampleIdx],mm_len)
+            mm_len = settings['matchmatrixRows']*settings['matchmatrixCols']
 
-            newQualityPlifsFeatures = newQualityPlifsFeatures.reshape(num_path[exampleIdx],run['totalQualSuppPoints'])
+            # check for correct reshaping !!
+            newSpliceAlign = numpy.mat(newSpliceAlign).reshape(num_path[exampleIdx],len(dna))
+            newWeightMatch = numpy.mat(newWeightMatch).reshape(num_path[exampleIdx],mm_len)
+
+            newQualityPlifsFeatures = numpy.mat(newQualityPlifsFeatures).reshape(num_path[exampleIdx],settings['totalQualSuppPoints'])
             # Calculate weights of the respective alignments. Note that we are calculating n-best alignments without 
             # hamming loss, so we have to keep track which of the n-best alignments correspond to the true one in order 
             # not to incorporate a true alignment in the
@@ -339,7 +322,7 @@ class QPalma:
                h, newSpliceAlign[pathNr,:].flatten().tolist()[0], don_supp,\
                acc_supp)
               
-               decodedQualityFeatures = zeros((run['totalQualSuppPoints'],1))
+               decodedQualityFeatures = zeros((settings['totalQualSuppPoints'],1))
                decodedQualityFeatures = newQualityPlifsFeatures[pathNr,:].T
                # Gewichte in restliche Zeilen der Matrix speichern
                allWeights[:,pathNr+1] = numpy.vstack([weightIntron, weightDon, weightAcc, newWeightMatch[pathNr,:].T, decodedQualityFeatures[:]])
@@ -358,12 +341,12 @@ class QPalma:
                   distinct_scores = True
 
                # Check wether scalar product + loss equals viterbi score
-               if not math.fabs(newDPScores[pathNr,0] - AlignmentScores[pathNr+1]) <= 1e-5:
+               if not math.fabs(newDPScores[pathNr] - AlignmentScores[pathNr+1]) <= 1e-5:
                   self.plog("Scalar prod. + loss not equals Viterbi output!\n")
                   pdb.set_trace()
 
                self.plog(" scalar prod (correct) : %f\n"%AlignmentScores[0])
-               self.plog(" scalar prod (pred.) : %f %f\n"%(newDPScores[pathNr,0],AlignmentScores[pathNr+1]))
+               self.plog(" scalar prod (pred.) : %f %f\n"%(newDPScores[pathNr],AlignmentScores[pathNr+1]))
 
                # if the pathNr-best alignment is very close to the true alignment consider it as true
                if norm( allWeights[:,0] - allWeights[:,pathNr+1] ) < 1e-5:
@@ -393,11 +376,7 @@ class QPalma:
 
                if False:
                   self.plog("Is considered as: %d\n" % true_map[1])
-
-                  #result_len = currentAlignment.getResultLength()
-
                   dna_array,est_array = currentAlignment.getAlignmentArraysNew()
-
                   _newSpliceAlign = newSpliceAlign[0].flatten().tolist()[0]
                   _newEstAlign = newEstAlign[0].flatten().tolist()[0]
 
@@ -413,8 +392,9 @@ class QPalma:
 
                # end of one example processing 
 
-            # call solver every nth example //added constraint
-            if exampleIdx != 0 and exampleIdx % numConstPerRound == 0:
+            # call solver every nth example / added constraint
+            #if exampleIdx != 0 and exampleIdx % numConstPerRound == 0:
+            if True:
                objValue,w,self.slacks = solver.solve()
                solver_call_ctr += 1
 
@@ -444,29 +424,27 @@ class QPalma:
                cPickle.dump(param,open('param_%d.pickle'%param_idx,'w+'))
                param_idx += 1
                [h,d,a,mmatrix,qualityPlifs] =\
-               set_param_palma(param,self.ARGS.train_with_intronlengthinformation,self.run)
-
-         ##############################################
-         # end of one iteration through all examples  #
-         ##############################################
+               set_param_palma(param,self.ARGS.train_with_intronlengthinformation,settings)
 
+         # end of one iteration through all examples
          self.plog("suboptimal rounds %d\n" %suboptimal_example)
 
          if self.noImprovementCtr == numExamples*2:
-            FinalizeTraining(param,'param_%d.pickle'%param_idx)
+            self.FinalizeTraining(param,'param_%d.pickle'%param_idx)
 
          iteration_nr += 1
 
       #
       # end of optimization 
       #  
-      FinalizeTraining(param,'param_%d.pickle'%param_idx)
+      self.FinalizeTraining(param,'param_%d.pickle'%param_idx)
 
 
-   def FinalizeTraining(self,vector,name):
+   def FinalizeTraining(self,param,name):
       self.plog("Training completed")
       cPickle.dump(param,open(name,'w+'))
-      self.logfh.close()
+      if self.logfh != None:
+         self.logfh.close()
    
 
 ###############################################################################
@@ -478,14 +456,13 @@ class QPalma:
 ###############################################################################
 
 
-   def init_prediction(self,dataset_fn,prediction_keys,param_fn,set_name):
+   def init_prediction(self,dataset_fn,prediction_keys,settings,set_name):
       """
       Performing a prediction takes...
       """
       self.set_name = set_name
 
-      #full_working_path = os.path.join(run['alignment_dir'],run['name'])
-      full_working_path = self.run['result_dir']
+      full_working_path = settings['prediction_dir']
 
       print 'full_working_path is %s' % full_working_path 
 
@@ -513,28 +490,20 @@ class QPalma:
       # we do not need the full dataset anymore
       del dataset
    
-      # Load parameter vector to predict with
-      param = cPickle.load(open(param_fn))
+      self.predict(prediction_set,settings)
 
-      self.predict(prediction_set,param)
 
-
-   def predict(self,prediction_set,param):
+   def predict(self,prediction_set,settings):
       """
       This method...
       """
 
-      if self.run['mode'] == 'normal':
-         self.use_quality_scores = False
-
-      elif self.run['mode'] == 'using_quality_scores':
-         self.use_quality_scores = True
-      else:
-         assert(False)
+      # Load parameter vector to predict with
+      param = cPickle.load(open(settings['prediction_param_fn']))
 
       # Set the parameters such as limits/penalties for the Plifs
       [h,d,a,mmatrix,qualityPlifs] =\
-      set_param_palma(param,self.ARGS.train_with_intronlengthinformation,self.run)
+      set_param_palma(param,self.ARGS.train_with_intronlengthinformation,settings)
 
       if not self.qpalma_debug_mode:
          self.plog('Starting prediction...\n')
@@ -560,7 +529,7 @@ class QPalma:
             if not self.qpalma_debug_mode:
                self.plog('Loading example id: %d...\n'% int(id))
 
-            if self.run['enable_quality_scores']:
+            if settings['enable_quality_scores']:
                quality = currentQualities[quality_index]
             else:
                quality = [40]*len(read)
@@ -569,9 +538,10 @@ class QPalma:
                currentDNASeq, currentAcc, currentDon = self.seqInfo.get_seq_and_scores(chromo,strand,genomicSeq_start,genomicSeq_stop)
             except:
                self.problem_ctr += 1
+               print sys.exc_info()
                continue
 
-            if not self.run['enable_splice_signals']:
+            if not settings['enable_splice_scores']:
                for idx,elem in enumerate(currentDon):
                   if elem != -inf:
                      currentDon[idx] = 0.0
@@ -581,7 +551,7 @@ class QPalma:
                      currentAcc[idx] = 0.0
 
             current_prediction = self.calc_alignment(currentDNASeq, read,\
-            quality, currentDon, currentAcc, d, a, h, mmatrix, qualityPlifs)
+            quality, currentDon, currentAcc, d, a, h, mmatrix, qualityPlifs,settings)
 
             current_prediction['id']         = id
             current_prediction['chr']        = chromo
@@ -606,12 +576,11 @@ class QPalma:
       self.logfh.close()
 
 
-   def calc_alignment(self, dna, read, quality, don_supp, acc_supp, d, a, h, mmatrix, qualityPlifs):
+   def calc_alignment(self, dna, read, quality, don_supp, acc_supp, d, a, h, mmatrix, qualityPlifs,settings):
       """
       Given two sequences and the parameters we calculate on alignment
       """
 
-      run = self.run
       donor, acceptor = compute_donacc(don_supp, acc_supp, d, a)
 
       if '-' in read:
@@ -623,9 +592,9 @@ class QPalma:
 
       newSpliceAlign, newEstAlign, newWeightMatch, newDPScores,\
       newQualityPlifsFeatures, dna_array, read_array =\
-      self.do_alignment(dna,read,quality,mmatrix,donor,acceptor,ps,qualityPlifs,1,True)
+      performAlignment(dna,read,quality,mmatrix,donor,acceptor,ps,qualityPlifs,1,True,settings)
 
-      mm_len = run['matchmatrixRows']*run['matchmatrixCols']
+      mm_len = settings['matchmatrixRows']*settings['matchmatrixCols']
 
       true_map    = [0]*2
       true_map[0] = 1