Mercurial > repos > weilong-guo > bs_seeker2

comparison BSseeker2/bs_seeker2-call_methylation.py @ 1:8b26adf64adc draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
V2.0.5
author weilong-guo
date Tue, 05 Nov 2013 01:55:39 -0500
parents e6df770c0e58
children
comparison
equal deleted inserted replaced
0:e6df770c0e58 1:8b26adf64adc
1 #!/usr/bin/python 1 #!/usr/bin/env python
2 2
3 from optparse import OptionParser, OptionGroup 3 from optparse import OptionParser, OptionGroup
4 from bs_utils.utils import * 4 from bs_utils.utils import *
5 5
6 try : 6 try :
69 (options, args) = parser.parse_args() 69 (options, args) = parser.parse_args()
70 70
71 71
72 # if no options were given by the user, print help and exit 72 # if no options were given by the user, print help and exit
73 if len(sys.argv) == 1: 73 if len(sys.argv) == 1:
74 print parser.print_help() 74 parser.print_help()
75 exit(0) 75 exit(0)
76 76
77 if options.version : 77 if options.version :
78 show_version() 78 show_version()
79 exit (-1) 79 exit (-1)
131 ATCG_rev[n] = 0 131 ATCG_rev[n] = 0
132 132
133 nuc, context, subcontext = context_calling(chrom_seq, col.pos) 133 nuc, context, subcontext = context_calling(chrom_seq, col.pos)
134 total_reads = 0 134 total_reads = 0
135 135
136
137
138 for pr in col.pileups: 136 for pr in col.pileups:
139 # print pr 137 # print pr
140 if (not pr.indel) : # skip indels 138 if (not pr.indel) : # skip indels
141 #if ( (options.RM_SX) and (pr.alignment.tags[1][1] == 1) ): 139 #if ( (options.RM_SX) and (pr.alignment.tags[1][1] == 1) ):
142 ##=== Fixed error reported by Roberto 140 ##=== Fixed error reported by Roberto
150 148
151 if pr.qpos >= len(pr.alignment.seq): 149 if pr.qpos >= len(pr.alignment.seq):
152 print 'WARNING: read %s has an invalid alignment. Discarding.. ' % pr.alignment.qname 150 print 'WARNING: read %s has an invalid alignment. Discarding.. ' % pr.alignment.qname
153 continue 151 continue
154 read_nuc = pr.alignment.seq[pr.qpos] 152 read_nuc = pr.alignment.seq[pr.qpos]
155 # print "read_nuc=", read_nuc
156 if pr.alignment.is_reverse: 153 if pr.alignment.is_reverse:
157 ATCG_rev[read_nuc] += 1 154 ATCG_rev[read_nuc] += 1
158 else: 155 else:
159 ATCG_fwd[read_nuc] += 1 156 ATCG_fwd[read_nuc] += 1
160 157
192 # print all_cytosines 189 # print all_cytosines
193 if nuc == 'C': 190 if nuc == 'C':
194 wiggle.write('%d\t%f\n' % (pos, meth_level)) 191 wiggle.write('%d\t%f\n' % (pos, meth_level))
195 else : 192 else :
196 wiggle.write('%d\t-%f\n' % (pos, meth_level)) 193 wiggle.write('%d\t-%f\n' % (pos, meth_level))
197 CGmap.write('%(chrom)s\t%(nuc)s\t%(pos)d\t%(context)s\t%(subcontext)s\t%(meth_level_string)s\t%(meth_cytosines)s\t%(all_cytosines)s\n' % locals()) 194
195 CGmap.write('%(chrom)s\t%(nuc)s\t%(pos)d\t%(context)s\t%(subcontext)s\t%(meth_level_string)s\t%(meth_cytosines)s\t%(all_cytosines)s\n' % locals())
198 ATCGmap.close() 196 ATCGmap.close()
199 CGmap.close() 197 CGmap.close()
200 wiggle.close() 198 wiggle.close()
201 199
202 logm('Wiggle: %s'% wiggle_fname) 200 logm('Wiggle: %s'% wiggle_fname)