---
title: "02 - Load into methylKit, QC, filter"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(eval = FALSE, message = FALSE, warning = FALSE)
library(methylKit); library(data.table); library(ggplot2)
# On osx-arm64 conda R, set before anything that may compile:
# Sys.setenv(R_MAKEVARS_USER = file.path(getwd(), ".rconfig/Makevars"))
```

# Read coverage files

```{r read}
meta <- fread("sample_metadata.csv")

obj <- methRead(
  as.list(meta$cov_file),
  sample.id = as.list(meta$sampleID),
  assembly  = "cgigas_roslin_v1",
  treatment = as.integer(meta$group_code),
  context   = "CpG",
  pipeline  = "bismarkCoverage",
  header    = FALSE,
  mincov    = 1
)
saveRDS(obj, "meta.rds")
```

# Per-sample QC

```{r qc}
qc <- rbindlist(lapply(seq_along(obj), function(i) {
  d <- getData(obj[[i]]); cov <- d$coverage; m <- 100 * d$numCs / cov
  data.table(sampleID = obj[[i]]@sample.id, nCpG = nrow(d),
             mean_cov = mean(cov), median_cov = median(cov),
             pct_cov_ge5  = 100 * mean(cov >= 5),
             pct_cov_ge10 = 100 * mean(cov >= 10),
             mean_meth = mean(m), median_meth = median(m))
}))
qc <- merge(qc, meta[, .(sampleID, ploidy, pH, group)], by = "sampleID")
fwrite(qc, "qc_per_sample.csv")
qc[]
```

Expected: ~9.5-10.9 M CpGs per sample, mean coverage 6.5-12.5x, 54-74% of CpGs
at >=5x, mean methylation ~10-12%. The methylation distribution is bimodal
(mosaic), typical of invertebrates (~83% of CpGs <=10% methylated, ~5% >=90%).

```{r qc_figure}
grp_cols <- c("2N_high"="#4E79A7","2N_low"="#A0CBE8","3N_high"="#E15759","3N_low"="#FF9D9A")
p <- ggplot(qc, aes(reorder(sampleID, group_code <- as.integer(factor(group))), mean_cov, fill = group)) +
  geom_col() + scale_fill_manual(values = grp_cols) +
  labs(x = NULL, y = "Mean CpG coverage (x)", fill = "Group") +
  coord_flip() + theme_bw()
ggsave("qc_coverage.png", p, width = 7, height = 6, dpi = 200)
```

# Filter and normalize

Remove CpGs below 5x and above the 99.9th coverage percentile, then
median-normalize coverage across samples.

```{r filter}
filt <- filterByCoverage(obj, lo.count = 5, lo.perc = NULL,
                         hi.count = NULL, hi.perc = 99.9)
filt <- normalizeCoverage(filt, method = "median")

retained <- data.table(
  sampleID = sapply(filt, function(x) x@sample.id),
  nCpG_filt = sapply(filt, function(x) nrow(getData(x))))
retained <- merge(meta[, .(sampleID, nCpG_raw = NA)], retained, by = "sampleID", all.y = TRUE)
fwrite(retained, "coverage_filter_retained.csv")

saveRDS(filt, "methylRaw_filtered.rds")   # checkpoint (large)
```

Expect 54-74% of CpGs retained per sample (5.2-8.0 M).