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.
Upstream: 13.4 (with
tiles_tested.bed) and 13.5. 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"
out_dir <- "../output/13-wgbs/13.8-dmr-go"
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"))
map_int(dmrs, length)## Qui_vs_Rio
## 304length(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 47 38 6219 203 0
## 2 Qui_vs_Rio_hyper 32 26 6219 143 0
## 3 Qui_vs_Rio_hypo 15 12 6219 18 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: 30 x 8
## set ontology term hits bg fold p fdr
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 Qui_vs_Rio_all BP lysosome localization 2 15 21.8 0.00363 1
## 2 Qui_vs_Rio_all BP vacuolar localization 2 15 21.8 0.00363 1
## 3 Qui_vs_Rio_all BP filopodium assembly 2 20 16.4 0.00645 1
## 4 Qui_vs_Rio_all BP positive regulation ~ 2 22 14.9 0.00778 1
## 5 Qui_vs_Rio_all BP alpha-amino acid met~ 3 80 6.14 0.0125 1
## 6 Qui_vs_Rio_all BP regulation of plasma~ 2 42 7.79 0.0268 1
## 7 Qui_vs_Rio_all BP regulation of cell p~ 2 43 7.61 0.0281 1
## 8 Qui_vs_Rio_all BP gonad development 2 50 6.55 0.0371 1
## 9 Qui_vs_Rio_all BP development of prima~ 2 50 6.55 0.0371 1
## 10 Qui_vs_Rio_all BP regulation of microt~ 2 50 6.55 0.0371 1
## # i 20 more rowsif (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)]))
}) %>%
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, protein_name)## # A tibble: 56 x 4
## gene_id meth_diff dmr_feature protein_name
## <chr> <dbl> <chr> <chr>
## 1 GN003827 -86.0 promoter BTB/POZ domain-containing protein KCTD7
## 2 GN020540 -72.0 exon Nuclear pore complex protein Nup107 (107 kDa ~
## 3 GN015239 61.3 exon <NA>
## 4 GN003810 59.5 exon Cytoplasmic tRNA 2-thiolation protein 2
## 5 GN020540 -57.0 promoter Nuclear pore complex protein Nup107 (107 kDa ~
## 6 GN000740 -55.5 intron <NA>
## 7 GN010445 51.0 exon <NA>
## 8 GN021273 47.1 intron <NA>
## 9 GN008799 44.5 promoter <NA>
## 10 GN018370 -44.2 promoter Pre-B-cell leukemia transcription factor 4 (H~
## # i 46 more rowssessionInfo()## 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