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