python/qpalma.py

   1 #!/usr/bin/env python
   2 # -*- coding: utf-8 -*-
   3
   4 ###########################################################
   5 #
   6 # This file containts the
   7 #
   8 ###########################################################
   9
  10 import sys
  11 import subprocess
  12 import scipy.io
  13 import pdb
  14
  15 from numpy.matlib import mat,zeros,ones,inf
  16 from numpy.linalg import norm
  17
  18 import QPalmaDP
  19
  20 from SIQP_CPX import SIQPSolver
  21
  22 from paths_load_data import *
  23 from computeSpliceWeights import *
  24 from set_param_palma import *
  25 from computeSpliceAlign import *
  26 from penalty_lookup_new import *
  27 from compute_donacc import *
  28 from TrainingParam import Param
  29
  30 import Configuration
  31
  32 class QPalma:
  33    """
  34    A training method for the QPalma project
  35    """
  36
  37    def __init__(self):
  38       self.ARGS = Param()
  39
  40       self.logfh = open('qpalma.log','w+')
  41
  42       gen_file= '%s/genome.config' % self.ARGS.basedir
  43
  44       cmd = ['']*4
  45       cmd[0] = 'addpath /fml/ag-raetsch/home/fabio/svn/tools/utils'
  46       cmd[1] = 'addpath /fml/ag-raetsch/home/fabio/svn/tools/genomes'
  47       cmd[2] = 'genome_info = init_genome(\'%s\')' % gen_file
  48       cmd[3] = 'save genome_info.mat genome_info'
  49       full_cmd = "matlab -nojvm -nodisplay -r \"%s; %s; %s; %s; exit\"" % (cmd[0],cmd[1],cmd[2],cmd[3])
  50
  51       obj = subprocess.Popen(full_cmd,shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
  52       out,err = obj.communicate()
  53       assert err == '', 'An error occured!\n%s'%err
  54
  55       ginfo = scipy.io.loadmat('genome_info.mat')
  56       self.genome_info = ginfo['genome_info']
  57
  58       self.plog('genome_info.basedir is %s\n'%self.genome_info.basedir)
  59
  60       self.C=1.0
  61
  62       self.numDonSuppPoints     = 30
  63       self.numAccSuppPoints     = 30
  64       self.numLengthSuppPoints  = 30
  65       self.sizeMMatrix          = 36
  66       self.numQualSuppPoints    = 16*0
  67
  68       self.numFeatures = self.numDonSuppPoints + self.numAccSuppPoints + self.numLengthSuppPoints\
  69                   + self.sizeMMatrix + self.numQualSuppPoints
  70
  71       self.plog('Initializing problem...\n')
  72
  73
  74    def plog(self,string):
  75       self.logfh.write(string)
  76
  77
  78    def run(self):
  79       # Load the whole dataset
  80       Sequences, Acceptors, Donors, Exons, Ests, Noises = paths_load_data('training',self.genome_info,self.ARGS)
  81
  82       #Sequences, Acceptors, Donors, Exons, Ests, QualityScores = paths_load_data('training',self.genome_info,self.ARGS)
  83
  84       # number of training instances
  85       N = len(Sequences)
  86       self.numExamples = N
  87       assert N == len(Acceptors) and N == len(Acceptors) and N == len(Exons)\
  88       and N == len(Ests), 'The Seq,Accept,Donor,.. arrays are of different lengths'
  89
  90       self.plog('Number of training examples: %d\n'% N)
  91
  92       iteration_steps = 200 ; #upper bound on iteration steps
  93
  94       remove_duplicate_scores = False
  95       print_matrix = False
  96       anzpath = 2
  97
  98       # Initialize parameter vector
  99       # param = numpy.matlib.rand(126,1)
 100       param = Configuration.fixedParam
 101
 102       # Set the parameters such as limits penalties for the Plifs
 103       [h,d,a,mmatrix] = set_param_palma(param,self.ARGS.train_with_intronlengthinformation)
 104
 105       # delete splicesite-score-information
 106       if not self.ARGS.train_with_splicesitescoreinformation:
 107          for i in range(len(Acceptors)):
 108             if Acceptors[i] > -20:
 109                Acceptors[i] = 1
 110             if Donors[i] >-20:
 111                Donors[i] = 1
 112
 113       # Initialize solver
 114       if not __debug__:
 115          solver = SIQPSolver(self.numFeatures,self.numExamples,self.C,self.logfh)
 116
 117       # stores the number of alignments done for each example (best path, second-best path etc.)
 118       num_path = [anzpath]*N
 119       # stores the gap for each example
 120       gap      = [0.0]*N
 121
 122       qualityMatrix = zeros((self.numQualSuppPoints,1))
 123
 124       #############################################################################################
 125       # Training
 126       #############################################################################################
 127       self.plog('Starting training...\n')
 128
 129       iteration_nr = 1
 130
 131       while True:
 132          print 'Iteration step: %d'% iteration_nr
 133
 134          for exampleIdx in range(self.numExamples):
 135             if exampleIdx% 1000 == 0:
 136                print 'Current example nr %d' % exampleIdx
 137
 138             dna = Sequences[exampleIdx]
 139             est = Ests[exampleIdx]
 140
 141             exons = Exons[exampleIdx]
 142             # NoiseMatrix = Noises[exampleIdx]
 143             don_supp = Donors[exampleIdx]
 144             acc_supp = Acceptors[exampleIdx]
 145
 146             # Berechne die Parameter des wirklichen Alignments (but with untrained d,a,h ...)
 147             trueSpliceAlign, trueWeightMatch = computeSpliceAlign(dna, exons)
 148
 149             # Calculate the weights
 150             trueWeightDon, trueWeightAcc, trueWeightIntron = computeSpliceWeights(d, a, h, trueSpliceAlign, don_supp, acc_supp)
 151             trueWeight = numpy.vstack([trueWeightIntron, trueWeightDon, trueWeightAcc, trueWeightMatch, qualityMatrix ])
 152
 153             currentPhi = zeros((self.numFeatures,1))
 154             currentPhi[0:30]     = mat(d.penalties[:]).reshape(30,1)
 155             currentPhi[30:60]    = mat(a.penalties[:]).reshape(30,1)
 156             currentPhi[60:90]    = mat(h.penalties[:]).reshape(30,1)
 157             currentPhi[90:126]   = mmatrix[:]
 158             currentPhi[126:]   = qualityMatrix[:]
 159
 160             # Calculate w'phi(x,y) the total score of the alignment
 161             trueAlignmentScore = (trueWeight.T * currentPhi)[0,0]
 162
 163             # The allWeights vector is supposed to store the weight parameter
 164             # of the true alignment as well as the weight parameters of the
 165             # num_path[exampleIdx] other alignments
 166             allWeights = zeros((self.numFeatures,num_path[exampleIdx]+1))
 167             allWeights[:,0] = trueWeight[:,0]
 168
 169             AlignmentScores = [0.0]*(num_path[exampleIdx]+1)
 170             AlignmentScores[0] = trueAlignmentScore
 171
 172             ################## Calculate wrong alignment(s) ######################
 173
 174             # Compute donor, acceptor with penalty_lookup_new
 175             # returns two double lists
 176             donor, acceptor = compute_donacc(don_supp, acc_supp, d, a)
 177
 178             #myalign wants the acceptor site on the g of the ag
 179             acceptor = acceptor[1:]
 180             acceptor.append(-inf)
 181
 182             dna = str(dna)
 183             est = str(est)
 184             dna_len = len(dna)
 185             est_len = len(est)
 186             ps = h.convert2SWIG()
 187
 188             matchmatrix = QPalmaDP.createDoubleArrayFromList(mmatrix.flatten().tolist()[0])
 189             mm_len = 36
 190
 191             d_len = len(donor)
 192             donor = QPalmaDP.createDoubleArrayFromList(donor)
 193             a_len = len(acceptor)
 194             acceptor = QPalmaDP.createDoubleArrayFromList(acceptor)
 195
 196             currentAlignment = QPalmaDP.Alignment()
 197             qualityMat = QPalmaDP.createDoubleArrayFromList(qualityMatrix)
 198             currentAlignment.setQualityMatrix(qualityMat,self.numQualSuppPoints)
 199
 200             print 'PYTHON: Calling myalign...'
 201             # calculates SpliceAlign, EstAlign, weightMatch, Gesamtscores, dnaest
 202             currentAlignment.myalign( num_path[exampleIdx], dna, dna_len,\
 203             est, est_len, ps, matchmatrix, mm_len, donor, d_len,\
 204             acceptor, a_len, remove_duplicate_scores, print_matrix)
 205             print 'PYTHON: After myalign call...'
 206
 207             newSpliceAlign       = QPalmaDP.createIntArrayFromList([0]*(dna_len*num_path[exampleIdx]))
 208             newEstAlign          = QPalmaDP.createIntArrayFromList([0]*(est_len*num_path[exampleIdx]))
 209             newWeightMatch       = QPalmaDP.createIntArrayFromList([0]*(mm_len*num_path[exampleIdx]))
 210             newAlignmentScores   = QPalmaDP.createDoubleArrayFromList([.0]*num_path[exampleIdx])
 211
 212             pdb.set_trace()
 213
 214             currentAlignment.getAlignmentResults(newSpliceAlign, newEstAlign, newWeightMatch, newAlignmentScores)
 215
 216             spliceAlign = zeros((num_path[exampleIdx]*dna_len,1))
 217             weightMatch = zeros((num_path[exampleIdx]*mm_len,1))
 218
 219             print 'spliceAlign'
 220             for i in range(dna_len*num_path[exampleIdx]):
 221                spliceAlign[i] = newSpliceAlign[i]
 222                print '%f' % (spliceAlign[i])
 223
 224             print 'weightMatch'
 225             for i in range(mm_len*num_path[exampleIdx]):
 226                weightMatch[i] = newWeightMatch[i]
 227                print '%f' % (weightMatch[i])
 228
 229             for i in range(num_path[exampleIdx]):
 230                AlignmentScores[i+1] = newAlignmentScores[i]
 231
 232             print AlignmentScores
 233
 234             spliceAlign = spliceAlign.reshape(num_path[exampleIdx],dna_len)
 235             weightMatch = weightMatch.reshape(num_path[exampleIdx],mm_len)
 236             # Calculate weights of the respective alignments Note that we are
 237             # calculating n-best alignments without any hamming loss, so we
 238             # have to keep track which of the n-best alignments correspond to
 239             # the true one in order not to incorporate a true alignment in the
 240             # constraints. To keep track of the true and false alignments we
 241             # define an array true_map with a boolean indicating the
 242             # equivalence to the true alignment for each decoded alignment.
 243             true_map = [0]*(num_path[exampleIdx]+1)
 244             true_map[0] = 1
 245             path_loss = [0]*(num_path[exampleIdx]+1)
 246
 247             for pathNr in range(num_path[exampleIdx]):
 248                #dna_numbers = dnaest{1,pathNr}
 249                #est_numbers = dnaest{2,pathNr}
 250
 251                weightDon, weightAcc, weightIntron = computeSpliceWeights(d, a, h, spliceAlign[pathNr,:].flatten().tolist()[0], don_supp, acc_supp)
 252
 253                # sum up positionwise loss between alignments
 254                for alignPosIdx in range(len(spliceAlign[pathNr,:])):
 255                   if spliceAlign[pathNr,alignPosIdx] != trueSpliceAlign[alignPosIdx]:
 256                      path_loss[pathNr+1] += 1
 257
 258                # Gewichte in restliche Zeilen der Matrix speichern
 259                wp = numpy.vstack([weightIntron, weightDon, weightAcc, weightMatch[pathNr,:].T, qualityMatrix ])
 260                allWeights[:,pathNr+1] = wp
 261
 262                hpen = mat(h.penalties).reshape(len(h.penalties),1)
 263                dpen = mat(d.penalties).reshape(len(d.penalties),1)
 264                apen = mat(a.penalties).reshape(len(a.penalties),1)
 265
 266                features = numpy.vstack([hpen , dpen , apen , mmatrix[:]])
 267                AlignmentScores[pathNr+1] = (allWeights[:,pathNr+1].T * features)[0,0]
 268
 269                # Check wether scalar product + loss equals viterbi score
 270                #assert math.fabs(newAlignmentScores[pathNr] - AlignmentScores[pathNr+1]) < 1e-6,\
 271                #'Scalar prod + loss is not equal Viterbi score. Respective values are %f, %f' % \
 272                #(newAlignmentScores[pathNr],AlignmentScores[pathNr+1])
 273
 274                #  # if the pathNr-best alignment is very close to the true alignment consider it as true
 275                if norm( allWeights[:,0] - allWeights[:,pathNr+1] ) < 1e-5:
 276                   true_map[pathNr+1] = 1
 277
 278                # the true label sequence should not have a larger score than the maximal one WHYYYYY?
 279
 280                # this means that all n-best paths are to close to each other
 281                # we have to extend the n-best search to a (n+1)-best
 282                if len([elem for elem in true_map if elem == 1]) == len(true_map):
 283                   num_path[exampleIdx] = num_path[exampleIdx]+1
 284
 285                # Choose true and first false alignment for extending A
 286                firstFalseIdx = -1
 287                for map_idx,elem in enumerate(true_map):
 288                   if elem == 0:
 289                      firstFalseIdx = map_idx
 290                      break
 291
 292                # if there is at least one useful false alignment add the
 293                # corresponding constraints to the optimization problem
 294                if firstFalseIdx != -1:
 295                   trueWeights        = allWeights[:,0]
 296                   firstFalseWeights  = allWeights[:,firstFalseIdx]
 297
 298                   # LMM.py code:
 299                   deltas  = firstFalseWeights - trueWeights
 300                   if not __debug__:
 301                      const_added = solver.addConstraint(deltas, exampleIdx)
 302                      objValue,w,self.slacks = solver.solve()
 303
 304                      sum_xis = 0
 305                      for elem in self.slacks:
 306                         sum_xis +=  elem
 307
 308             if exampleIdx==0:
 309                break
 310
 311          iteration_nr += 1
 312          break
 313
 314       self.logfh.close()
 315       print 'Training completed'
 316
 317 if __name__ == '__main__':
 318    qpalma = QPalma()
 319    qpalma.run()