13.8 - GO enrichment of DMR genes
2026-10-10
Are genes near population differentially methylated regions (DMRs)
enriched for particular functions? DMRs come from 13.4 (1
kb tiles, SNP-masked); GO terms come from the Swiss-Prot best hit of
each GN gene (13.5).
Gene sets.
- DMR genes: genes whose body or
promoter_bpupstream overlaps a DMR, the same promoter definition as the13.4feature annotation. - Background: genes whose body or promoter overlaps a tile
that 13.4 actually tested (
tiles_tested.bed), so the comparison is against genes that could have produced a DMR, not the whole genome. Both sets are restricted to genes with at least one GO term.
Test. One-sided Fisher’s exact test per GO term,
with each gene’s terms propagated to all ancestor terms (the topGO
“classic” method, written out here because topGO is not in the lab
container). Terms with fewer than min_term or more than
max_term background genes are skipped; BH correction is per
ontology.
BH runs over every term in the size range; only terms with at least
min_hits DMR genes are then reported.
Caveats. Most DMRs are intergenic, so the DMR gene
sets are small and the test has little power. GO terms come from each
gene’s best Swiss-Prot hit, so tandem paralogs that share a hit (and
often sit in the same DMR) count as independent genes and can make one
locus look like an enriched term; check the genes column.
Treat individual terms as leads, not results.
Restricting to clean DMRs. Render with
params = list(dmr_class = "no genetic signal") to keep only
DMRs that the 13.9 bisulfite SNP check found free of
population-differentiated variants, CpG-destroying SNPs and repeat
signal. Other classes (untested, repeat-like,
genetic: CpG-destroying, genetic: other) work
the same way. Results then go to 13.8-dmr-go.<class>/
so the all-DMR run is kept. The background is unchanged (all genes on
tested tiles); since a DMR can only be classed clean where every oyster
has enough depth to test for SNPs, the clean set leans toward
well-covered regions, which the background does not.
Upstream: 13.4 (with
tiles_tested.bed) and 13.5; 13.9
when dmr_class is set. Paths are relative to
code/, like the other notebooks; runs in the lab R
container.
library(GenomicRanges)
library(rtracklayer)
library(GO.db)
library(tidyverse)
# tidyverse masks the range versions of reduce()/setdiff(); none are used here.dm_dir <- "../output/13-wgbs/13.4-methylation"
annot_file <- "../output/13-wgbs/13.5-gene-annotation/GN_gene_annotation.tsv"
genome_dir <- "../data/13-wgbs-genome"
snp_dir <- "../output/13-wgbs/13.9-bs-snp-check" # batch_DMR_<contrast>.tsv
dmr_class <- params$dmr_class # "all" = no 13.9 filter
out_dir <- paste0("../output/13-wgbs/13.8-dmr-go",
if (dmr_class == "all") "" else
paste0(".", gsub("[^A-Za-z0-9]+", "-", dmr_class)))
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
promoter_bp <- 2000
min_term <- 5 # background genes per term
max_term <- 500
min_hits <- 2 # DMR genes per term to report
fdr <- 0.05Genes, DMRs, tested tiles
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, "GN.gene.gff3"))
seqlevels(gff, pruning.mode = "coarse") <- seqlevels(si)
seqinfo(gff) <- si
genes <- gff[gff$type == "gene"]
names(genes) <- genes$ID
# Gene body plus upstream promoter, strand-aware.
gene_ext <- trim(resize(genes, width(genes) + promoter_bp, fix = "end"))
read_dm <- function(f) {
read_tsv(f, show_col_types = FALSE) %>%
{ GRanges(.$seqnames, IRanges(.$start, .$end), meth.diff = .$meth.diff,
feature = .$feature) }
}
dmr_files <- list.files(dm_dir, pattern = "^DMR_.*_vs_.*\\.tsv$", full.names = TRUE)
dmrs <- set_names(map(dmr_files, read_dm),
sub("^DMR_(.*)\\.tsv$", "\\1", basename(dmr_files)))
tiles <- import(file.path(dm_dir, "tiles_tested.bed"))
# 13.9 class per DMR (matched on coordinates), and the optional class filter.
dmrs <- imap(dmrs, function(d, cn) {
f <- file.path(snp_dir, paste0("batch_DMR_", cn, ".tsv"))
if (!file.exists(f)) {
if (dmr_class != "all") stop("dmr_class is set but ", f, " is missing; run 13.9 first")
d$snp_class <- NA_character_
return(d)
}
cls <- read_tsv(f, show_col_types = FALSE) %>%
transmute(key = paste(seqnames, start, end), class)
d$snp_class <- cls$class[match(paste(seqnames(d), start(d), end(d)), cls$key)]
stopifnot(!anyNA(d$snp_class)) # 13.9 must have been run on these DMRs
if (dmr_class == "all") d else d[d$snp_class == dmr_class]
})
message("DMR class filter: ", dmr_class)
map_int(dmrs, length)## Qui_vs_Rio
## 82length(tiles)## [1] 446760GO annotation
go_ids holds the Swiss-Prot best hit’s GO terms. Each
gene gets those terms plus all their ancestors in the same ontology.
annot <- read_tsv(annot_file, show_col_types = FALSE)
gene_go <- annot %>%
select(gene_id, go_ids) %>%
filter(!is.na(go_ids)) %>%
separate_rows(go_ids, sep = ";\\s*") %>%
rename(go_id = go_ids) %>%
filter(str_detect(go_id, "^GO:\\d{7}$")) %>%
distinct()
go_ont <- suppressMessages(AnnotationDbi::select(
GO.db, keys = unique(gene_go$go_id), columns = c("ONTOLOGY", "TERM"), keytype = "GOID"))
ancestors <- list(BP = as.list(GOBPANCESTOR), MF = as.list(GOMFANCESTOR),
CC = as.list(GOCCANCESTOR))
gene_go_prop <- gene_go %>%
inner_join(go_ont, by = c("go_id" = "GOID")) %>%
mutate(anc = map2(go_id, ONTOLOGY, ~ c(.x, setdiff(ancestors[[.y]][[.x]], "all")))) %>%
select(gene_id, ontology = ONTOLOGY, anc) %>%
unnest(anc) %>%
rename(go_id = anc) %>%
distinct()
term_names <- suppressMessages(AnnotationDbi::select(
GO.db, keys = unique(gene_go_prop$go_id), columns = "TERM", keytype = "GOID")) %>%
rename(go_id = GOID, term = TERM)
tibble(genes_with_go = n_distinct(gene_go$gene_id),
direct_terms = n_distinct(gene_go$go_id),
terms_after_propagation = n_distinct(gene_go_prop$go_id))## # A tibble: 1 x 3
## genes_with_go direct_terms terms_after_propagation
## <int> <int> <int>
## 1 9606 12401 17129Enrichment
genes_hit <- function(gr) names(genes)[unique(subjectHits(findOverlaps(gr, gene_ext, ignore.strand = TRUE)))]
universe <- intersect(genes_hit(tiles), unique(gene_go_prop$gene_id))
go_fisher <- function(set, universe) {
set <- intersect(set, universe)
gp <- filter(gene_go_prop, gene_id %in% universe)
gp %>%
group_by(ontology, go_id) %>%
summarise(bg = n(), hits = sum(gene_id %in% set),
genes = paste(sort(gene_id[gene_id %in% set]), collapse = ";"),
.groups = "drop") %>%
filter(bg >= min_term, bg <= max_term) %>%
mutate(n_set = length(set), n_universe = length(universe),
expected = bg * n_set / n_universe,
fold = hits / expected,
p = map2_dbl(hits, bg, ~ fisher.test(
matrix(c(.x, length(set) - .x, .y - .x,
length(universe) - length(set) - .y + .x), 2),
alternative = "greater")$p.value)) %>%
group_by(ontology) %>%
mutate(fdr = p.adjust(p, "BH")) %>% # over every term in the size range,
ungroup() %>% # before dropping low-hit terms
filter(hits >= min_hits) %>%
left_join(term_names, by = "go_id") %>%
arrange(ontology, p) %>%
select(ontology, go_id, term, hits, bg, expected, fold, p, fdr, genes)
}
results <- list()
gene_sets <- list()
for (cn in names(dmrs)) {
d <- dmrs[[cn]]
sets <- list(all = genes_hit(d),
hyper = genes_hit(d[d$meth.diff > 0]), # higher in the 2nd population
hypo = genes_hit(d[d$meth.diff < 0]))
for (dir in names(sets)) {
key <- paste(cn, dir, sep = "_")
gene_sets[[key]] <- sets[[dir]]
results[[key]] <- go_fisher(sets[[dir]], universe)
write_tsv(results[[key]], file.path(out_dir, paste0("GO_", key, ".tsv")))
}
}
summary_tbl <- imap_dfr(gene_sets, ~ tibble(
set = .y,
dmr_genes = length(.x),
dmr_genes_with_go = length(intersect(.x, universe)),
universe = length(universe),
terms_reported = nrow(results[[.y]]),
terms_fdr = sum(results[[.y]]$fdr < fdr)))
write_tsv(summary_tbl, file.path(out_dir, "GO_summary.tsv"))
summary_tbl## # A tibble: 3 x 6
## set dmr_genes dmr_genes_with_go universe terms_reported terms_fdr
## <chr> <int> <int> <int> <int> <int>
## 1 Qui_vs_Rio_all 8 7 6219 2 0
## 2 Qui_vs_Rio_hyper 4 4 6219 0 0
## 3 Qui_vs_Rio_hypo 4 3 6219 0 0Top terms
top <- imap_dfr(results, ~ .x %>% group_by(ontology) %>% slice_min(p, n = 10, with_ties = FALSE) %>%
ungroup() %>% mutate(set = .y)) %>%
filter(str_ends(set, "_all"))
top %>% select(set, ontology, term, hits, bg, fold, p, fdr)## # A tibble: 2 x 8
## set ontology term hits bg fold p fdr
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 Qui_vs_Rio_all BP nitrogen compound tran~ 2 487 3.65 0.0988 1
## 2 Qui_vs_Rio_all CC supramolecular complex 2 336 5.29 0.0510 1if (nrow(top) > 0) {
top %>%
mutate(term = fct_reorder(str_trunc(term, 60), -log10(p))) %>%
ggplot(aes(-log10(p), term, size = hits, colour = fdr < fdr)) +
geom_point() +
facet_grid(ontology ~ ., scales = "free_y", space = "free_y") +
scale_colour_manual(values = c(`FALSE` = "grey50", `TRUE` = "firebrick"),
name = paste("FDR <", fdr)) +
labs(x = "-log10 Fisher p", y = NULL,
title = "GO terms among DMR genes (body + promoter)") +
theme_bw()
}
DMR genes
Every gene near a DMR with its best Swiss-Prot hit, for reading through directly.
dmr_gene_tbl <- imap_dfr(dmrs, function(d, cn) {
ov <- findOverlaps(d, gene_ext, ignore.strand = TRUE)
tibble(contrast = cn, gene_id = names(genes)[subjectHits(ov)],
dmr = paste0(seqnames(d), ":", start(d), "-", end(d))[queryHits(ov)],
meth_diff = d$meth.diff[queryHits(ov)],
dmr_feature = as.character(d$feature[queryHits(ov)]),
snp_class = d$snp_class[queryHits(ov)])
}) %>%
left_join(select(annot, gene_id, entry_name, protein_name, evalue), by = "gene_id") %>%
arrange(contrast, desc(abs(meth_diff)))
write_tsv(dmr_gene_tbl, file.path(out_dir, "DMR_genes.tsv"))
dmr_gene_tbl %>% select(gene_id, meth_diff, dmr_feature, snp_class, protein_name)## # A tibble: 9 x 5
## gene_id meth_diff dmr_feature snp_class protein_name
## <chr> <dbl> <chr> <chr> <chr>
## 1 GN020540 -31.5 intron no genetic signal Nuclear pore complex protein~
## 2 GN017845 31.4 intron no genetic signal Protein farnesyltransferase ~
## 3 GN020540 -30.7 intron no genetic signal Nuclear pore complex protein~
## 4 GN011877 -28.3 promoter no genetic signal <NA>
## 5 GN020180 -26.1 promoter no genetic signal Testis-expressed protein 9
## 6 GN018537 -24.5 exon no genetic signal Probable methylmalonate-semi~
## 7 GN007260 22.2 exon no genetic signal Neutrophil cytosol factor 2 ~
## 8 GN020788 20.4 promoter no genetic signal Monocarboxylate transporter ~
## 9 GN012742 20.2 intron no genetic signal 3'-5' exoribonuclease HELZ2 ~sessionInfo()## R version 4.4.3 (2025-02-28)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
## [4] dplyr_1.1.4 purrr_1.0.4 readr_2.1.5
## [7] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.2
## [10] tidyverse_2.0.0 GO.db_3.20.0 AnnotationDbi_1.68.0
## [13] Biobase_2.66.0 rtracklayer_1.66.0 GenomicRanges_1.58.0
## [16] GenomeInfoDb_1.42.3 IRanges_2.40.1 S4Vectors_0.44.0
## [19] BiocGenerics_0.52.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 farver_2.1.2
## [3] blob_1.2.4 Biostrings_2.74.1
## [5] bitops_1.0-9 fastmap_1.2.0
## [7] RCurl_1.98-1.19 GenomicAlignments_1.42.0
## [9] XML_3.99-0.18 digest_0.6.37
## [11] timechange_0.3.0 lifecycle_1.0.4
## [13] KEGGREST_1.46.0 RSQLite_2.3.9
## [15] magrittr_2.0.3 compiler_4.4.3
## [17] rlang_1.1.5 sass_0.4.9
## [19] tools_4.4.3 utf8_1.2.4
## [21] yaml_2.3.10 knitr_1.50
## [23] labeling_0.4.3 S4Arrays_1.6.0
## [25] bit_4.6.0 curl_6.2.2
## [27] DelayedArray_0.32.0 abind_1.4-8
## [29] BiocParallel_1.40.2 withr_3.0.2
## [31] grid_4.4.3 colorspace_2.1-1
## [33] scales_1.3.0 SummarizedExperiment_1.36.0
## [35] cli_3.6.4 rmarkdown_2.29
## [37] crayon_1.5.3 generics_0.1.3
## [39] httr_1.4.7 tzdb_0.5.0
## [41] rjson_0.2.23 DBI_1.2.3
## [43] cachem_1.1.0 zlibbioc_1.52.0
## [45] parallel_4.4.3 XVector_0.46.0
## [47] restfulr_0.0.17 matrixStats_1.5.0
## [49] vctrs_0.6.5 Matrix_1.7-3
## [51] jsonlite_2.0.0 hms_1.1.3
## [53] bit64_4.6.0-1 jquerylib_0.1.4
## [55] glue_1.8.0 codetools_0.2-20
## [57] stringi_1.8.7 gtable_0.3.6
## [59] BiocIO_1.16.0 UCSC.utils_1.2.0
## [61] munsell_0.5.1 pillar_1.10.2
## [63] htmltools_0.5.8.1 GenomeInfoDbData_1.2.13
## [65] R6_2.6.1 vroom_1.6.5
## [67] evaluate_1.0.3 lattice_0.22-7
## [69] png_0.1-8 Rsamtools_2.22.0
## [71] memoise_2.0.1 bslib_0.9.0
## [73] SparseArray_1.6.2 xfun_0.52
## [75] MatrixGenerics_1.18.1 pkgconfig_2.0.3