--- title: "13.4 - WGBS DNA methylation across three populations" output: html_document: toc: true toc_float: true date: "2026-10-08" params: ref: "merged" # or "hapa": the WGBS_REF=hapa alignment, with 13.11 lifted genes --- Population differences in CpG methylation in *Ostrea chilensis*, from whole-genome bisulfite sequencing of oysters from Pumalin, Quilhua, and Rio Pudeto. **Sample set.** 15 libraries were delivered (5 per site), but several FASTQ pairs share reads with others (`output/13-wgbs/raw_read_overlap.tsv`; see README Known issues). After removing duplicates, `data/13-wgbs-samples.tsv` holds 9 independent oysters: Pum1, Qui1-4, Rio1-4. Only populations with at least `min_pop_n` oysters are tested, so Pumalin (n = 1) is left out of this notebook until its libraries are re-delivered; the contrasts below are built from whichever populations remain. **Upstream:** `code/13-wgbs-submit.sh` (13.1 genome prep → 13.2 trim/align/extract → 13.3 QC summary) must be finished. This notebook reads the strand-merged CpG files `output/13-wgbs/methylation/*.CpG_report.merged_CpG_evidence.cov.gz`. **Reference:** `merged_out.fasta` + `GN.gene.gff3`, the same as the current RNA-seq pipeline (`09`), so differentially methylated genes and DE genes share IDs. Render with `params = list(ref = "hapa")` to run on the `WGBS_REF=hapa` alignment to `Och_HapA_assembly.fa` instead: it reads `output/13-wgbs-hapa/` (methylation calls, alignment summary, the HapA `13.7` SNP mask), annotates with the GN genes lifted onto HapA by `13.11` (same gene IDs), and writes `output/13-wgbs-hapa/13.4-methylation`. DMR/DML coordinates are then on HapA, not the merged reference. **Where to run:** on a compute node inside the lab R container, e.g. ``` salloc -A srlab -p cpu-g2-mem2x -c 16 --mem=150G -t 12:00:00 apptainer exec /gscratch/srlab/containers/srlab-R4.4-bioinformatics-container-c3d3116.sif R ``` Paths are relative to `code/`, like the other notebooks. **Depth caveat.** At ~23-38 M read pairs per sample, raw coverage of the 1.18 Gb genome is only ~6-9x before mapping and deduplication. Single-CpG tests will therefore be underpowered and cover a fraction of sites. The analysis is layered so that the region (tile) and gene-body levels, which pool CpGs, carry most of the inference; treat single-CpG results as supporting. ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.width = 8, fig.height = 5) ``` ```{r packages} library(methylKit) library(GenomicRanges) library(rtracklayer) library(limma) library(pheatmap) library(tidyverse) # tidyverse is loaded last, so it masks methylKit::unite (tidyr), and the range # versions of reduce (purrr) and setdiff (dplyr/lubridate); those calls are namespaced. ``` ```{r parameters} ref <- match.arg(params$ref, c("merged", "hapa")) ref_cfg <- list( merged = list(out_root = "../output/13-wgbs", genome_dir = "../data/13-wgbs-genome", gff = "GN.gene.gff3", assembly = "Och_merged_out"), hapa = list(out_root = "../output/13-wgbs-hapa", genome_dir = "../data/13-wgbs-genome-hapa", gff = "GN.gene.hapa.gff3", assembly = "Och_HapA") # gff from 13.11 )[[ref]] out_root <- ref_cfg$out_root out_dir <- file.path(out_root, "13.4-methylation") meth_dir <- file.path(out_root, "methylation") genome_dir <- ref_cfg$genome_dir dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) min_cov <- 5 # per-CpG minimum (strand-merged); 10 is typical but costs most sites here max_cov_pctl <- 99.9 # drop PCR/repeat pile-ups snp_mask <- file.path(out_root, "13.7-snps/wgbs_snps.q20.CpG.bed") # from 13.7; skipped if absent min_pop_n <- 3L # populations with fewer oysters are not tested min_per_group <- 3L # CpG must pass in >= 3 oysters (of 4) in every population tile_size <- 1000 # bp, non-overlapping windows for DMRs tile_min_cpgs <- 5 # CpGs required per tile per sample min_gene_cpgs <- 5 # CpGs required per gene for gene-body methylation diff_cutoff <- 20 # % methylation difference for DML/DMR calls q_cutoff <- 0.01 cores <- 8 ``` # Samples Groups come from the explicit table `data/13-wgbs-samples.tsv`, never from file order, so a re-sorted directory can't silently swap populations. ```{r samples} samples_all <- read_tsv("../data/13-wgbs-samples.tsv", show_col_types = FALSE) pop_n <- count(samples_all, population) untested <- filter(pop_n, n < min_pop_n) if (nrow(untested) > 0) message("Not tested (fewer than ", min_pop_n, " oysters): ", paste0(untested$population, " (n = ", untested$n, ")", collapse = ", ")) samples <- samples_all %>% filter(!population %in% untested$population) %>% mutate(population = factor(population, levels = intersect(c("Pum", "Qui", "Rio"), population)), treatment = as.integer(population) - 1L, # 0, 1, ... for methylKit file = file.path(meth_dir, paste0(sample, ".CpG_report.merged_CpG_evidence.cov.gz"))) stopifnot(all(file.exists(samples$file)), nlevels(samples$population) >= 2) samples %>% select(sample, population, site) ``` # Alignment QC From 13.3. CHH methylation is the conversion check: invertebrate non-CpG methylation is near zero, so CHH % ≈ 100 − conversion efficiency. Flag any sample above ~1%. ```{r align-qc} aln <- read_tsv(file.path(out_root, "alignment_summary.tsv"), show_col_types = FALSE) %>% inner_join(select(samples_all, sample, population), by = "sample") # clean samples only aln aln %>% pivot_longer(c(mapping_pct, dup_pct, CpG_meth_pct, CHH_meth_pct)) %>% ggplot(aes(sample, value, fill = population)) + geom_col() + facet_wrap(~name, scales = "free_y") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) ``` # Read CpG calls `bismarkCoverage` format: chr, start, end, % methylation, count methylated, count unmethylated (1-based). Strands are already merged, so no `destrand` later. ```{r read} meth_raw <- methRead( location = as.list(samples$file), sample.id = as.list(samples$sample), treatment = samples$treatment, assembly = ref_cfg$assembly, pipeline = "bismarkCoverage", context = "CpG", mincov = 1 ) ``` ## Mask CpG SNPs A C->T SNP at a CpG (G->A on the other strand) reads as an unmethylated CpG, so a SNP that differs in frequency between populations can look like a methylation difference. CpGs where either base carries a SNP called from these same bisulfite BAMs (`13.7-wgbs-snps.sh`, QUAL >= 20) are removed from every sample before any filtering or testing. Strand-merged calls sit at the C of the CpG, so each site spans `start` to `start + 1`. ```{r snp-mask} if (file.exists(snp_mask)) { mask <- import(snp_mask, format = "bed") n_before <- map_int(meth_raw, ~ nrow(getData(.x))) meth_raw <- methylRawList( map(meth_raw, function(x) { d <- getData(x) hit <- overlapsAny(GRanges(d$chr, IRanges(d$start, d$start + 1)), mask) methylKit::select(x, which(!hit)) }), treatment = samples$treatment) n_after <- map_int(meth_raw, ~ nrow(getData(.x))) tibble(sample = samples$sample, cpgs = n_before, masked = n_before - n_after, pct_masked = round(100 * (n_before - n_after) / n_before, 3)) } else { warning("No SNP mask at ", snp_mask, " - CpG SNPs are NOT removed. Run 13.7 first.") } ``` ```{r coverage-stats, fig.height=8} # Coverage distribution per sample, before filtering. cov_summary <- map_dfr(seq_along(meth_raw), function(i) { d <- getData(meth_raw[[i]]) tibble(sample = samples$sample[i], cpgs_1x = nrow(d), cpgs_5x = sum(d$coverage >= 5), cpgs_10x = sum(d$coverage >= 10), median_cov = median(d$coverage), # coverage-weighted global CpG methylation global_meth_pct = 100 * sum(d$numCs) / sum(d$coverage)) }) write_tsv(cov_summary, file.path(out_dir, "cpg_coverage_summary.tsv")) cov_summary ggplot(left_join(cov_summary, samples, by = "sample"), aes(population, global_meth_pct, colour = population)) + geom_jitter(width = 0.1, size = 3) + labs(y = "Global CpG methylation (%)", title = "Global methylation by population") + theme_bw() ``` # Filter, normalize, unite ```{r filter-unite} meth_filt <- filterByCoverage(meth_raw, lo.count = min_cov, hi.perc = max_cov_pctl) meth_norm <- normalizeCoverage(meth_filt, method = "median") meth_unite <- methylKit::unite(meth_norm, destrand = FALSE, min.per.group = min_per_group, mc.cores = cores) message(nrow(meth_unite), " CpGs pass in >= ", min_per_group, " oysters per population") pm <- percMethylation(meth_unite) # CpG x sample, NAs where a sample was filtered ``` ## Distribution Invertebrate methylomes are typically mosaic/bimodal: most CpGs near 0% or near 100%, with methylation concentrated in gene bodies. ```{r distribution} as_tibble(pm) %>% pivot_longer(everything(), names_to = "sample", values_to = "pct") %>% filter(!is.na(pct)) %>% left_join(samples, by = "sample") %>% ggplot(aes(pct, colour = population, group = sample)) + geom_density() + labs(x = "% methylation per CpG") + theme_bw() ``` ## Sample structure PCA on CpGs covered in every tested oyster, using the most variable sites. ```{r pca} pm_complete <- pm[complete.cases(pm), ] top <- head(order(apply(pm_complete, 1, var), decreasing = TRUE), 50000) pca <- prcomp(t(pm_complete[top, ])) ve <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1) as_tibble(pca$x[, 1:3], rownames = "sample") %>% left_join(samples, by = "sample") %>% ggplot(aes(PC1, PC2, colour = population, label = sample)) + geom_point(size = 3) + geom_text(vjust = -0.8, size = 3, show.legend = FALSE) + labs(x = paste0("PC1 (", ve[1], "%)"), y = paste0("PC2 (", ve[2], "%)"), subtitle = paste(nrow(pm_complete), "complete CpGs; top", length(top), "by variance")) + theme_bw() pheatmap(cor(pm_complete[top, ], method = "spearman"), annotation_col = data.frame(population = samples$population, row.names = samples$sample), main = "Spearman correlation, top variable CpGs") ``` # Genomic features Annotation is built directly from `GN.gene.gff3`. Each CpG is assigned to one feature by priority exon > intron > promoter (2 kb upstream of a gene) > intergenic. ```{r features} seqlens <- read_tsv(file.path(genome_dir, "genome.chrom.sizes"), col_names = c("chr", "len"), show_col_types = FALSE) si <- Seqinfo(seqlens$chr, seqlens$len) gff <- import(file.path(genome_dir, ref_cfg$gff)) seqlevels(gff, pruning.mode = "coarse") <- seqlevels(si) # align to FASTA order seqinfo(gff) <- si # lengths needed for trim()/gaps() genes <- gff[gff$type == "gene"] names(genes) <- genes$ID exons <- IRanges::reduce(gff[gff$type == "exon"], ignore.strand = TRUE) introns <- BiocGenerics::setdiff(IRanges::reduce(genes, ignore.strand = TRUE), exons) upstream <- trim(promoters(genes, upstream = 2000, downstream = 0)) promoters <- BiocGenerics::setdiff(IRanges::reduce(upstream, ignore.strand = TRUE), IRanges::reduce(genes, ignore.strand = TRUE)) intergenic <- gaps(IRanges::reduce(c(granges(genes), granges(promoters)), ignore.strand = TRUE)) intergenic <- intergenic[strand(intergenic) == "*"] feature_sets <- list(exon = exons, intron = introns, promoter = promoters, intergenic = intergenic) annotate_feature <- function(gr) { out <- rep(NA_character_, length(gr)) for (f in names(feature_sets)) { hit <- is.na(out) & overlapsAny(gr, feature_sets[[f]], ignore.strand = TRUE) out[hit] <- f } factor(out, levels = names(feature_sets)) } cpg_gr <- as(meth_unite, "GRanges") cpg_gr$feature <- annotate_feature(cpg_gr) cpg_gr$mean_meth <- rowMeans(pm, na.rm = TRUE) as_tibble(mcols(cpg_gr)) %>% ggplot(aes(feature, mean_meth)) + geom_violin(fill = "grey80") + stat_summary(fun = mean, geom = "point", colour = "red") + labs(y = "Mean % methylation across samples") + theme_bw() table(cpg_gr$feature) ``` # Differential methylation Every pairwise contrast among the tested populations, named as in the RNA-seq DE (`Pum_vs_Qui`, `Pum_vs_Rio`, `Qui_vs_Rio`). Positive `meth.diff` = higher in the second population. ```{r contrasts} pairs <- combn(levels(samples$population), 2, simplify = FALSE) contrasts <- set_names(pairs, map_chr(pairs, paste, collapse = "_vs_")) names(contrasts) # methylKit's binomial test with overdispersion correction, on a 2-population subset. run_pair <- function(obj, groups) { keep <- samples %>% filter(population %in% groups) sub <- reorganize(obj, sample.ids = keep$sample, treatment = as.integer(keep$population == groups[2])) calculateDiffMeth(sub, overdispersion = "MN", test = "Chisq", mc.cores = cores) } write_hits <- function(diff, name, level) { hits <- getMethylDiff(diff, difference = diff_cutoff, qvalue = q_cutoff, type = "all") gr <- as(hits, "GRanges") gr$feature <- annotate_feature(gr) write_tsv(as_tibble(as.data.frame(gr)), file.path(out_dir, paste0(level, "_", name, ".tsv"))) export.bed(gr, file.path(out_dir, paste0(level, "_", name, ".bed"))) gr } ``` ## Single CpGs (DMLs) ```{r dml} dml <- imap(contrasts, function(g, name) { d <- run_pair(meth_unite, g) write_hits(d, name, "DML") }) map_int(dml, length) ``` ## 1 kb tiles (DMRs) Tiling pools reads across neighbouring CpGs, which is the better-powered test at this depth. ```{r dmr} tiles <- tileMethylCounts(meth_norm, win.size = tile_size, step.size = tile_size, cov.bases = tile_min_cpgs, mc.cores = cores) tiles_unite <- methylKit::unite(tiles, destrand = FALSE, min.per.group = min_per_group, mc.cores = cores) message(nrow(tiles_unite), " tiles tested") # Tested tiles = the background for DMR gene-set tests (13.8). export.bed(as(tiles_unite, "GRanges"), file.path(out_dir, "tiles_tested.bed")) dmr <- imap(contrasts, function(g, name) { d <- run_pair(tiles_unite, g) write_hits(d, name, "DMR") }) dmr_summary <- imap_dfr(dmr, ~ tibble(contrast = .y, n_dmr = length(.x), hyper = sum(.x$meth.diff > 0), hypo = sum(.x$meth.diff < 0))) %>% left_join(imap_dfr(dmr, ~ as_tibble(table(feature = .x$feature)) %>% mutate(contrast = .y)) %>% pivot_wider(names_from = feature, values_from = n), by = "contrast") write_tsv(dmr_summary, file.path(out_dir, "DMR_summary.tsv")) dmr_summary ``` # Gene-body methylation Per gene and oyster: methylated / total calls summed over all CpGs in the gene (a coverage-weighted mean). Genes need at least `min_gene_cpgs` CpGs with data in every sample. Tested with limma on M-values (logit), the same framework and contrasts as the RNA-seq DE in `09`, so results can be joined on gene ID. ```{r gene-body} cov_mat <- as.matrix(getData(meth_unite)[, meth_unite@coverage.index]) cs_mat <- as.matrix(getData(meth_unite)[, meth_unite@numCs.index]) colnames(cov_mat) <- colnames(cs_mat) <- meth_unite@sample.ids cov_mat[is.na(cov_mat)] <- 0 cs_mat[is.na(cs_mat)] <- 0 ov <- findOverlaps(cpg_gr, genes, ignore.strand = TRUE) gene_id <- names(genes)[subjectHits(ov)] sum_by_gene <- function(m) rowsum(m[queryHits(ov), , drop = FALSE], gene_id) gene_cov <- sum_by_gene(cov_mat) gene_cs <- sum_by_gene(cs_mat) gene_ncpg <- sum_by_gene((cov_mat > 0) * 1) keep <- apply(gene_ncpg, 1, min) >= min_gene_cpgs gene_frac <- (gene_cs[keep, ] + 0.5) / (gene_cov[keep, ] + 1) # offset avoids logit(0/1) gene_M <- log2(gene_frac / (1 - gene_frac)) message(sum(keep), " genes with >= ", min_gene_cpgs, " covered CpGs in every oyster") write_tsv(as_tibble(100 * gene_frac, rownames = "gene_id"), file.path(out_dir, "gene_body_methylation_pct.tsv")) # Gene-body methylation is typically bimodal in bivalves (methylated vs unmethylated genes). tibble(mean_pct = rowMeans(100 * gene_frac)) %>% ggplot(aes(mean_pct)) + geom_histogram(bins = 50) + labs(x = "Mean gene-body methylation (%)", y = "Genes") + theme_bw() ``` ```{r gene-body-limma} group <- samples$population[match(colnames(gene_M), samples$sample)] design <- model.matrix(~ 0 + group); colnames(design) <- levels(group) cm <- makeContrasts(contrasts = map_chr(contrasts, ~ paste(.x[2], "-", .x[1])), levels = design) colnames(cm) <- names(contrasts) fit <- eBayes(contrasts.fit(lmFit(gene_M, design), cm), trend = TRUE) # Omnibus: any difference among the tested populations (with two populations # this is just their single contrast). any_pop <- topTable(fit, coef = if (ncol(cm) == 1) 1 else NULL, number = Inf, sort.by = if (ncol(cm) == 1) "P" else "F") %>% as_tibble(rownames = "gene_id") write_tsv(any_pop, file.path(out_dir, "DMG_any_population_F.tsv")) dmg <- map(colnames(cm), function(cn) { tt <- topTable(fit, coef = cn, number = Inf) %>% as_tibble(rownames = "gene_id") # logFC here is a difference in M-values; also report it in % methylation. tt <- tt %>% mutate(delta_pct = 100 * (rowMeans(gene_frac[gene_id, group == sub(".*_vs_", "", cn)]) - rowMeans(gene_frac[gene_id, group == sub("_vs_.*", "", cn)]))) write_tsv(tt, file.path(out_dir, paste0("DMG_", cn, ".tsv"))) filter(tt, adj.P.Val < 0.05, abs(delta_pct) >= 10) }) %>% set_names(colnames(cm)) map_int(dmg, nrow) sig <- any_pop %>% filter(adj.P.Val < 0.05) %>% pull(gene_id) if (length(sig) >= 2) { pheatmap(100 * gene_frac[head(sig, 100), ], scale = "row", show_rownames = length(sig) <= 50, annotation_col = data.frame(population = group, row.names = colnames(gene_M)), main = "Top differentially methylated genes (F-test)") } ``` # Next steps - **Expression link:** join `DMG_*.tsv` to the `09` limma tables on `gene_id` (GN gene IDs; the `09` StringTie gene IDs map to these via the reference-guided merge) and test whether gene-body methylation level or change predicts expression. Note the RNA-seq oysters are different individuals from the WGBS oysters, so this is a population-level, not individual-level, comparison. - **Genetic confounding:** CpG SNPs from 13.7 are masked above, but at ~6-9x per oyster many SNPs (especially rare ones) go uncalled, so a call with no SNP is not proof of a clean site. Calls that are near 0% in one population and 100% in the other remain the ones to check. - **Functional enrichment:** GO enrichment of DMGs against the background of all tested genes (`any_pop$gene_id`), using `SPannot-GO.tsv` via the `03` join. ```{r session} sessionInfo() ```