numberOfExons = (exons.shape)[0] # how many rows ?
exonSizes = [-1]*numberOfExons
+ assert numberOfExons == 3
+
for idx in range(numberOfExons):
exonSizes[idx] = exons[idx,1] - exons[idx,0]
# 1:donorpos, 3:intron 2:acceptorpos, 0:exon, 4: dangling end
SpliceAlign = []
- if exons[0,0] > 1:
+ if exons[0,0] > 0:
SpliceAlign.extend([4]*(exons[0,0]))
- #for i=1:nr_exons
- # exon_length = exon_sizes(i); %exons(i,2) is begin of intron
- # SpliceAlign = [SpliceAlign, zeros(1,exon_length)] ;
- # if i~= nr_exons,
- # intron_length = exons(i+1,1) - exons(i,2) ;
- # SpliceAlign = [SpliceAlign, 1, ones(1,intron_length-2)*3, 2] ;
- # end
- #end
-
for idx in range(numberOfExons):
exonLength = exonSizes[idx]
SpliceAlign.extend([0]*exonLength)
if idx < numberOfExons-1:
intronLength = exons[idx+1,0] - exons[idx,1]
- SpliceAlign.extend([1])
- SpliceAlign.extend([3]*((intronLength-2)))
- SpliceAlign.extend([2])
-
+ SpliceAlign.extend([1]+[3]*(intronLength-2)+[2])
+
if len(dna) > exons[2,1]:
- SpliceAlign.extend([4]*(len(dna)+1-exons[2,1]))
+ SpliceAlign.extend([4]*(len(dna)-exons[2,1]))
assert len(SpliceAlign) == len(dna), pdb.set_trace()
Lower = len([elem for elem in plf.limits if elem <= value])
Upper = Lower+1 ; # x-werte bleiben fest
- Lower = Lower - 1
- Upper = Upper - 1
+ #print Lower,Upper,len(plf.limits)
if Lower == 0:
currentWeight[0] += 1
currentWeight[Upper] = currentWeight[Upper] + weightup
currentWeight[Lower] = currentWeight[Lower] + weightlow
- #print value,plf.limits[Lower],plf.limits[Upper]
- #print weightup,weightlow,currentWeight[Upper],currentWeight[Lower]
+ print value,plf.limits[Lower],plf.limits[Upper]
+ print weightup,weightlow,currentWeight[Upper],currentWeight[Lower]
return currentWeight
DonorScores = [elem for pos,elem in enumerate(don_supp) if SpliceAlign[pos] == 1]
assert not ( -inf in DonorScores )
- print 'donor'
+ #print 'donor'
weightDon = calculateWeights(d,DonorScores)
####################################################################################
AcceptorScores = [elem for pos,elem in enumerate(acc_supp) if pos > 0 and SpliceAlign[pos-1] == 2]
assert not ( -inf in AcceptorScores )
- print 'acceptor'
+ #print 'acceptor'
weightAcc = calculateWeights(a,AcceptorScores)
####################################################################################
def compute_donacc(donor_supp, acceptor_supp, d, a):
assert(len(donor_supp)==len(acceptor_supp))
+ size = len(donor_supp)
- #donor = zeros((len(donor_supp),1))
- #acceptor= zeros((len(acceptor_supp),1))
- donor = [0.0]*len(donor_supp)
- acceptor= [0.0]*len(acceptor_supp)
+ donor = [0.0]*size
+ acceptor = [0.0]*size
- for idx in range(len(donor_supp)):
+ for idx in range(size):
if isinf(donor_supp[idx]):
donor[idx] = donor_supp[idx]
else:
# Exons = TestExon ; % exon boundaries
# Ests = TestEsts ; % est sequences
+ # Lower all indices by one to convert matlab
+ # to python indices
+
+ Exons -= 1
+
return Sequences, Acceptors, Donors, Exons, Ests, Noises
# -*- coding: utf-8 -*-
import math
+import pdb
def penalty_lookup_new(penalty_struct, value):
# Berechne die Parameter des wirklichen Alignments (but with untrained d,a,h ...)
trueSpliceAlign, trueWeightMatch = computeSpliceAlign(dna, exons)
+
+ pdb.set_trace()
# Calculate
trueWeightDon, trueWeightAcc, trueWeightIntron = computeSpliceWeights(d, a, h, trueSpliceAlign, don_supp, acc_supp)
-
# Reshape currentW param
currentW = zeros((numFeatures,1))
currentW[numDonSuppPoints+numAccSuppPoints+numLengthSuppPoints+sizeMMatrix:numDonSuppPoints+numAccSuppPoints+numLengthSuppPoints+sizeMMatrix+numQualSuppPoints,0] = qualityMatrix[:,0]
currentPhi = zeros((numFeatures,1))
- currentPhi[0:30] = d.penalties[:]
- currentPhi[30:60] = a.penalties[:]
- currentPhi[60:90] = h.penalties[:]
+ currentPhi[0:30] = mat(d.penalties[:]).reshape(30,1)
+ currentPhi[30:60] = mat(a.penalties[:]).reshape(30,1)
+ currentPhi[60:90] = mat(h.penalties[:]).reshape(30,1)
currentPhi[90:126] = mmatrix[:]
currentPhi[126:] = qualityMatrix[:]
dna_len = len(dna)
est_len = len(est)
+ ps = h.convert2SWIG()
+
matchmatrix = QPalmaDP.createDoubleArrayFromList(mmatrix.flatten().tolist()[0])
mm_len = 36
- donor = QPalmaDP.createDoubleArrayFromList([1,2.0,4])
- d_len = 3
- acceptor = QPalmaDP.createDoubleArrayFromList([.1,2,4])
- a_len = 3
+
+ d_len = len(donor)
+ donor = QPalmaDP.createDoubleArrayFromList(donor)
+ a_len = len(acceptor)
+ acceptor = QPalmaDP.createDoubleArrayFromList(acceptor)
+
remove_duplicate_scores = False
print_matrix = False
currentAlignment = QPalmaDP.Alignment()
-
qualityMat = QPalmaDP.createDoubleArrayFromList(qualityMatrix)
currentAlignment.setQualityMatrix(qualityMat,numQualSuppPoints)
- ps = h.convert2SWIG()
+ print 'PYTHON: Calling myalign...'
# returns [SpliceAlign, EstAlign, weightMatch, Gesamtscores, dnaest]
currentAlignment.myalign( nr_paths, dna, dna_len,\
est, est_len, ps, matchmatrix, mm_len, donor, d_len,\
acceptor, a_len, remove_duplicate_scores, print_matrix)
- newSpliceAlign = QPalmaDP.createIntArrayFromList([0]*(dna_len*nr_paths))
- newEstAlign = QPalmaDP.createIntArrayFromList([0]*(est_len*nr_paths))
- newWeightMatch = QPalmaDP.createIntArrayFromList([0]*(mm_len*nr_paths))
- newAlignmentScores = QPalmaDP.createDoubleArrayFromList([.0]*nr_paths)
+ #newSpliceAlign = QPalmaDP.createIntArrayFromList([0]*(dna_len*nr_paths))
+ #newEstAlign = QPalmaDP.createIntArrayFromList([0]*(est_len*nr_paths))
+ #newWeightMatch = QPalmaDP.createIntArrayFromList([0]*(mm_len*nr_paths))
+ #newAlignmentScores = QPalmaDP.createDoubleArrayFromList([.0]*nr_paths)
- currentAlignment.getAlignmentResults(newSpliceAlign, newEstAlign, newWeightMatch, newAlignmentScores)
+ #currentAlignment.getAlignmentResults(newSpliceAlign, newEstAlign, newWeightMatch, newAlignmentScores)
- print newSpliceAlign
+ #print newSpliceAlign
print 'after myalign call...'
- if iteration_nr == 2:
+ if exampleId == 2:
break
iteration_nr += 1
ps = QPalmaDP.penalty_struct()
ps.limits = QPalmaDP.createDoubleArrayFromList(self.limits)
- ps.penalties = QPalmaDP.createDoubleArrayFromList([elem[0] for elem in self.penalties.tolist()])
+ ps.penalties = QPalmaDP.createDoubleArrayFromList(self.penalties)
ps.max_len = self.max_len
ps.min_len = self.min_len
# Gapfunktion
####################
h.limits = logspace(math.log(min_intron_len,10),math.log(max_intron_len,10),30)
- h.penalties = param[0:30]
+ h.penalties = param[0:30].flatten().tolist()[0]
#h.transform = '+1'
h.transform = ''
h.name = 'h'
# Donorfunktion
####################
d.limits = linspace(min_svm_score,max_svm_score,30)
- d.penalties = param[30:60]
+ d.penalties = param[30:60].flatten().tolist()[0]
#d.transform = '+1'
d.transform = ''
d.name = 'd'
# Acceptorfunktion
####################
a.limits = linspace(min_svm_score,max_svm_score,30)
- a.penalties = param[60:90]
+ a.penalties = param[60:90].flatten().tolist()[0]
#a.transform = '+1'
a.transform = ''
a.name = 'a'