Mercurial > repos > eiriche > bsmap
comparison methratio.patch @ 27:549ff234f35e draft
Uploaded
author | eiriche |
---|---|
date | Mon, 03 Dec 2012 02:52:27 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
26:701d6bb0e28d | 27:549ff234f35e |
---|---|
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.') |