Hi,

I've got bisulfite-sequencing data for two differentiation stages. The raw data was mapped using Bismark. For each CpG site, the methylation ratio was marked as "A/B" (A methylated reads vs B unmethylated reads for this site).

Now I want to compare the overall methylation of a certain region. Assuming this region contains 3 CpG sites, the methylated/unmethylated ratio of each site in the two stages are 45/60, 4/10, 5/15 and 14/65, 3/9, 7/20 respectively.

I think there are two ways to calculate the overall methylation rate of the region:

**Method 1: total amount of methylated reads / total number of reads within the region**

methylation rate of stage 1 = (45+4+5)/(45+60+4+10+5+15) = **54/139** = 0.39

methylation rate of stage 2 = (14+3+7)/(14+65+3+9+7+20) = **24/118** = 0.20

**Method 2: average the methylation rate of each CpG site within the region**

methylation rate of stage 1 = (45/(45+60) + 4/(4+10) + 5/(5+15))/3 = 0.32

methylation rate of stage 2 = (14/(14+65) + 3/(3+9) + 7/(7+20))/3 = 0.22

For **method 1**, I can then calculate the confidence inteval and the P value to know whether **54 out of 139** is significantly different from **24 out of 118**.

However, the methylation rate calculated by **method 1** is biased toward the site with higher reads coverage.

**Method 2** seems to be more robust to indicate the actual methylation rate of the region.

But I don't know which statistics test should I use for method 2 to know the significance.

Please help.

Thanks in advance.
From context I gather that you lack biological replicates, which is unfortunate but still rather common. As you mentioned, "method 1", which is basically summing across CpGs and performing a Fisher's test is less than ideal, since it will miss interesting cases (using your earlier nomenclature, think of 50/0, 0/50 in sample1 and 25/25, 25/25 in sample 2, they're very different but this method would miss that).

Method 2 lends itself to a weighted paired t-test, which you can do in R. This method is, of course, also not ideal, since you give a weight to each CpG, but if coverage is drastically different between the two samples at a CpG then you start running into problems (there are ways around that, but they'd be moderately annoying).

What would likely be a nice method in this case would be to perform a Fisher's test on the individual CpGs between the samples and then spatially correct the p-values (the original paper from Benjamini & Heller is <a href="http://www.tandfonline.com/doi/abs/10.1198/016214507000000941#preview">here</a> but you can find an implementation in the <a href="http://www.bioconductor.org/packages/release/bioc/html/BiSeq.html">BiSeq bioconductor package</a>, which you could also just directly use to make your life easier since it's actually targeted at RRBS or other targeted BS-seq datasets). You could also use the SLIM model for p-value adjustment, such as is used in <a href="https://code.google.com/p/methylkit/">methylKit</a>. Somewhat similarly, you could also first use a smoothing method (see <a href="http://bioconductor.org/packages/release/bioc/html/bsseq.html">BSseq</a> as an example) and then do a paired t-test, though I wouldn't recomment that for your use case.

There are some additional possibilities, but I'd recommend just using BiSeq and being done with it.