+ fixed minor bugs
[qpalma.git] / qpalma / palma_utils.py
1 #
2 # This program is free software; you can redistribute it and/or modify
3 # it under the terms of the GNU General Public License as published by
4 # the Free Software Foundation; either version 2 of the License, or
5 # (at your option) any later version.
6 #
7 # Written (W) 2006-2007 Uta Schulze, Gunnar Raetsch
8 # Copyright (C) 2006-2007 Max-Planck-Society
9 #
10
11 def set_param(model):
12 import numpy
13
14 class plf: #means piecewise linear function
15 len = 0
16 limits = []
17 penalties = []
18 transform = ''
19 name = ''
20 max_len = 0
21 min_len = 0
22
23 h = plf()
24 h.len = int(model.intron_len_bins) ;
25 h.limits = model.intron_len_limits ;
26 h.penalties = model.intron_len_penalties ;
27 h.name = 'h'
28 h.max_len = int(model.intron_len_max)
29 h.min_len = int(model.intron_len_min)
30 h.transform = model.intron_len_transform
31
32 d = plf()
33 d.len = int(model.donor_bins)
34 d.limits = model.donor_limits
35 d.penalties = model.donor_penalties
36 d.name = 'd'
37 d.max_len = 100
38 d.min_len = -100
39
40 a = plf()
41 a.len = int(model.acceptor_bins)
42 a.limits = model.acceptor_limits
43 a.penalties = model.acceptor_penalties
44 a.name = 'a'
45 a.max_len = 100
46 a.min_len = -100
47
48 mmatrix = model.substitution_matrix
49 return h,d,a,mmatrix
50 #<------------------------------------------------------------------------------>
51
52
53 #<------------------------------------------------------------------------------>
54 # gives back y-value acc, to "support"-value (piecewise linear function)
55 #<------------------------------------------------------------------------------>
56 def penalty_lookup(plf, supp_value):
57 import numpy
58 limits = plf.limits
59 penalties = plf.penalties
60 transform = plf.transform
61
62 if transform == '':
63 supp_value = supp_value
64 elif transform == 'log':
65 print 'log'
66 supp_value = numpy.log(supp_value)
67 elif transform == 'log(+3)':
68 print 'log(+3)'
69 supp_value = numpy.log(supp_value+3)
70 elif transform == 'log(+1)':
71 print 'log(+1)'
72 supp_value = numpy.log(supp_value+1)
73 elif transform == '(+3)':
74 print '(+3)'
75 supp_value = supp_value+3
76 elif transform == 'mod3':
77 print 'mod3'
78 supp_value = numpy.mod(supp_value,3)
79 else:
80 raise NameError, 'unknown penalty transform'
81
82 #print 'supp_value = %i\n' %supp_value
83
84 if supp_value <= limits[0]:
85 score = penalties[0]
86 #print 'supp_value <= limits[0]'
87 elif supp_value >= limits[plf.len-1]:
88 score = penalties[plf.len-1]
89 #print 'supp_value >= limits[len(limits)-1]'
90 else:
91 for cnt in range(plf.len):
92 if supp_value <= limits[cnt]:
93 upper_idx = cnt
94 break
95 lower_idx = upper_idx-1
96 #print '%1.3f <= %1.3f <= %1.3f\n' %(limits[lower_idx], supp_value, limits[upper_idx])
97 score = ((supp_value-limits[lower_idx])*penalties[upper_idx] + (limits[upper_idx]-supp_value)*penalties[lower_idx])/(limits[upper_idx]-limits[lower_idx])
98
99 return score
100 #<------------------------------------------------------------------------------>
101
102
103 #<------------------------------------------------------------------------------>
104 # counts matches, mismatches, gaps on DNA (means bases inserted in query, qGapb)
105 # and gaps on EST (means baes inserted in transcript, tGapb)
106 #
107 # matchmatrix (columnwise)
108 # - A C G T N (dna)
109 # - x x x x x x
110 # A x x x x x x
111 # C x x x x x x
112 # G x x x x x x
113 # T x x x x x x
114 # N x x x x x x
115 # (est)
116 #<------------------------------------------------------------------------------>
117 def evaluate_mmatrix(mmatrix):
118 import numpy
119 #7: [a,a], 14:[c,c], 21:[g,g], 28:[t,t], 35:[n.n]
120 match = mmatrix[7] + mmatrix[14] + mmatrix[21] + mmatrix[28] + mmatrix[35]
121 qGapb = numpy.sum(mmatrix[1:6]) #gaps on DNA
122 tGapb = mmatrix[6] + mmatrix[12] + mmatrix[18] + mmatrix[24] + mmatrix[30] #gaps on EST
123 mism = numpy.sum(mmatrix) - (match + qGapb + tGapb)
124
125 if 0:
126 print 'matches: %d' %match
127 print 'mismatches: %d' %mism
128 print 'qGapb (gaps on DNA): %d' %qGapb
129 print 'tGapb (gaps on EST): %d' %tGapb
130
131 return match, mism, qGapb, tGapb
132 #<------------------------------------------------------------------------------>
133
134
135 #<------------------------------------------------------------------------------>
136 # qStart, qEnd, tStart, tEnd, blocks, qStarts, tStarts
137 #<------------------------------------------------------------------------------>
138 def get_splice_info(splice_align, est_align):
139 import sys
140 #problem: because of the dangling ends, alignment can start with an intron
141
142 #print splice_align
143 #print est_align
144
145 dna_len = len(splice_align)
146 est_len = len(est_align)
147 qStart = 1 #(default)
148 qEnd = est_len #(default)
149 tStart = 1 #(default)
150 tEnd = dna_len #(default)
151 qStarts = []
152 qEnds = []
153 qExonSizes = []
154 tStarts = []
155 tEnds = []
156 tExonSizes = []
157
158 #<---------------------------------------------------------->
159 #query: dangling start
160 if est_align[0]==4:
161 for i in range(est_len):
162 if est_align[i]!=4:
163 qStart = i+1 #+1 bc index starts with zero
164 break
165
166 #query: dangling end
167 if est_align[est_len-1]==4:
168 list_idx = range(est_len)
169 list_idx.reverse()
170 for i in list_idx:
171 if est_align[i]!=4:
172 qEnd = i+1 #last "4", +1 bc index starts with zero
173 break
174 #qStarts: list of exon starting points on EST
175 #qEnds: list of exon ending points on EST
176 #(2 in est_align means: first nucleotideof new exon)
177
178 index_count = 0
179 exon_label = est_align[qStart-1]
180 qStarts.append(qStart)
181 for i in range(qStart-1,qEnd):
182 if (est_align[i]==1) and (exon_label==2): #new exon labeled "1"
183 exon_label = 1
184 qStarts.append(i+1) #+1 bc index starts with zero, start of this exon saved
185 qEnds.append(i-1+1) #end of last exon labeled "2" saved
186 qExonSizes.append(qEnds[index_count] - qStarts[index_count] +1) #length of last exon
187 index_count+=1
188 elif (est_align[i]==2) and (exon_label==1): #new exon labeled "2"
189 exon_label = 2
190 qStarts.append(i+1) #+1 bc index starts with zero, start of this exon saved
191 qEnds.append(i-1+1) #end of last exon labeled "2" saved
192 qExonSizes.append(qEnds[index_count] - qStarts[index_count] +1) #length of last exon
193 index_count+=1
194 if (i == qEnd-1): #end of last exon
195 qEnds.append(i+1) #end of last exon saved
196 qExonSizes.append(qEnds[index_count] - qStarts[index_count] +1) #length of last exon
197
198 assert(len(qStarts)==len(qEnds))
199
200 if 0:
201 print 'qStart: %d' %qStart
202 print 'qEnd: %d' %qEnd
203 print qStarts
204 print qEnds
205 print qExonSizes
206 #<---------------------------------------------------------->
207
208 #<---------------------------------------------------------->
209 #transcript: dangling start
210 if splice_align[0] == 4:
211 for i in range(dna_len):
212 if splice_align[i]!=4:
213 tStart = i+1 #+1 bc index starts with zero
214 break
215 #if splice_align[tStart-1]==1: #alignment starts with intron
216 #print 'alignment starts with intron'
217
218 #transcript: dangling end
219 if splice_align[dna_len-1]==4:
220 list_idx = range(dna_len)
221 list_idx.reverse()
222 for i in list_idx:
223 if splice_align[i]!=4:
224 tEnd = i+1 #+1 bc index starts with zero
225 break
226 #if splice_align[tEnd-1]==2: #alignment ends with intron
227 #print 'alignment ends with intron'
228
229 #tStarts: list of exon starting points on DNA
230 #tEnds: list of exon ending points on DNA
231
232 index_count = 0
233 act_state = 3 #3 means intron, 0 means exon
234 for i in range(tStart-1,tEnd):
235 if (splice_align[i]==0) and (act_state==3): #first base of exon
236 tStarts.append(i+1) #+1 bc index starts with zero
237 act_state=0 #exon
238 elif (splice_align[i]==1) and (act_state==0): #exon ended one base ago
239 tEnds.append(i-1+1) #+1 bc index starts with zero
240 act_state=3 #intron
241 tExonSizes.append(tEnds[index_count] - tStarts[index_count] +1) #length
242 index_count+=1
243 if (i==tEnd-1) and (act_state==0): #end of alignment reached and still in exon
244 tEnds.append(i+1) #+1 bc index starts with zero
245 tExonSizes.append(tEnds[index_count] - tStarts[index_count] +1) #length
246
247 #<---------------------------------------------------------->
248
249
250
251 if (len(tStarts)!=len(tEnds)):
252 print len(tStarts)
253 print len(tEnds)
254 print 'tStart: %d' %tStart
255 print 'tEnd: %d' %tEnd
256 print tStarts
257 print tEnds
258 sys.exit(2)
259
260 #print 'exons on query: %d, exons on target: %d' %(len(qStarts),len(tStarts))
261 num_exons = len(tStarts)
262 #print num_exons
263 #print len(qStarts)
264
265 if 0:
266 print 'tStart: %d' %tStart
267 print 'tEnd: %d' %tEnd
268 print tStarts
269 print tEnds
270 print tExonSizes
271
272 # assert(num_exons==len(qStarts))
273
274 return qStart, qEnd, tStart, tEnd, num_exons, qExonSizes, qStarts, qEnds, tExonSizes, tStarts, tEnds
275 #<------------------------------------------------------------------------------>
276
277 """ read fasta as dictionary """
278 def read_fasta(f):
279 import sys
280 fasta=dict()
281 for s in f.readlines():
282 if s.startswith('>'):
283 key=s[1:-1]
284 print key
285 if fasta.has_key(key):
286 sys.stderr.write("ERROR: duplicated sequence name '%s'\n" %(key))
287 sys.exit(2)
288 fasta[key]=""
289 else:
290 fasta[key]+=s[:-1]
291
292 return fasta
293
294 #<------------------------------------------------------------------------------>
295 #
296 # DNA and EST with gaps and est with intron
297 # len(dna) == len(est)
298 # identity: list of length exon_numbers: #matches / length(exon_with_gaps)
299 #<------------------------------------------------------------------------------>
300 def comp_identity(DNA, EST, strand, qStart, tStart, translation, comparison_char):
301 def set_idx(i,qNum,tNum):
302 qIDX[i] = qNum
303 tIDX[i] = tNum
304 return qIDX, tIDX
305
306 gap_char = translation[0] #e.g. "-"
307 in_char = translation[6] #e.g. "*", intron char
308 equ_char = comparison_char[0] #match char e.g. "|"
309 comp_in_char = comparison_char[1] #intron char when comparing
310 comparison = [' '] *len(EST)
311 qIDX = [-1] *len(EST)
312 tIDX = [-1] *len(DNA)
313 exon_idx = 0
314 matches = []
315 gaps = []
316 exon_length = []
317 identity = []
318 exonQuality = [] # matches - 2*gaps
319 ssQuality = []
320 state = 3 #intron:3, exon:0
321 EST_started = 0
322
323 for i in range(len(DNA)):
324 if (i>0 and i<len(DNA)-1):
325 if (EST[i] == in_char): #intron
326 (qIDX, tIDX) = set_idx(i,qIDX[i-1],tIDX[i-1]+1)
327 comparison[i] = comp_in_char
328 if (state==0): #exon finished
329 identity.append(float(matches[exon_idx])/float(exon_length[exon_idx])) #save identity for finished exon
330 exonQuality.append(float(matches[exon_idx]-2*gaps[exon_idx])/float(exon_length[exon_idx])) #save identity for finished exon
331 exon_idx+=1 #next exon starts sooner or later
332 state = 3 #now intron
333 last_ssQuality = 0.0 ;
334 last_ssLen = 0 ;
335 for k in range(5):
336 if i-k-1<0: break
337 last_ssLen+=1
338 if DNA[i-k-1]==EST[i-k-1]: last_ssQuality+=1.0
339 else:
340 if (DNA[i-k-1]=='N' and EST[i-k-1]!=gap_char) or (DNA[i-k-1]!=gap_char and EST[i-k-1]=='N'): last_ssQuality+=0.5
341 #print k,DNA[i-k-1],EST[i-k-1],last_ssQuality
342 else: #exon
343 if (state==0): #still in same exon
344 if EST_started and DNA[i]!=gap_char: # only count DNA length
345 exon_length[exon_idx] += 1 #one bp added
346 if EST[i]!=gap_char: EST_started=1 ;
347 #print i,EST[i],EST_started,exon_length[exon_idx]
348 if (DNA[i]==EST[i]):
349 matches[exon_idx]+=1 #one match added
350 else:
351 if (DNA[i]=='N' and EST[i]!=gap_char) or (EST[i]=='N' and DNA[i]!=gap_char): matches[exon_idx]+=0.5
352 else: #new exon
353 for k in range(5):
354 if i+k>=len(DNA) or i+k>len(EST): break
355 last_ssLen+=1 ;
356 if DNA[i+k]==EST[i+k]: last_ssQuality+=1.0 ;
357 else:
358 if (DNA[i+k]=='N' and EST[i+k]!=gap_char) or (DNA[i+k]!=gap_char and EST[i+k]=='N'): last_ssQuality+=0.5
359 #print k,DNA[i+k],EST[i+k],last_ssQuality
360 ssQuality.append(last_ssQuality/last_ssLen)
361 if DNA[i]!=gap_char: exon_length.append(1)
362 else: exon_length.append(0)
363 state = 0 #now exon
364 if (DNA[i]==EST[i]):
365 gaps.append(0)
366 matches.append(1)
367 else:
368 if (DNA[i]=='N' and EST[i]!=gap_char) or (EST[i]=='N' and DNA[i]!=gap_char): matches.append(0.5); gaps.append(0)
369 else: matches.append(0); gaps.append(1)
370
371 if (DNA[i]==EST[i]) and (DNA[i]!='-'): #match
372 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1]+1)
373 comparison[i] = equ_char
374 elif (DNA[i]==EST[i]) and (DNA[i]=='-'): #strange match
375 (qIDX, tIDX) = set_idx(i,qIDX[i-1],tIDX[i-1])
376 elif (EST[i]==gap_char): #gap on est
377 (qIDX, tIDX) = set_idx(i,qIDX[i-1],tIDX[i-1]+1)
378 elif (DNA[i]==gap_char): #gap on dna
379 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1])
380 else: #mismatch
381 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1]+1)
382
383 elif (i==0):
384 if (EST[i] == in_char): #alignment starts with intron
385 (qIDX, tIDX) = set_idx(i,0,1)
386 comparison[i] = comp_in_char
387 else: #alignment starts with exon
388 state = 0 #go into exon_state
389 if DNA[i]!=gap_char: exon_length.append(1)
390 else: exon_length.append(0)
391 if (DNA[i]==EST[i]) and (DNA[i]!='-'): #match
392 #(qIDX, tIDX) = set_idx(i,1,1)
393 EST_started = 1
394 (qIDX, tIDX) = set_idx(i,qStart,tStart)
395 matches.append(1)
396 gaps.append(0)
397 comparison[i] = equ_char
398 elif (DNA[i]==EST[i]) and (DNA[i]=='-'): #strange match
399 (qIDX, tIDX) = set_idx(i,qStart,tStart)
400 matches.append(0)
401 gaps.append(1)
402 elif (EST[i]==gap_char): #gap on est
403 #(qIDX, tIDX) = set_idx(i,0,1)
404 (qIDX, tIDX) = set_idx(i,qStart,tStart)
405 matches.append(0)
406 gaps.append(1)
407 elif (DNA[i]==gap_char): #gap on dna
408 #(qIDX, tIDX) = set_idx(i,1,0)
409 (qIDX, tIDX) = set_idx(i,qStart,tStart)
410 matches.append(0)
411 gaps.append(1)
412 else: #mismatch
413 #(qIDX, tIDX) = set_idx(i,1,1)
414 (qIDX, tIDX) = set_idx(i,qStart,tStart)
415 gaps.append(0)
416 if (DNA[i]=='N' or EST[i]=='N'): matches.append(0.5)
417 else: matches.append(0.5)
418 elif (i==len(DNA)-1):
419 if (EST[i] == in_char): #alignment ends with intron, with exons everything ok
420 assert(state==3) #ok, because intron longer than 4
421 comparison[i] = comp_in_char
422 (qIDX, tIDX) = set_idx(i,qIDX[i-1],tIDX[i-1]+1)
423 else: #alignment ends with exon, exon finished
424 if (state==0): #finish exon
425 if DNA[i]!=gap_char: exon_length[exon_idx] += 1 #one bp added
426 else: exon_length[exon_idx] += 1 #one bp added
427 if (DNA[i]==EST[i]) and (DNA[i]!=gap_char): #match
428 matches[exon_idx]+=1 #one match added
429 comparison[i] = equ_char
430 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1]+1)
431 elif (DNA[i]==EST[i]) and (DNA[i]==gap_char): #strange match
432 (qIDX, tIDX) = set_idx(i,qIDX[i-1],tIDX[i-1])
433 elif (EST[i]==gap_char): #gap on est
434 (qIDX, tIDX) = set_idx(i,qIDX[i-1],tIDX[i-1]+1)
435 elif (DNA[i]==gap_char): #gap on dna
436 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1])
437 else: #mismatch
438 if (DNA[i]=='N' or EST[i]=='N'): matches[exon_idx]+=0.5
439 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1]+1)
440 identity.append(float(matches[exon_idx])/float(exon_length[exon_idx])) #save identity for finished exon
441 exonQuality.append(float(matches[exon_idx]-2*gaps[exon_idx])/float(exon_length[exon_idx])) #save identity for finished exon
442 else: #last exon has length 1
443 if DNA[i]!=gap_char: exon_length.append(1)
444 else: exon_length.append(1)
445 if (DNA[i]==EST[i]): #match
446 matches.append(1) #one match added
447 gaps.append(0)
448 comparison[i] = equ_char
449 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1]+1)
450 elif (EST[i]==gap_char): #gap on est
451 gaps.append(1)
452 matches.append(0)
453 (qIDX, tIDX) = set_idx(i,qIDX[i-1],tIDX[i-1]+1)
454 elif (DNA[i]==gap_char): #gap on dna
455 gaps.append(1)
456 matches.append(0)
457 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1])
458 else: #mismatch
459 gaps.append(0)
460 if (DNA[i]=='N' or EST[i]=='N'): matches.append(0.5)
461 else: matches.append(0)
462 (qIDX, tIDX) = set_idx(i,qIDX[i-1]+1,tIDX[i-1]+1)
463 identity.append(float(matches[exon_idx])/float(exon_length[exon_idx])) #save identity for finished exon
464 exonQuality.append(float(matches[exon_idx]-2*gaps[exon_idx])/float(exon_length[exon_idx])) #save identity for finished exon
465 else:
466 print 'this "else" should not happen'
467
468 comparison = reduce(lambda x,y: x+y, comparison)
469 #print qIDX
470 #print tIDX
471 #print exon_length
472 #print ssQuality, exonQuality, identity
473
474 return exon_length, identity, ssQuality, exonQuality, comparison, qIDX, tIDX
475 #<------------------------------------------------------------------------------>
476
477
478 #<------------------------------------------------------------------------------>
479 #
480 #<------------------------------------------------------------------------------>
481 def reverse_complement(DNA):
482 import string
483 import sys
484 dict = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'N': 'N'}
485
486 DNA = DNA.upper()
487 reverse_DNA = ''
488 cnt_list = range(len(DNA))
489 cnt_list.reverse()
490 for i in cnt_list:
491 if not dict.has_key(DNA[i]):
492 print 'unknown letter', DNA[i], '-> translating to N'
493 reverse_DNA = reverse_DNA + 'N' ;
494 else:
495 reverse_DNA = reverse_DNA + dict[DNA[i]]
496 return reverse_DNA
497 #<------------------------------------------------------------------------------>
498
499
500 #<------------------------------------------------------------------------------>
501 #
502 #<------------------------------------------------------------------------------>
503 def parse_wublast(output, e_value_bound):
504 import numpy
505 import sys
506
507 strand_matrix='- +'
508
509 class wublast_info:
510 score = []
511 q_name = []
512 t_name = []
513 q_strand = []
514 t_strand = []
515 q_start = []
516 q_end = []
517 t_start = []
518 t_end = []
519
520 est_entry = 0
521 still_in_group = 0
522
523 est_wbi = wublast_info()
524 for line in output:
525 fields = line.strip().split('\t')
526 #print line, still_in_group
527 if float(fields[2]) <= e_value_bound:
528 if len(est_wbi.q_name)>0:
529 flag = (fields[0] == est_wbi.q_name[len(est_wbi.t_name)-1]) and (fields[1] == est_wbi.t_name[len(est_wbi.t_name)-1]) and (int(fields[3]) == group_size)
530 else:
531 flag = True
532
533 if (still_in_group == 0) or (not flag):
534 if (not still_in_group == 0):
535 sys.stderr.write('blastn output format problem near line:\n'+line)
536
537 # start a new group
538 est_wbi.score.append(float(fields[4]))
539 est_wbi.q_name.append(fields[0])
540 est_wbi.t_name.append(fields[1])
541 est_wbi.q_strand.append(strand_matrix[int(fields[16])+1])
542 est_wbi.t_strand.append(strand_matrix[int(fields[19])+1])
543 if strand_matrix[int(fields[16])+1] == '+':
544 est_wbi.q_start.append(int(fields[17]))
545 est_wbi.q_end.append(int(fields[18]))
546 else:
547 est_wbi.q_start.append(int(fields[18]))
548 est_wbi.q_end.append(int(fields[17]))
549 est_wbi.t_start.append(int(fields[20]))
550 est_wbi.t_end.append(int(fields[21]))
551 still_in_group = int(fields[3])-1
552 #print 'still_in_group', still_in_group
553 group_size = int(fields[3])
554 else:
555 # extend a group
556 assert(fields[0] == est_wbi.q_name[len(est_wbi.t_name)-1])
557 assert(fields[1] == est_wbi.t_name[len(est_wbi.t_name)-1])
558 assert(int(fields[3]) == group_size)
559 est_wbi.score[len(est_wbi.score)-1] += float(fields[4])
560 if strand_matrix[int(fields[16])+1] == '+':
561 est_wbi.q_start[len(est_wbi.q_start)-1] = numpy.minimum(est_wbi.q_start[len(est_wbi.q_start)-1],int(fields[17]))
562 est_wbi.q_end[len(est_wbi.q_end)-1] =numpy.maximum(est_wbi.q_end[len(est_wbi.q_end)-1],int(fields[18]))
563 else:
564 est_wbi.q_start[len(est_wbi.q_start)-1] =numpy.minimum(est_wbi.q_start[len(est_wbi.q_start)-1],int(fields[18]))
565 est_wbi.q_end[len(est_wbi.q_end)-1] =numpy.maximum(est_wbi.q_end[len(est_wbi.q_end)-1],int(fields[17]))
566 est_wbi.t_start[len(est_wbi.t_start)-1] =numpy.minimum(est_wbi.t_start[len(est_wbi.t_start)-1],int(fields[20]))
567 est_wbi.t_end[len(est_wbi.t_end)-1] =numpy.maximum(est_wbi.t_end[len(est_wbi.t_end)-1],int(fields[21]))
568 still_in_group -= 1
569
570 if 0:
571 print est_wbi.q_name
572 print est_wbi.t_name
573 print est_wbi.q_strand
574 print est_wbi.t_strand
575 print est_wbi.q_start
576 print est_wbi.q_end
577 print est_wbi.t_start
578 print est_wbi.t_end
579
580 return est_wbi
581 #<------------------------------------------------------------------------------>