---
title: "06 - Differential methylation (oil vs. control)"
output: html_document
---

Test for CpGs (DMLs) and 1 kb regions (DMRs) that differ in methylation between oil-exposed
and control gill, using methylKit on the deduplicated, strand-merged CpG calls from step 5.

Decisions carried over from steps 4–5 (see `plan.md`):

- **HB30 is dropped.** The comparison is 3 control (NB3, NB6, NB11) vs. 2 oil (HB2, HB16).
- **Deduplicated calls are used.** Duplicates agree with the deduplicated calls but would inflate test confidence.
- **Global methylation differs a lot between samples in the same group** (36–72.5% among controls),
  almost certainly for technical reasons. Check the PCA and the permutation test below before trusting any DML.

- Input: `data/methylation/dedup/SRR*.CpG_report.merged_CpG_evidence.cov.gz` (from `05-methylation-calling`)
- Output: `output/06-differential-methylation/`

Uses methylKit, data.table, ggplot2 and pheatmap from the lab R container.
Run on a compute node; the tests use `n_cores` cores.

```{r setup, include=FALSE}
# Paths are relative to code/, so chunks work the same run one at a time or knitted
knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
library(methylKit)
library(data.table)
library(ggplot2)
library(pheatmap)
```

# Settings

```{r settings}
out_dir <- "../output/06-differential-methylation"
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

n_cores   <- 8
min_covs  <- c(5, 10)   # DML coverage cutoffs, required in every sample
max_perc  <- 99.9       # drop CpGs above this coverage percentile in each sample (PCR / repeat pileups)
q_cutoff  <- 0.01
diff_cutoff <- 25       # percentage points

tile_size       <- 1000
tile_min_cpgs   <- 3    # covered CpGs a tile needs in a sample
tile_min_cov    <- 10   # total reads a tile needs in a sample
```

# Samples

```{r samples}
sample_sheet <- fread("../data/sample_sheet.csv")
samples <- sample_sheet[sample != "HB30"][order(treatment, sample)]
samples[, treatment_code := as.integer(treatment == "oil")]   # methylKit: 0 = control, 1 = oil
samples[, file := file.path("../data/methylation/dedup",
                            paste0(run, ".CpG_report.merged_CpG_evidence.cov.gz"))]
stopifnot(all(file.exists(samples$file)), nrow(samples) == 5)
samples[, .(run, sample, treatment, treatment_code)]
```

# Read CpG calls

The merged files are in Bismark coverage format (chromosome, start, end, % methylated, count methylated, count unmethylated),
with the two strands of each CpG already added together, so no de-stranding is needed.
Read everything at 1× so the same objects serve both the single-CpG tests and the tiles.

```{r read}
raw <- methRead(
  location   = as.list(samples$file),
  sample.id  = as.list(samples$sample),
  assembly   = "C_virginica-3.0",
  treatment  = samples$treatment_code,
  context    = "CpG",
  pipeline   = "bismarkCoverage",
  mincov     = 1
)
raw
```

# Single CpGs

For each coverage cutoff: drop CpGs below it and above the 99.9th coverage percentile in each sample,
equalize median coverage across samples, and keep CpGs that pass in all five.

```{r unite}
united <- lapply(setNames(min_covs, paste0(min_covs, "x")), function(min_cov) {
  filtered <- filterByCoverage(raw, lo.count = min_cov, hi.perc = max_perc)
  unite(normalizeCoverage(filtered, method = "median"), destrand = FALSE, mc.cores = n_cores)
})
sapply(united, nrow)
```

## Global methylation

Mean methylation per sample over the CpGs shared by all five.

```{r global-meth}
global <- rbindlist(lapply(names(united), function(cov) {
  pm <- percMethylation(united[[cov]])
  data.table(min_cov = cov, sample = colnames(pm), mean_meth = round(colMeans(pm), 1),
             median_meth = round(apply(pm, 2, median), 1))
}))
global <- merge(samples[, .(sample, treatment)], global, by = "sample")
setorder(global, min_cov, treatment, sample)
fwrite(global, file.path(out_dir, "global-methylation.csv"))
global
```

