---
title: "05 - Differential methylation (DMLs): main contrasts"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(eval = FALSE, message = FALSE, warning = FALSE)
library(methylKit); library(data.table)
```

# Contrast helper

`calculateDiffMeth` uses logistic regression with overdispersion correction.
meth.diff = treatment (group 1) - control (group 0). DML threshold:
**q < 0.01 and |meth.diff| >= 15%**.

```{r helper}
run_contrast <- function(filt, meta, idx, treat, label,
                         qval = 0.01, mdiff = 15, mc.cores = 4) {
  sub <- reorganize(filt, sample.ids = as.list(meta$sampleID[idx]), treatment = treat)
  m   <- unite(sub, destrand = FALSE)
  dm  <- calculateDiffMeth(m, mc.cores = mc.cores)
  all <- as.data.table(getData(dm))
  dml <- all[qvalue < qval & abs(meth.diff) >= mdiff]
  fwrite(all, paste0("diffmeth_", label, "_all.csv.gz"))
  fwrite(dml, paste0("dml_", label, ".csv"))
  saveRDS(dm, paste0("res_", label, ".rds"))
  data.table(contrast = label, n_tested = nrow(all), n_DML = nrow(dml),
             hyper = sum(dml$meth.diff > 0), hypo = sum(dml$meth.diff < 0))
}
```

# Run the three contrasts

```{r contrasts}
filt <- readRDS("methylRaw_filtered.rds")
meta <- fread("sample_metadata.csv")

# (i) pH within diploids: samples 1-12, treat 0 = high, 1 = low
r1 <- run_contrast(filt, meta, idx = 1:12,
                   treat = c(rep(0,6), rep(1,6)), label = "pH_within_2N")

# (ii) pH within triploids: samples 13-24
r2 <- run_contrast(filt, meta, idx = 13:24,
                   treat = c(rep(0,6), rep(1,6)), label = "pH_within_3N")

# (iii) ploidy main effect: all 24, 0 = 2N, 1 = 3N (uses the pre-united object)
meth <- readRDS("meth_united.rds")
meth@treatment <- c(rep(0,12), rep(1,12))
dm3 <- calculateDiffMeth(meth, mc.cores = 4)
all3 <- as.data.table(getData(dm3))
dml3 <- all3[qvalue < 0.01 & abs(meth.diff) >= 15]
fwrite(all3, "diffmeth_ploidy_main_all.csv.gz"); fwrite(dml3, "dml_ploidy_main.csv")
r3 <- data.table(contrast = "ploidy_main", n_tested = nrow(all3), n_DML = nrow(dml3),
                 hyper = sum(dml3$meth.diff > 0), hypo = sum(dml3$meth.diff < 0))

dml_summary <- rbind(r1, r2, r3)
fwrite(dml_summary, "dml_summary.csv")
dml_summary[]
```

Expected: pH within 2N ~16,036 DMLs (mostly hyper); pH within 3N ~25,986
(mostly hypo - opposite direction); ploidy main ~8,216. The direction flip
between ploidies is the key result.