|
20
|
1 '''
|
|
|
2 usage:
|
|
|
3
|
|
|
4 python endbias.py utr5-coverage utr3-coverage outputfile
|
|
|
5 '''
|
|
|
6 import sys,math
|
|
|
7
|
|
|
8 def getCoverage(filename):
|
|
|
9 f = open(filename)
|
|
|
10 coverage = {}
|
|
|
11 for line in f:
|
|
|
12 flds = line.strip().split('\t')
|
|
|
13 score = float(flds[4])
|
|
|
14 name = (flds[0].split('utr'))[0].strip('_')
|
|
|
15 if coverage.has_key(name):
|
|
|
16 if score > coverage[name]:
|
|
|
17 coverage[name] = score
|
|
|
18 else:
|
|
|
19 coverage[name] = score
|
|
|
20 return coverage
|
|
|
21
|
|
|
22 def endBias(filename,utr5,utr3):
|
|
|
23 out = open(filename,'w')
|
|
|
24 for txpt in utr5.keys():
|
|
|
25 if utr3.has_key(txpt):
|
|
|
26 out.write('\t'.join([txpt,str(utr5[txpt]),str(utr3[txpt]),str(math.log((1+utr5[txpt])/(1+utr3[txpt]),2))])+'\n')
|
|
|
27 out.close()
|
|
|
28
|
|
|
29
|
|
|
30 utr5 = getCoverage(sys.argv[1])
|
|
|
31 utr3 = getCoverage(sys.argv[2])
|
|
|
32 endBias(sys.argv[3],utr5,utr3)
|
|
|
33
|
|
|
34 '''
|
|
|
35
|
|
|
36 utr5 = getCoverage('hmga2-utr5.coverage')
|
|
|
37 utr3 = getCoverage('hmga2-utr3.coverage')
|
|
|
38 logratio, cov5,cov3= endBias(utr5,utr3)
|
|
|
39 2**pylab.median(logratio.values())
|
|
|
40
|
|
|
41 log2utr5 = pylab.log2(pylab.array(cov5)+1)
|
|
|
42 log2utr3 = pylab.log2(pylab.array(cov3)+1)
|
|
|
43
|
|
|
44 pylab.plot(log2utr5,log2utr3,'bo')
|
|
|
45
|
|
|
46 pylab.show()
|
|
|
47
|
|
|
48 utr5 = getCoverage('control-utr5.coverage')
|
|
|
49 utr3 = getCoverage('control-utr3.coverage')
|
|
|
50 logratio, cov5,cov3= endBias(utr5,utr3)
|
|
|
51 2**pylab.median(logratio.values())
|
|
|
52 '''
|
