I put a script here that uses the UCSC mysql server to grab the exons for your gene and then calculate the coverage from the cpgIslandExt table

Usage is

$ python exon-cpg.py [genome] [gene]

e.g.

$ python exon-cpg.py hg19 gata3

gives:

exon-start  exon-end    coverage    cpg:pct,...
8096666 8096854 1.0 CpG::60.1
8097249 8097859 1.0 CpG::60.1
8100267 8100804 0.715083798883  CpG::71.4
8105955 8106101 0   
8111435 8111561 0   
8115701 8117164 0

where the final column is a comma-delimited list of cpg-name:pct-gc

You can verify this is correct by looking at Gata3 in the browser

1. Find your specific gene in NCBI's Gene database (e.g. GNAS complex locus) and open gene description page.
2. Find NCBI Reference Sequence (RefSeq) section, find Genomic subsection and open gene DNA sequence in GenBank format from Nucleotide DB GNAS example.
3. Find GI in the Version section (281182534 GNAS)
4. Find the Features section and save all coordinates of exons.
5. Get all exons sequences with Biopython and my code example.
6. Compute percentage of CG pairs with this

Python code:

exons = [exon_seq1, exon_seq2, exon_seq3]
for exon in exons:
    exon = exon.lower()
    CpG_number = exon.count("cg")
    CpG_proc = CpG_number * 2 * 100.0 / len(exon)