---
title: "04 - Global methylation comparison"
output: html_document
---

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

# Weighted global methylation per sample

Weighted global methylation = 100 * sum(methylated) / sum(coverage) over the
filtered CpGs (weights each CpG by its coverage).

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

gm <- rbindlist(lapply(seq_along(filt), function(i) {
  d <- getData(filt[[i]])
  data.table(sampleID = filt[[i]]@sample.id,
             global_meth = 100 * sum(d$numCs) / sum(d$coverage))
}))
gm <- merge(gm, meta[, .(sampleID, ploidy, pH, group)], by = "sampleID")
fwrite(gm, "global_methylation.csv")
```

# Two-way ANOVA: the interaction

```{r anova}
gm[, ploidy := factor(ploidy, levels = c("2N","3N"))]
gm[, pH     := factor(pH,     levels = c("high","low"))]
fit <- aov(global_meth ~ ploidy * pH, data = gm)
summary(fit)
# Expect: interaction p ~ 0.018 (significant); ploidy p ~ 0.13; pH p ~ 0.63
gm[, .(mean = mean(global_meth), sd = sd(global_meth)), by = group]
# Crossover: low pH RAISES methylation in 2N but LOWERS it in 3N
```

```{r plot}
grp_cols <- c("2N_high"="#4E79A7","2N_low"="#A0CBE8","3N_high"="#E15759","3N_low"="#FF9D9A")
agg <- gm[, .(m = mean(global_meth), se = sd(global_meth)/sqrt(.N)), by = .(ploidy, pH)]
p <- ggplot(agg, aes(pH, m, color = ploidy, group = ploidy)) +
  geom_line() + geom_point(size = 3) +
  geom_errorbar(aes(ymin = m - se, ymax = m + se), width = 0.1) +
  labs(y = "Weighted global methylation (%)",
       subtitle = "ANOVA interaction p = 0.018") + theme_bw()
ggsave("global_methylation.png", p, width = 7, height = 5, dpi = 200)
```