## Do samples group by treatment?

PCA and clustering on per-CpG methylation (5× set). If samples separate by global methylation level
(NB6 low, NB3 high) rather than by treatment, most of the variation is technical.

```{r pca, fig.width = 6, fig.height = 5}
pm <- percMethylation(united[["5x"]])
# Per-CpG methylation of the 5x set, used by step 7 to describe the methylation landscape
fwrite(cbind(as.data.table(getData(united[["5x"]])[, c("chr", "start", "end")]), round(pm, 1)),
       file.path(out_dir, "cpg-methylation-5x.csv.gz"))
pca <- prcomp(t(pm), center = TRUE, scale. = FALSE)
pct <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1)
pca_dt <- data.table(sample = rownames(pca$x), PC1 = pca$x[, 1], PC2 = pca$x[, 2])
pca_dt <- merge(pca_dt, samples[, .(sample, treatment)], by = "sample")
pca_dt <- merge(pca_dt, global[min_cov == "5x", .(sample, mean_meth)], by = "sample")

p <- ggplot(pca_dt, aes(PC1, PC2, colour = treatment, label = paste0(sample, " (", mean_meth, "%)"))) +
  geom_point(size = 3) +
  geom_text(vjust = -1, size = 3.5, show.legend = FALSE) +
  scale_x_continuous(expand = expansion(mult = 0.25)) +
  scale_y_continuous(expand = expansion(mult = 0.25)) +
  labs(x = paste0("PC1 (", pct[1], "%)"), y = paste0("PC2 (", pct[2], "%)"), colour = "Treatment",
       title = "PCA of per-CpG methylation (>=5x in all samples)",
       subtitle = "Labels show each sample's mean methylation") +
  theme_bw()
ggsave(file.path(out_dir, "pca-5x.png"), p, width = 6, height = 5, dpi = 150)
p
```

```{r correlation, fig.width = 6, fig.height = 5}
cor_mat <- cor(pm)
fwrite(as.data.table(cor_mat, keep.rownames = "sample"), file.path(out_dir, "sample-correlation-5x.csv"))
annotation <- data.frame(treatment = samples$treatment, row.names = samples$sample)
pheatmap(cor_mat, annotation_row = annotation, annotation_col = annotation,
         display_numbers = TRUE, number_format = "%.2f",
         clustering_method = "ward.D2", main = "Sample correlation (>=5x CpGs)",
         filename = file.path(out_dir, "sample-correlation-5x.png"), width = 6, height = 5)
pheatmap(cor_mat, annotation_row = annotation, annotation_col = annotation,
         display_numbers = TRUE, number_format = "%.2f",
         clustering_method = "ward.D2", main = "Sample correlation (>=5x CpGs)")
```

## Test each CpG

Logistic regression with overdispersion correction (`overdispersion = "MN"`, `test = "Chisq"`),
which methylKit recommends for small groups with biological replicates.
A DML needs q < `q_cutoff` and a difference of at least `diff_cutoff` percentage points (oil − control).

