annotate methratio.patch @ 27:549ff234f35e draft

Uploaded
author eiriche
date Mon, 03 Dec 2012 02:52:27 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27
549ff234f35e Uploaded
eiriche
parents:
diff changeset
1 --- methratio.py.orig 2012-04-10 20:54:39.000000000 +0200
549ff234f35e Uploaded
eiriche
parents:
diff changeset
2 +++ methratio.py 2012-08-16 10:54:42.000000000 +0200
549ff234f35e Uploaded
eiriche
parents:
diff changeset
3 @@ -1,4 +1,5 @@
549ff234f35e Uploaded
eiriche
parents:
diff changeset
4 -import sys, time, os, array, optparse
549ff234f35e Uploaded
eiriche
parents:
diff changeset
5 +#!/usr/bin/python
549ff234f35e Uploaded
eiriche
parents:
diff changeset
6 +import sys, time, os, array, optparse, re
549ff234f35e Uploaded
eiriche
parents:
diff changeset
7 usage = "usage: %prog [options] BSMAP_MAPPING_FILES"
549ff234f35e Uploaded
eiriche
parents:
diff changeset
8 parser = optparse.OptionParser(usage=usage)
549ff234f35e Uploaded
eiriche
parents:
diff changeset
9
549ff234f35e Uploaded
eiriche
parents:
diff changeset
10 @@ -128,26 +129,32 @@
549ff234f35e Uploaded
eiriche
parents:
diff changeset
11
549ff234f35e Uploaded
eiriche
parents:
diff changeset
12 disp('writing %s ...' % options.outfile)
549ff234f35e Uploaded
eiriche
parents:
diff changeset
13 ss = {'C': '+', 'G': '-'}
549ff234f35e Uploaded
eiriche
parents:
diff changeset
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' }
549ff234f35e Uploaded
eiriche
parents:
diff changeset
15 fout = open(options.outfile, 'w')
549ff234f35e Uploaded
eiriche
parents:
diff changeset
16 -z95, z95sq = 1.96, 1.96 * 1.96
549ff234f35e Uploaded
eiriche
parents:
diff changeset
17 -fout.write('chr\tpos\tstrand\tcontext\tratio\ttotal_C\tmethy_C\tCI_lower\tCI_upper\n')
549ff234f35e Uploaded
eiriche
parents:
diff changeset
18 +
549ff234f35e Uploaded
eiriche
parents:
diff changeset
19 +fout.write('chr\tpos\tstrand\tcontext\tcoverage\tmethylated_C\tperc_methy_C\n')
549ff234f35e Uploaded
eiriche
parents:
diff changeset
20 nc, nd, dep0 = 0, 0, options.min_depth
549ff234f35e Uploaded
eiriche
parents:
diff changeset
21 for cr in sorted(depth.keys()):
549ff234f35e Uploaded
eiriche
parents:
diff changeset
22 depthcr, methcr, refcr = depth[cr], meth[cr], ref[cr]
549ff234f35e Uploaded
eiriche
parents:
diff changeset
23 for i, d in enumerate(depthcr):
549ff234f35e Uploaded
eiriche
parents:
diff changeset
24 + if i < 2: continue
549ff234f35e Uploaded
eiriche
parents:
diff changeset
25 if d < dep0: continue
549ff234f35e Uploaded
eiriche
parents:
diff changeset
26 nc += 1
549ff234f35e Uploaded
eiriche
parents:
diff changeset
27 nd += d
549ff234f35e Uploaded
eiriche
parents:
diff changeset
28 m = methcr[i]
549ff234f35e Uploaded
eiriche
parents:
diff changeset
29 if m == 0 and not options.meth0: continue
549ff234f35e Uploaded
eiriche
parents:
diff changeset
30 ratio = float(m) / d
549ff234f35e Uploaded
eiriche
parents:
diff changeset
31 - seq = refcr[i-2:i+3]
549ff234f35e Uploaded
eiriche
parents:
diff changeset
32 strand = ss[refcr[i]]
549ff234f35e Uploaded
eiriche
parents:
diff changeset
33 - pmid = ratio + z95sq / (2 * d)
549ff234f35e Uploaded
eiriche
parents:
diff changeset
34 - sd = z95 * ((ratio*(1-ratio)/d + z95sq/(4*d*d)) ** 0.5)
549ff234f35e Uploaded
eiriche
parents:
diff changeset
35 - norminator = 1 + z95sq / d
549ff234f35e Uploaded
eiriche
parents:
diff changeset
36 - CIl, CIu = (pmid - sd) / norminator, (pmid + sd) / norminator
549ff234f35e Uploaded
eiriche
parents:
diff changeset
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))
549ff234f35e Uploaded
eiriche
parents:
diff changeset
38 + if strand == '-':
549ff234f35e Uploaded
eiriche
parents:
diff changeset
39 + seq = refcr[i-2:i+1]
549ff234f35e Uploaded
eiriche
parents:
diff changeset
40 + else:
549ff234f35e Uploaded
eiriche
parents:
diff changeset
41 + seq = refcr[i:i+3]
549ff234f35e Uploaded
eiriche
parents:
diff changeset
42 +
549ff234f35e Uploaded
eiriche
parents:
diff changeset
43 + for k, v in trint.items():
549ff234f35e Uploaded
eiriche
parents:
diff changeset
44 + if re.match(k,seq):
549ff234f35e Uploaded
eiriche
parents:
diff changeset
45 + context=v
549ff234f35e Uploaded
eiriche
parents:
diff changeset
46 +
549ff234f35e Uploaded
eiriche
parents:
diff changeset
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))
549ff234f35e Uploaded
eiriche
parents:
diff changeset
48
549ff234f35e Uploaded
eiriche
parents:
diff changeset
49 fout.close()
549ff234f35e Uploaded
eiriche
parents:
diff changeset
50 disp('done.')