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 here but you can find an implementation in the BiSeq bioconductor package, 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 methylKit. Somewhat similarly, you could also first use a smoothing method (see BSseq 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.