Mercurial > repos > xuebing > sharplabtool

comparison mytools/endbias.py @ 9:87eb5c5ddfe9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 20:01:43 -0500
parents f0dc65e7f6c0
children
comparison
equal deleted inserted replaced
8:361ec1c0479d 9:87eb5c5ddfe9
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 '''