---
title: "06 - Ploidy x pH interaction test (per-CpG GLM)"
output: html_document
---

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

methylKit has no interaction term, so we fit a per-CpG binomial GLM
`cbind(methylated, unmethylated) ~ ploidy * pH` on the common CpGs and extract
the Wald p-value for the interaction coefficient (`ploidy3N:pHlow`),
BH-adjusted. Effect size = difference of level-differences from group weighted
means: (3N low - 3N high) - (2N low - 2N high). Threshold **q < 0.05 and
|effect| >= 15%**.

```{r prep}
meth <- readRDS("meth_united.rds")
meta <- fread("sample_metadata.csv")
d <- as.data.table(getData(meth))
cs  <- as.matrix(d[, grep("^numCs",     names(d)), with = FALSE])   # methylated
cov <- as.matrix(d[, grep("^coverage",  names(d)), with = FALSE])
ts  <- cov - cs                                                     # unmethylated
ploidy <- factor(meta$ploidy, levels = c("2N","3N"))
pH     <- factor(meta$pH,     levels = c("high","low"))
```

```{r effect_size}
grp <- interaction(ploidy, pH)
gm <- sapply(levels(grp), function(g){
  k <- which(grp == g); tc <- rowSums(cs[,k,drop=FALSE]); tt <- rowSums(ts[,k,drop=FALSE])
  100 * tc / (tc + tt)
})
d2N <- gm[,"2N.low"] - gm[,"2N.high"]
d3N <- gm[,"3N.low"] - gm[,"3N.high"]
int_effect <- d3N - d2N
```

```{r glm_parallel}
fit_one <- function(i){
  y <- cbind(cs[i,], ts[i,])
  fit <- tryCatch(suppressWarnings(glm(y ~ ploidy * pH, family = binomial())),
                  error = function(e) NULL)
  if (is.null(fit)) return(NA_real_)
  cf <- summary(fit)$coefficients
  ir <- grep(":", rownames(cf), value = TRUE)   # the interaction term
  if (length(ir) == 0) return(NA_real_)
  cf[ir[1], "Pr(>|z|)"]
}
n <- nrow(cs)
chunks <- split(seq_len(n), cut(seq_len(n), 12, labels = FALSE))
pv <- unlist(mclapply(chunks, function(ix) vapply(ix, fit_one, numeric(1)),
                      mc.cores = 6), use.names = FALSE)   # ~4-5 min
qv <- p.adjust(pv, "BH")

res <- cbind(d[, .(chr, start, end, strand)],
             data.table(d2N = d2N, d3N = d3N, int_effect = int_effect,
                        pvalue = pv, qvalue = qv))
res[, sig := qvalue < 0.05 & abs(int_effect) >= 15]
fwrite(res, "diffmeth_interaction_all.csv.gz")
fwrite(res[sig == TRUE], "dml_interaction.csv")
sum(res$sig)   # expect ~17,053
```

```{r classify}
sigres <- res[sig == TRUE]
sigres[, pattern := fifelse(d2N > 0 & d3N < 0, "2N up / 3N down",
                     fifelse(d2N < 0 & d3N > 0, "2N down / 3N up",
                     fifelse(sign(d2N)==sign(d3N) & abs(d3N)>abs(d2N), "3N amplified",
                     fifelse(sign(d2N)==sign(d3N), "2N amplified", "other"))))]
sigres[, .N, by = pattern][order(-N)]
# ~96% are directional reversals (opposite sign between ploidies)
```