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.') |