python/computeSpliceWeights.py

   1 #!/usr/bin/env python
   2 # -*- coding: utf-8 -*-
   3
   4 from numpy import inf
   5 from numpy.matlib import zeros
   6 import pdb
   7
   8 # 3 * 30 supporting points: x_1 ... x_30 => 3 * 30 parameters (1 .. 90): y_1 ... y_30
   9 # piecewise linear function:
  10 # take score from SVM vektor (don_supp: case 1, acc_supp: case 2) and compute length of intron: case 3
  11 # these are our values x
  12 #
  13 #        | y_1                                          if x <= x_1
  14 #        |
  15 #        |        x_i+1 - x              x - x_i
  16 # f(x) = | y_i * -----------  +  y_i+1 * -----------    if x_i <= x <= x_i+1
  17 #        |       x_i+1 - x_i             x_i+1 - x_i
  18 #        |
  19 #        | y_30                                         if x_n <= x
  20 #
  21 #  y_i and y_i+1 parameters, so the fractions are saved in the weight vectors!
  22
  23 def calculateWeights(plf, scores):
  24    currentWeight = zeros((30,1))
  25
  26    for k in range(len(scores)):
  27       value = scores[k]
  28       Lower = len([elem for elem in plf.limits if elem <= value])
  29
  30
  31       if Lower == 0:
  32          currentWeight[0] += 1
  33       elif Lower == len(plf.limits):
  34          currentWeight[-1] += 1
  35       else:
  36          # because we count from 0 in python
  37          Lower -= 1
  38          Upper = Lower+1 ; # x-werte bleiben fest
  39          #print value,Lower,Upper
  40          weightup  = 1.0*(value - plf.limits[Lower]) / (plf.limits[Upper] - plf.limits[Lower])
  41          weightlow = 1.0*(plf.limits[Upper] - value) / (plf.limits[Upper] - plf.limits[Lower])
  42          currentWeight[Upper] = currentWeight[Upper] + weightup
  43          currentWeight[Lower] = currentWeight[Lower] + weightlow
  44
  45          #print plf.limits[Lower],plf.limits[Upper]
  46          #print weightup,weightlow,currentWeight[Upper],currentWeight[Lower]
  47
  48    return currentWeight
  49
  50 def computeSpliceWeights(d, a, h, SpliceAlign, don_supp, acc_supp):
  51    ####################################################################################
  52    # 1. Donor: In don_supp stehen Werte der SVM., in SpliceAlign die Einsen
  53    ####################################################################################
  54
  55    # Picke die Positionen raus, an denen eine Donorstelle ist
  56    DonorScores = [elem for pos,elem in enumerate(don_supp) if SpliceAlign[pos] == 1]
  57    assert not ( -inf in DonorScores )
  58
  59    #print 'donor'
  60    weightDon = calculateWeights(d,DonorScores)
  61
  62    ####################################################################################
  63    # 2. Acceptor: In acc_supp stehen Werte der SVM., in SpliceAlign die Einsen
  64    ####################################################################################
  65
  66    #Den Vektor Acceptorstellen durchgehen und die Gewichtsvektoren belasten:
  67    AcceptorScores = [elem for pos,elem in enumerate(acc_supp) if pos > 0 and SpliceAlign[pos-1] == 2]
  68    assert not ( -inf in AcceptorScores )
  69
  70    #print 'acceptor'
  71    weightAcc = calculateWeights(a,AcceptorScores)
  72
  73    ####################################################################################
  74    # 3. Intron length: SpliceAlign: Gaps zaehlen und auf Gapgewichte addieren
  75    ####################################################################################
  76
  77    intron_starts = []
  78    intron_ends   = []
  79    for pos,elem in enumerate(SpliceAlign):
  80       if elem == 1:
  81          intron_starts.append(pos)
  82
  83       if elem == 2:
  84          intron_ends.append(pos)
  85
  86    assert len(intron_starts) == len(intron_ends)
  87
  88    for i in range(len(intron_starts)):
  89       assert intron_starts[i] < intron_ends[i]
  90
  91    values = [0.0]*len(intron_starts)
  92    for pos in range(len(intron_starts)):
  93       values[pos] = intron_ends[pos] - intron_starts[pos] + 1
  94
  95    weightIntron = calculateWeights(h,values)
  96
  97    return weightDon, weightAcc, weightIntron