+ minor changes in qpalma_main to take into accout that jobs may be rescheduled
[qpalma.git] / scripts / Evaluation.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 import cPickle
5 import sys
6 import pydb
7 import pdb
8 import os
9 import os.path
10 import math
11
12
13 def createTable(results):
14 """
15 This function gets the results from the evaluation and creates a tex table.
16 """
17
18 fh = open('result_table.tex','w+')
19 lines = ['\begin{tabular}{|c|c|c|r|}', '\hline',\
20 'Quality & Splice & Intron & \multicolumn{1}{c|}{Error on Positions} & \multicolumn{1}{c|}{Error on Scores} & \\',\
21 'information & site pred. & length & \multicolumn{1}{c|}{rate}\\',\
22 '\hline',\
23
24 #'- & - & - & %2.3f \%\\' % ( results['---'][0], results['---'][1] ),\
25 #'+ & - & - & %2.3f \%\\' % ( results['+--'][0], results['+--'][1] ),\
26
27 #'\hline',\
28
29 #'- & + & - & %2.3f \%\\' % ( results['-+-'][0], results['-+-'][1] ),\
30 #'+ & + & - & %2.3f \%\\' % ( results['++-'][0], results['++-'][1] ),\
31
32 '\hline\hline',\
33
34 '- & - & + & %2.3f & %2.3f \\%%\\\\' % ( results['--+'][0], results['--+'][1] ),\
35 '+ & - & + & %2.3f & %2.3f \\%%\\\\' % ( results['+-+'][0], results['+-+'][1] ),\
36
37 '\hline',\
38
39 '- & + & + & %2.3f & %2.3f \\%%\\\\' % ( results['-++'][0], results['-++'][1] ),\
40 '+ & + & + & %2.3f & %2.3f \\%%\\\\' % ( results['+++'][0], results['+++'][1] ),\
41
42 '\hline',\
43 '\end{tabular}']
44
45
46 lines = [l+'\n' for l in lines]
47 for l in lines:
48 fh.write(l)
49 fh.close()
50
51
52 def compare_scores_and_labels(scores,labels):
53 """
54 Iterate through all predictions. If we find a correct prediction check
55 whether this correct prediction scores higher than the incorrect
56 predictions for this example.
57 """
58
59 for currentPos,currentElem in enumerate(scores):
60 if labels[currentPos] == True:
61 for otherPos,otherElem in enumerate(scores):
62 if otherPos == currentPos:
63 continue
64
65 if labels[otherPos] == False and otherElem > currentElem:
66 return False
67
68 return True
69
70
71 def evaluate_unmapped_example(current_prediction):
72 predExons = current_prediction['predExons']
73 trueExons = current_prediction['trueExons']
74
75 result = compare_exons(predExons,trueExons)
76 return result
77
78
79 def evaluate_example(current_prediction):
80 label = False
81 label = current_prediction['label']
82
83 pred_score = current_prediction['DPScores'].flatten().tolist()[0][0]
84
85 # if the read was mapped by vmatch at an incorrect position we only have to
86 # compare the score
87 if label == False:
88 return label,False,pred_score
89
90 predExons = current_prediction['predExons']
91 trueExons = current_prediction['trueExons']
92
93 predPositions = [elem + current_prediction['alternative_start_pos'] for elem in predExons]
94 truePositions = [elem + current_prediction['start_pos'] for elem in trueExons.flatten().tolist()[0]]
95
96 #pdb.set_trace()
97
98 pos_comparison = (predPositions == truePositions)
99
100 #if label == True and pos_comparison == False:
101 # pdb.set_trace()
102
103 return label,pos_comparison,pred_score
104
105
106 def compare_exons(predExons,trueExons):
107 e1_b_off,e1_e_off,e2_b_off,e2_e_off = 0,0,0,0
108
109 if len(predExons) == 4:
110 e1_begin,e1_end = predExons[0],predExons[1]
111 e2_begin,e2_end = predExons[2],predExons[3]
112 else:
113 return False
114
115 e1_b_off = int(math.fabs(e1_begin - trueExons[0,0]))
116 e1_e_off = int(math.fabs(e1_end - trueExons[0,1]))
117
118 e2_b_off = int(math.fabs(e2_begin - trueExons[1,0]))
119 e2_e_off = int(math.fabs(e2_end - trueExons[1,1]))
120
121 if e1_b_off == 0 and e1_e_off == 0 and e2_b_off == 0\
122 and e2_e_off == 0:
123 return True
124
125 return False
126
127 def prediction_on(filename):
128 allPredictions = cPickle.load(open(filename))
129
130 exon1Begin = []
131 exon1End = []
132 exon2Begin = []
133 exon2End = []
134 allWrongExons = []
135 allDoubleScores = []
136
137 gt_correct_ctr = 0
138 pos_correct_ctr = 0
139 pos_incorrect_ctr = 0
140 score_correct_ctr = 0
141 score_incorrect_ctr = 0
142
143 total_vmatch_instances_ctr = 0
144 true_vmatch_instances_ctr = 0
145
146 for current_example_pred in allPredictions:
147 gt_example = current_example_pred[0]
148 gt_score = gt_example['DPScores'].flatten().tolist()[0][0]
149 gt_correct = evaluate_unmapped_example(gt_example)
150
151 current_scores = []
152 current_labels = []
153 current_scores.append(gt_score)
154 current_labels.append(gt_correct)
155
156 if gt_correct:
157 gt_correct_ctr += 1
158
159 for elem_nr,current_pred in enumerate(current_example_pred[1:]):
160 current_label,comparison_result,current_score = evaluate_example(current_pred)
161
162 # if vmatch found the right read pos we check for right exons
163 # boundaries
164 if current_label:
165 if comparison_result:
166 pos_correct_ctr += 1
167 else:
168 pos_incorrect_ctr += 1
169
170 true_vmatch_instances_ctr += 1
171
172 current_scores.append(current_score)
173 current_labels.append(current_label)
174
175 total_vmatch_instances_ctr += 1
176
177 # check whether the correct predictions score higher than the incorrect
178 # ones
179 cmp_res = compare_scores_and_labels(current_scores,current_labels)
180 if cmp_res:
181 score_correct_ctr += 1
182 else:
183 score_incorrect_ctr += 1
184
185 # now that we have evaluated all instances put out all counters and sizes
186 print 'Total num. of examples: %d' % len(allPredictions)
187 print 'Number of correct ground truth examples: %d' % gt_correct_ctr
188 print 'Total num. of true vmatch instances %d' % true_vmatch_instances_ctr
189 print 'Correct pos: %d, incorrect pos: %d' %\
190 (pos_correct_ctr,pos_incorrect_ctr)
191 print 'Total num. of vmatch instances %d' % total_vmatch_instances_ctr
192 print 'Correct scores: %d, incorrect scores: %d' %\
193 (score_correct_ctr,score_incorrect_ctr)
194
195 pos_error = 1.0 * pos_incorrect_ctr / true_vmatch_instances_ctr
196 score_error = 1.0 * score_incorrect_ctr / total_vmatch_instances_ctr
197
198 print pos_error,score_error
199
200 return (pos_error,score_error)
201
202 def collect_prediction(current_dir,run_name):
203 """
204 Given the toplevel directory this function takes care that for each distinct
205 experiment the training and test predictions are evaluated.
206
207 """
208 train_suffix = '_allPredictions_TRAIN'
209 test_suffix = '_allPredictions_TEST'
210
211 jp = os.path.join
212 b2s = ['-','+']
213
214 currentRun = cPickle.load(open(jp(current_dir,'run_object.pickle')))
215 QFlag = currentRun['enable_quality_scores']
216 SSFlag = currentRun['enable_splice_signals']
217 ILFlag = currentRun['enable_intron_length']
218 currentRunId = '%s%s%s' % (b2s[QFlag],b2s[SSFlag],b2s[ILFlag])
219
220 filename = jp(current_dir,run_name)+train_suffix
221 print 'Prediction on: %s' % filename
222 train_result = prediction_on(filename)
223
224 filename = jp(current_dir,run_name)+test_suffix
225 print 'Prediction on: %s' % filename
226 test_result = prediction_on(filename)
227
228 return train_result,test_result,currentRunId
229
230
231 def perform_prediction(current_dir,run_name):
232 """
233 This function takes care of starting the jobs needed for the prediction phase
234 of qpalma
235 """
236 cmd = 'echo /fml/ag-raetsch/home/fabio/svn/projects/QPalma/scripts/doPrediction.sh %s | qsub -l h_vmem=12.0G -cwd -j y -N \"%s.log\"'%(current_dir,run_name)
237 #cmd = './doPrediction.sh %s 1>%s.out 2>%s.err' %(current_dir,run_name,run_name)
238 #print cmd
239 os.system(cmd)
240
241
242 def forall_experiments(current_func,tl_dir):
243 """
244 Given the toplevel directoy this function calls for each subdir the
245 function given as first argument. Which are at the moment:
246
247 - perform_prediction, and
248 - collect_prediction.
249
250 """
251
252 dir_entries = os.listdir(tl_dir)
253 dir_entries = [os.path.join(tl_dir,de) for de in dir_entries]
254 run_dirs = [de for de in dir_entries if os.path.isdir(de)]
255
256 all_results = {}
257
258 for current_dir in run_dirs:
259 run_name = current_dir.split('/')[-1]
260 current_func(current_dir,run_name)
261
262 #train_result,test_result,currentRunId = current_func(current_dir,run_name)
263 #all_results[currentRunId] = test_result
264
265 #createTable(all_results)
266
267
268 if __name__ == '__main__':
269 dir = sys.argv[1]
270 assert os.path.exists(dir), 'Error directory does not exist!'
271
272 #forall_experiments(perform_prediction,dir)
273 forall_experiments(collect_prediction,dir)
274
275 """
276 def evaluateExample(dna,est,exons,SpliceAlign,newEstAlign):
277 newExons = []
278 oldElem = -1
279 SpliceAlign = SpliceAlign.flatten().tolist()[0]
280 SpliceAlign.append(-1)
281 for pos,elem in enumerate(SpliceAlign):
282 if pos == 0:
283 oldElem = -1
284 else:
285 oldElem = SpliceAlign[pos-1]
286
287 if oldElem != 0 and elem == 0: # start of exon
288 newExons.append(pos)
289
290 if oldElem == 0 and elem != 0: # end of exon
291 newExons.append(pos)
292
293 e1_b_off,e1_e_off,e2_b_off,e2_e_off = 0,0,0,0
294
295 if len(newExons) == 4:
296 e1_begin,e1_end = newExons[0],newExons[1]
297 e2_begin,e2_end = newExons[2],newExons[3]
298 else:
299 return None,None,None,None,newExons
300
301 e1_b_off = int(math.fabs(e1_begin - exons[0,0]))
302 e1_e_off = int(math.fabs(e1_end - exons[0,1]))
303
304 e2_b_off = int(math.fabs(e2_begin - exons[1,0]))
305 e2_e_off = int(math.fabs(e2_end - exons[1,1]))
306
307 pdb.set_trace()
308
309 return e1_b_off,e1_e_off,e2_b_off,e2_e_off,newExons
310 """
311
312 """
313 def evaluatePositions(eBegin,eEnd):
314 eBegin_pos = [elem for elem in eBegin if elem == 0]
315 eBegin_neg = [elem for elem in eBegin if elem != 0]
316 eEnd_pos = [elem for elem in eEnd if elem == 0]
317 eEnd_neg = [elem for elem in eEnd if elem != 0]
318
319 mean_eBegin_neg = 0
320 for idx in range(len(eBegin_neg)):
321 mean_eBegin_neg += eBegin_neg[idx]
322
323 try:
324 mean_eBegin_neg /= 1.0*len(eBegin_neg)
325 except:
326 mean_eBegin_neg = -1
327
328 mean_eEnd_neg = 0
329 for idx in range(len(eEnd_neg)):
330 mean_eEnd_neg += eEnd_neg[idx]
331
332 try:
333 mean_eEnd_neg /= 1.0*len(eEnd_neg)
334 except:
335 mean_eEnd_neg = -1
336
337 return eBegin_pos,eBegin_neg,eEnd_pos,eEnd_neg,mean_eBegin_neg,mean_eEnd_neg
338 """