936a163cf35662dcfdb556c36c38c0e9adb8ef5b
[qpalma.git] / tests / test_qpalma_prediction.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 import cPickle
5 import math
6 import numpy
7 import os.path
8 import pdb
9 import unittest
10
11 from qpalma_main import QPalma
12 from Utils import print_prediction
13
14 from createAlignmentFileFromPrediction import alignment_reconstruct
15
16 from qpalma.sequence_utils import DataAccessWrapper,SeqSpliceInfo
17
18 jp = os.path.join
19
20
21 class TestQPalmaPrediction(unittest.TestCase):
22 """
23 This class...
24 """
25
26 def _setUp(self):
27 data_fn = '/fml/ag-raetsch/home/fabio/tmp/vmatch_evaluation/spliced_1/dataset/prediction_set.pickle'
28
29 self.prediction_set = cPickle.load(open(data_fn))
30
31
32
33 def setUp(self):
34 print
35 self.prediction_set = {}
36
37 # xxxxxxxxxxx_______________________________________________________________________xxxxxxxxxxxxxxxxxxxxxxxxxxx
38 # 'tattttggaaggtatttcatctctccgacattgctctcaacactgtcccctccaatgcctagtccttttatttttttcttagttccaattcccttaaatacatctcacagtcttcttcttcttctcgattgcagtagc'
39
40 read = 'tattttggaagttccaattcccttaatacatctcacag'
41 currentQualities = [[30]*len(read)]
42
43 id = 1
44 chromo = 3
45 strand = '-'
46
47 tsize = 23470805
48
49 genomicSeq_start = tsize - (2038+10)
50 genomicSeq_stop = tsize - 1900 + 10
51 print 'Position: ',
52 print genomicSeq_start,genomicSeq_stop
53
54 currentSeqInfo = id,chromo,strand,genomicSeq_start,genomicSeq_stop
55
56 example = (currentSeqInfo,read,currentQualities)
57 self.prediction_set[id] = [example]
58
59 # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxx
60 # 'gctactgcaatcgagaagaagaagaagactgtgagatgtatttaagggaattggaactaagaaaaaaataaaaggactaggcattggaggggacagtgttgagagcaatgtcggagagatgaaataccttccaaaata'
61
62 read = 'gctactgcaatcgagaagaagaagaagactatgaaata'
63 currentQualities = [[40]*30+[22]*8]
64
65 id = 2
66 chromo = 3
67 strand = '+'
68
69 tsize = 23470805
70
71 genomicSeq_start = tsize - (2038+10)
72 genomicSeq_stop = tsize - 1900 + 10
73 print 'Position: ',
74 print genomicSeq_start,genomicSeq_stop
75
76 currentSeqInfo = id,chromo,strand,genomicSeq_start,genomicSeq_stop
77
78 example = (currentSeqInfo,read,currentQualities)
79 self.prediction_set[id] = [example]
80
81
82 def testAlignments(self):
83 run_dir = '/fml/ag-raetsch/home/fabio/tmp/newest_run/alignment/saved_run'
84
85 run = cPickle.load(open(jp(run_dir,'run_obj.pickle')))
86 run['name'] = 'test_run'
87 run['result_dir'] = '/fml/ag-raetsch/home/fabio/tmp/sandbox/testcases'
88
89 param_fn = jp(run_dir,'param_526.pickle')
90 param = cPickle.load(open(param_fn))
91
92 print self.prediction_set
93 for example_key in self.prediction_set.keys():
94 print 'Current example %d' % example_key
95
96 for example in self.prediction_set[example_key]:
97 print example
98 print 'size'
99 print len(example)
100
101
102 # fetch the data needed
103 g_dir = run['dna_flat_files'] #'/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
104 acc_dir = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/acc'
105 don_dir = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/don'
106
107 g_fmt = 'chr%d.dna.flat'
108 s_fmt = 'contig_%d%s'
109
110 num_chromo = 6
111
112 accessWrapper = DataAccessWrapper(g_dir,acc_dir,don_dir,g_fmt,s_fmt)
113 seqInfo = SeqSpliceInfo(accessWrapper,range(1,num_chromo))
114
115
116
117 qp = QPalma(True)
118 #qp.init_prediction(run,set_name)
119 allPredictions = qp.predict(run,self.prediction_set,param,seqInfo )
120 for current_prediction in allPredictions:
121 align_str = print_prediction(current_prediction)
122 print align_str
123
124 id = current_prediction['id']
125 seq = current_prediction['read']
126 dna = current_prediction['dna']
127 chromo = current_prediction['chr']
128 strand = current_prediction['strand']
129 start_pos = current_prediction['start_pos']
130 predExons = current_prediction['predExons']
131
132 numExons = int(math.ceil(len(predExons) / 2))
133
134 print alignment_reconstruct(current_prediction,numExons)
135 #print id,start_pos,predExons
136
137
138 if __name__ == '__main__':
139 suite = unittest.TestLoader().loadTestsFromTestCase(TestQPalmaPrediction)
140 unittest.TextTestRunner(verbosity=2).run(suite)