```{r dml}
bed_out <- function(x, file) {
  # BED is 0-based, half-open; methylKit coordinates are 1-based
  fwrite(data.table(chr = x$chr, start = x$start - 1, end = x$end,
                    name = paste0("meth.diff=", round(x$meth.diff, 1), ";q=", signif(x$qvalue, 3))),
         file, sep = "\t", col.names = FALSE)
}

diff_cpg <- lapply(united, calculateDiffMeth, overdispersion = "MN", test = "Chisq", mc.cores = n_cores)

dml_counts <- rbindlist(lapply(names(diff_cpg), function(cov) {
  d <- diff_cpg[[cov]]
  all_tested <- getData(d)
  dml <- getData(getMethylDiff(d, difference = diff_cutoff, qvalue = q_cutoff, type = "all"))

  fwrite(all_tested, file.path(out_dir, paste0("cpg-tests-", cov, ".csv.gz")))
  fwrite(dml, file.path(out_dir, paste0("dml-", cov, ".csv")))
  bed_out(dml, file.path(out_dir, paste0("dml-", cov, ".bed")))
  bed_out(all_tested, file.path(out_dir, paste0("tested-cpgs-", cov, ".bed")))   # background for step 7

  data.table(min_cov = cov, CpGs_tested = nrow(all_tested),
             DMLs = nrow(dml), hyper_in_oil = sum(dml$meth.diff > 0), hypo_in_oil = sum(dml$meth.diff < 0),
             pct_of_tested = round(100 * nrow(dml) / nrow(all_tested), 3))
}))
fwrite(dml_counts, file.path(out_dir, "dml-counts.csv"))
dml_counts
```

```{r pvalue-hist, fig.width = 8, fig.height = 3.5}
pvals <- rbindlist(lapply(names(diff_cpg), function(cov) data.table(min_cov = cov, pvalue = getData(diff_cpg[[cov]])$pvalue)))
p <- ggplot(pvals, aes(pvalue)) +
  geom_histogram(breaks = seq(0, 1, 0.02)) +
  facet_wrap(~ min_cov, scales = "free_y") +
  labs(x = "p-value", y = "CpGs", title = "p-value distributions",
       subtitle = "A real effect shows as a spike near 0 on top of a flat background") +
  theme_bw()
ggsave(file.path(out_dir, "pvalue-histograms.png"), p, width = 8, height = 3.5, dpi = 150)
p
```

## Permutation check

With 3 vs. 2 samples there are only 10 ways to pick which two samples are "oil".
Re-run the test for every split and count DMLs. If the real split (HB2 + HB16) does not stand out from the other nine,
the DMLs above are no more than the technical noise between samples produces by chance.
Uses the 5× set. With only 10 splits, the smallest possible permutation p-value is 0.1.

Each split's DMLs are also written to `permutations/dml-5x_<oil samples>.bed`, so step 7 can check
whether the real list's genomic distribution differs from what random splits give.

```{r permutation}
perm_dir <- file.path(out_dir, "permutations")
dir.create(perm_dir, showWarnings = FALSE)
split_name <- function(oil_pair) paste(oil_pair, collapse = "+")

perm_counts <- rbindlist(lapply(combn(samples$sample, 2, simplify = FALSE), function(oil_pair) {
  u <- united[["5x"]]
  u@treatment <- as.integer(u@sample.ids %in% oil_pair)
  d <- calculateDiffMeth(u, overdispersion = "MN", test = "Chisq", mc.cores = n_cores)
  dml <- getData(getMethylDiff(d, difference = diff_cutoff, qvalue = q_cutoff, type = "all"))
  bed_out(dml, file.path(perm_dir, paste0("dml-5x_", split_name(oil_pair), ".bed")))
  data.table(oil_samples = paste(oil_pair, collapse = " + "),
             real_split  = setequal(oil_pair, samples[treatment == "oil", sample]),
             DMLs = nrow(dml), hyper = sum(dml$meth.diff > 0), hypo = sum(dml$meth.diff < 0))
}))
setorder(perm_counts, -DMLs)
perm_counts[, rank := seq_len(.N)]
fwrite(perm_counts, file.path(out_dir, "permutation-dml-counts-5x.csv"))

real_rank <- perm_counts[real_split == TRUE, rank]
cat("Real split ranks", real_rank, "of", nrow(perm_counts),
    "by DML count (permutation p =", real_rank / nrow(perm_counts), ")\n")
perm_counts
```

# Regions (1 kb tiles)

