|
27
|
1 --- methratio.py.orig 2012-04-10 20:54:39.000000000 +0200
|
|
|
2 +++ methratio.py 2012-08-16 10:54:42.000000000 +0200
|
|
|
3 @@ -1,4 +1,5 @@
|
|
|
4 -import sys, time, os, array, optparse
|
|
|
5 +#!/usr/bin/python
|
|
|
6 +import sys, time, os, array, optparse, re
|
|
|
7 usage = "usage: %prog [options] BSMAP_MAPPING_FILES"
|
|
|
8 parser = optparse.OptionParser(usage=usage)
|
|
|
9
|
|
|
10 @@ -128,26 +129,32 @@
|
|
|
11
|
|
|
12 disp('writing %s ...' % options.outfile)
|
|
|
13 ss = {'C': '+', 'G': '-'}
|
|
|
14 +trint = {'CG[ACGT]$':'CPG','C[ACT]G$':'CHG', r"C[ACT][ACT]":'CHH', r"[ACGT]CG":'CPG',r"C[AGT]G":'CHG',r"[AGT][AGT]G":'CHH' }
|
|
|
15 fout = open(options.outfile, 'w')
|
|
|
16 -z95, z95sq = 1.96, 1.96 * 1.96
|
|
|
17 -fout.write('chr\tpos\tstrand\tcontext\tratio\ttotal_C\tmethy_C\tCI_lower\tCI_upper\n')
|
|
|
18 +
|
|
|
19 +fout.write('chr\tpos\tstrand\tcontext\tcoverage\tmethylated_C\tperc_methy_C\n')
|
|
|
20 nc, nd, dep0 = 0, 0, options.min_depth
|
|
|
21 for cr in sorted(depth.keys()):
|
|
|
22 depthcr, methcr, refcr = depth[cr], meth[cr], ref[cr]
|
|
|
23 for i, d in enumerate(depthcr):
|
|
|
24 + if i < 2: continue
|
|
|
25 if d < dep0: continue
|
|
|
26 nc += 1
|
|
|
27 nd += d
|
|
|
28 m = methcr[i]
|
|
|
29 if m == 0 and not options.meth0: continue
|
|
|
30 ratio = float(m) / d
|
|
|
31 - seq = refcr[i-2:i+3]
|
|
|
32 strand = ss[refcr[i]]
|
|
|
33 - pmid = ratio + z95sq / (2 * d)
|
|
|
34 - sd = z95 * ((ratio*(1-ratio)/d + z95sq/(4*d*d)) ** 0.5)
|
|
|
35 - norminator = 1 + z95sq / d
|
|
|
36 - CIl, CIu = (pmid - sd) / norminator, (pmid + sd) / norminator
|
|
|
37 - fout.write('%s\t%d\t%c\t%s\t%.3f\t%d\t%d\t%.3f\t%.3f\n' % (cr, i+1, strand, seq, ratio, d, m, CIl, CIu))
|
|
|
38 + if strand == '-':
|
|
|
39 + seq = refcr[i-2:i+1]
|
|
|
40 + else:
|
|
|
41 + seq = refcr[i:i+3]
|
|
|
42 +
|
|
|
43 + for k, v in trint.items():
|
|
|
44 + if re.match(k,seq):
|
|
|
45 + context=v
|
|
|
46 +
|
|
|
47 + fout.write('%s\t%d\t%c\t%s\t%d\t%d\t%.2f\n' % (cr, i+1, strand, context, d, m, m*100/d))
|
|
|
48
|
|
|
49 fout.close()
|
|
|
50 disp('done.')
|