Pools CpGs in non-overlapping 1 kb windows, which recovers much of the sparse data that single-CpG
tests throw away. A tile counts in a sample if it has at least `tile_min_cpgs` covered CpGs and
`tile_min_cov` reads in total; it is tested if it passes in all five samples.

```{r tiles}
tiles <- tileMethylCounts(raw, win.size = tile_size, step.size = tile_size,
                          cov.bases = tile_min_cpgs, mc.cores = n_cores)
tiles <- filterByCoverage(tiles, lo.count = tile_min_cov, hi.perc = max_perc)
tiles_united <- unite(normalizeCoverage(tiles, method = "median"), destrand = FALSE, mc.cores = n_cores)
diff_tiles <- calculateDiffMeth(tiles_united, overdispersion = "MN", test = "Chisq", mc.cores = n_cores)

tiles_tested <- getData(diff_tiles)
dmr <- getData(getMethylDiff(diff_tiles, difference = diff_cutoff, qvalue = q_cutoff, type = "all"))

fwrite(tiles_tested, file.path(out_dir, "tile-tests-1kb.csv.gz"))
fwrite(dmr, file.path(out_dir, "dmr-1kb.csv"))
bed_out(dmr, file.path(out_dir, "dmr-1kb.bed"))
bed_out(tiles_tested, file.path(out_dir, "tested-tiles-1kb.bed"))   # background for step 7

dmr_counts <- data.table(tiles_tested = nrow(tiles_tested), DMRs = nrow(dmr),
                         hyper_in_oil = sum(dmr$meth.diff > 0), hypo_in_oil = sum(dmr$meth.diff < 0))
fwrite(dmr_counts, file.path(out_dir, "dmr-counts.csv"))
dmr_counts
```

## Permutation check for regions

Same idea as for single CpGs: re-test the tiles for all 10 ways of picking two "oil" samples
and see where the real split (HB2 + HB16) ranks by DMR count.
Each split's DMRs go to `permutations/dmr-1kb_<oil samples>.bed`.

```{r tile-permutation}
tile_perm_counts <- rbindlist(lapply(combn(samples$sample, 2, simplify = FALSE), function(oil_pair) {
  u <- tiles_united
  u@treatment <- as.integer(u@sample.ids %in% oil_pair)
  d <- calculateDiffMeth(u, overdispersion = "MN", test = "Chisq", mc.cores = n_cores)
  hits <- getData(getMethylDiff(d, difference = diff_cutoff, qvalue = q_cutoff, type = "all"))
  bed_out(hits, file.path(perm_dir, paste0("dmr-1kb_", split_name(oil_pair), ".bed")))
  data.table(oil_samples = paste(oil_pair, collapse = " + "),
             real_split  = setequal(oil_pair, samples[treatment == "oil", sample]),
             DMRs = nrow(hits), hyper = sum(hits$meth.diff > 0), hypo = sum(hits$meth.diff < 0))
}))
setorder(tile_perm_counts, -DMRs)
tile_perm_counts[, rank := seq_len(.N)]
fwrite(tile_perm_counts, file.path(out_dir, "permutation-dmr-counts-1kb.csv"))

real_rank <- tile_perm_counts[real_split == TRUE, rank]
cat("Real split ranks", real_rank, "of", nrow(tile_perm_counts),
    "by DMR count (permutation p =", real_rank / nrow(tile_perm_counts), ")\n")
tile_perm_counts
```

# Summary

```{r session}
sessionInfo()
```

Things to check before moving on to annotation (step 7):

- **PCA / correlation:** do HB2 and HB16 sit together, apart from the controls? Or do samples line up by mean methylation?
- **p-value histograms:** a spike near 0 means real differences; a flat or hump-shaped histogram means little signal.
- **Permutation ranks (CpGs and tiles):** the real split should have more hits than most of the other nine. Rank 1 of 10 is the best possible (p = 0.1).
- **Hyper vs. hypo balance:** a strong skew in one direction can come from the global methylation differences rather than oil.