13.6 - Genome-wide DNA methylation landscape
2026-10-08
General characterization of CpG methylation in Ostrea
chilensis: how much of the genome is methylated, how methylation is
distributed across genomic features and along chromosomes, what
gene-body methylation looks like, and whether gene methylation tracks
gene function. Population differences are in 13.4; this
notebook describes the shared methylome, pooling the 9 independent
oysters in data/13-wgbs-samples.tsv (duplicated libraries
excluded; see README Known issues).
Upstream:
code/13-wgbs-submit.sh(13.1 → 13.2 → 13.3) finished: reads the strand-merged CpG filesoutput/13-wgbs/methylation/*.CpG_report.merged_CpG_evidence.cov.gz.code/13.5-wgbs-gene-annotation.shfinished: readsoutput/13-wgbs/13.5-gene-annotation/GN_gene_annotation.tsv(Swiss-Prot best hit + GO per GN gene). Only the last section needs it.
Reference: merged_out.fasta +
GN.gene.gff3, as in 09 and 13.4.
The GFF has gene / mRNA / exon / CDS only, with every exon identical to
a CDS: there are no UTRs, so “exon” here means coding
exon and 5’/3’ UTR methylation can’t be separated from the flanks.
Pooling. At ~6–9x raw coverage per oyster, single-sample per-CpG estimates are noisy. For the landscape, counts are summed across all 9 oysters per CpG (a species-level methylome at roughly 9x the depth of one oyster). Per-oyster values are kept for global, per-feature, and per-gene summaries, which pool many CpGs and so are stable within a single animal; these show whether the patterns are consistent across individuals and populations.
Where to run: a compute node inside the lab R container (~100 GB RAM), e.g.
salloc -A srlab -p cpu-g2-mem2x -c 8 --mem=150G -t 8:00:00
apptainer exec /gscratch/srlab/containers/srlab-R4.4-bioinformatics-container-c3d3116.sif RPaths are relative to code/, like the other
notebooks.
# tidyverse first, so the Bioconductor set operations (reduce, setdiff, union, ...)
# mask the dplyr/purrr versions rather than the other way round.
library(tidyverse)
library(data.table)
library(GenomicRanges)
library(rtracklayer)
library(Biostrings)
library(Rsamtools)out_dir <- "../output/13-wgbs/13.6-landscape"
meth_dir <- "../output/13-wgbs/methylation"
genome_dir <- "../data/13-wgbs-genome"
annot_file <- "../output/13-wgbs/13.5-gene-annotation/GN_gene_annotation.tsv"
repeat_gff <- NULL # RepeatMasker/EDTA GFF for this assembly, if one becomes available
snp_bed <- "../output/13-wgbs/13.7-snps/wgbs_snps.q20.CpG.bed" # from 13.7; skipped if absent
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
# CpG level (pooled across oysters)
min_cov_pooled <- 10 # summed reads per CpG
max_cov_pctl <- 0.999 # drop pile-ups (collapsed repeats)
min_samples <- 5L # CpG must have data in >= 5 of 9 oysters
meth_breaks <- c(-Inf, 10, 50, Inf) # % methylation
meth_labels <- c("unmethylated (<10%)", "sparse (10-50%)", "methylated (>=50%)")
# features
flank_bp <- 2000 # promoter (upstream) and downstream flank; same promoter as 13.4
window_size <- 1e5 # chromosome-scale windows
# genes
min_gene_cpgs <- 5 # covered CpGs needed for a gene-body estimate
min_gene_calls <- 30 # per-oyster reads summed over the gene, for per-oyster values
n_body_bins <- 20 # metagene: scaled gene body
n_flank_bins <- 10 # metagene: each 2 kb flank
# function
min_go_genes <- 10
max_go_genes <- 500Samples
samples <- read_tsv("../data/13-wgbs-samples.tsv", show_col_types = FALSE) %>%
mutate(population = factor(population, levels = c("Pum", "Qui", "Rio")),
file = file.path(meth_dir,
paste0(sample, ".CpG_report.merged_CpG_evidence.cov.gz")))
stopifnot(all(file.exists(samples$file)))Genome and features
Five non-overlapping compartments, so every base and CpG belongs to exactly one:
| feature | definition |
|---|---|
| exon | coding exon (no UTRs annotated) |
| intron | gene body minus exons |
| promoter | 2 kb upstream of a TSS, outside any gene |
| downstream | 2 kb downstream of a TES, outside any gene or promoter |
| intergenic | everything else |
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)
fa <- FaFile(file.path(genome_dir, "merged_out.fasta"))
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
exon_gr <- gff[gff$type == "exon"]
stopifnot(all(lengths(exon_gr$Parent) == 1))
exon_gr$tx <- as.character(unlist(exon_gr$Parent))
exon_gr$gene_id <- sub("\\.[0-9]+$", "", exon_gr$tx)
gene_red <- reduce(genes, ignore.strand = TRUE)
exons <- reduce(exon_gr, ignore.strand = TRUE)
introns <- setdiff(gene_red, exons)
promoters <- setdiff(reduce(trim(promoters(genes, upstream = flank_bp, downstream = 0)),
ignore.strand = TRUE), gene_red)
downstream <- setdiff(reduce(trim(flank(genes, flank_bp, start = FALSE)), ignore.strand = TRUE),
union(gene_red, promoters))
intergenic <- gaps(reduce(c(gene_red, promoters, downstream)))
intergenic <- intergenic[strand(intergenic) == "*"]
feature_levels <- c("exon", "intron", "promoter", "downstream", "intergenic")
features <- c(exons, introns, promoters, downstream, intergenic)
features$feature <- factor(rep(feature_levels,
lengths(list(exons, introns, promoters, downstream, intergenic))),
levels = feature_levels)
stopifnot(!any(countOverlaps(features, features) > 1)) # disjoint
stopifnot(sum(as.numeric(width(features))) == sum(as.numeric(seqlengths(si))))
annotate_feature <- function(gr) {
features$feature[findOverlaps(gr, features, ignore.strand = TRUE, select = "first")]
}
# Every CpG dinucleotide in the assembly (1-based position of the C, as in Bismark).
genome_seq <- scanFa(fa)
names(genome_seq) <- sub(" .*", "", names(genome_seq))
genome_cpg <- do.call(c, lapply(seqlevels(si), function(chr) {
GRanges(chr, IRanges(start(matchPattern("CG", genome_seq[[chr]])), width = 1),
seqinfo = si)
}))
genome_cpg$feature <- annotate_feature(genome_cpg)
rm(genome_seq); invisible(gc())
message(length(genome_cpg), " CpGs in the assembly")# Observed/expected CpG per feature: CG * length / (C * G). Germline methylation
# deaminates mC -> T over evolutionary time, so long-term methylated sequence is CpG-depleted.
composition <- map_dfr(feature_levels, function(f) {
gr <- features[features$feature == f]
seq <- getSeq(fa, gr)
bf <- colSums(letterFrequency(seq, c("C", "G")))
tibble(feature = f,
bp = sum(as.numeric(width(gr))),
n_cpg = sum(as.numeric(vcountPattern("CG", seq))),
gc_pct = 100 * sum(bf) / bp,
cpg_oe = n_cpg * bp / (bf[["C"]] * bf[["G"]]))
}) %>%
mutate(pct_genome = 100 * bp / sum(bp),
cpg_per_kb = 1000 * n_cpg / bp)
composition## # A tibble: 5 x 7
## feature bp n_cpg gc_pct cpg_oe pct_genome cpg_per_kb
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 exon 25932309 724981 43.9 0.581 2.20 28.0
## 2 intron 109651400 1756427 34.4 0.542 9.29 16.0
## 3 promoter 39443476 892857 35.5 0.718 3.34 22.6
## 4 downstream 37404907 745204 35.1 0.647 3.17 19.9
## 5 intergenic 968488830 21980792 35.9 0.706 82.0 22.7Read CpG calls
Each oyster is read once. Per-oyster summaries (global, per feature, per gene) are taken on the way through; per-CpG counts are summed into one pooled table.
SNP mask. CpGs carrying a SNP on their C or G (from
13.7, biscuit genotyping of the bisulfite BAMs) are dropped
as each oyster is read, since a C→T SNP reads as an unmethylated CpG.
Every summary below is therefore SNP-masked; assembly CpG counts
(genome_cpg) are not. If the 13.7 file is missing, nothing
is masked.
read_cov <- function(f) {
fread(cmd = paste("gzip -dc", shQuote(f)), header = FALSE,
select = c(1, 2, 5, 6), col.names = c("chr", "pos", "meth", "unmeth"))
}
snp_gr <- NULL
if (file.exists(snp_bed)) {
snp_gr <- import(snp_bed, format = "bed") # one base: the C or G hit
message(length(snp_gr), " CpG-SNP sites from ", snp_bed)
} else {
message("No SNP mask (", snp_bed, " not found); run code/13.7-wgbs-snps.sh to add one")
}
pooled <- NULL
ps_global <- ps_feature <- ps_gene <- vector("list", nrow(samples))
n_masked <- set_names(integer(nrow(samples)), samples$sample)
for (i in seq_len(nrow(samples))) {
s <- samples$sample[i]
d <- read_cov(samples$file[i])[chr %in% seqlevels(si)]
if (!is.null(snp_gr)) { # each CpG spans its C (pos) and G (pos + 1)
hit <- overlapsAny(GRanges(d$chr, IRanges(d$pos, width = 2)), snp_gr, ignore.strand = TRUE)
n_masked[s] <- sum(hit)
d <- d[!hit]
}
d[, cov := meth + unmeth]
gr <- GRanges(d$chr, IRanges(d$pos, width = 1), seqinfo = si)
d[, feature := annotate_feature(gr)]
ps_global[[i]] <- d[, .(sample = s, cpgs_1x = .N, cpgs_5x = sum(cov >= 5),
median_cov = median(cov),
global_meth_pct = 100 * sum(meth) / sum(cov),
methylated_cpg_pct = 100 * mean(meth[cov >= 5] / cov[cov >= 5] >= 0.5))]
ps_feature[[i]] <- d[, .(sample = s, meth = sum(meth), cov = sum(cov), n_cpg = .N),
by = feature]
ov <- findOverlaps(gr, genes, ignore.strand = TRUE)
ps_gene[[i]] <- data.table(gene_id = names(genes)[subjectHits(ov)],
meth = d$meth[queryHits(ov)], cov = d$cov[queryHits(ov)]
)[, .(sample = s, meth = sum(meth), cov = sum(cov), n_cpg = .N),
by = gene_id]
pooled <- rbind(pooled, d[, .(chr, pos, meth, unmeth, n = 1L)]
)[, .(meth = sum(meth), unmeth = sum(unmeth), n = sum(n)), by = .(chr, pos)]
message(s, ": ", nrow(d), " CpGs; pooled table now ", nrow(pooled))
}
ps_global <- rbindlist(ps_global) %>% left_join(samples, by = "sample")
ps_feature <- rbindlist(ps_feature) %>% left_join(samples, by = "sample")
ps_gene <- rbindlist(ps_gene)
if (!is.null(snp_gr))
write_tsv(tibble(sample = names(n_masked), cpgs_masked = n_masked),
file.path(out_dir, "snp_mask_summary.tsv"))
rm(d, gr, ov); invisible(gc())pooled[, cov := meth + unmeth]
cov_max <- quantile(pooled$cov, max_cov_pctl)
pooled_f <- pooled[cov >= min_cov_pooled & cov <= cov_max & n >= min_samples]
cpg <- GRanges(pooled_f$chr, IRanges(pooled_f$pos, width = 1), seqinfo = si,
meth = pooled_f$meth, cov = pooled_f$cov,
pct = 100 * pooled_f$meth / pooled_f$cov)
cpg$feature <- annotate_feature(cpg)
cpg$class <- cut(cpg$pct, meth_breaks, meth_labels, right = FALSE)
tibble(step = c("CpGs in assembly", "covered in >= 1 oyster", "pass pooled filters"),
n = c(length(genome_cpg), nrow(pooled), length(cpg))) %>%
mutate(pct_of_assembly = round(100 * n / length(genome_cpg), 1))## # A tibble: 3 x 3
## step n pct_of_assembly
## <chr> <int> <dbl>
## 1 CpGs in assembly 26104181 100
## 2 covered in >= 1 oyster 18689605 71.6
## 3 pass pooled filters 14433110 55.3# Pooled per-CpG methylation track for IGV (alongside the 13.2 BAMs).
bw <- granges(cpg); bw$score <- cpg$pct
export.bw(sort(bw), file.path(out_dir, "pooled_CpG_methylation.bw"))
fwrite(pooled_f[, .(chr, start = pos, end = pos + 1L, pct = round(100 * meth / cov, 2),
meth, unmeth, n_samples = n)],
file.path(out_dir, "pooled_CpG_methylation.tsv.gz"), sep = "\t")
rm(pooled, pooled_f); invisible(gc())Global methylation
Per oyster
CpG methylation in invertebrates is moderate overall (bivalves typically ~10–20% of CpGs methylated), unlike the ~70–80% of vertebrates.
aln <- read_tsv("../output/13-wgbs/alignment_summary.tsv", show_col_types = FALSE)
global <- ps_global %>%
left_join(select(aln, sample, CHG_meth_pct, CHH_meth_pct), by = "sample") %>%
select(sample, population, cpgs_1x, cpgs_5x, median_cov,
global_meth_pct, methylated_cpg_pct, CHG_meth_pct, CHH_meth_pct)
write_tsv(global, file.path(out_dir, "global_methylation_per_oyster.tsv"))
global## sample population cpgs_1x cpgs_5x median_cov global_meth_pct
## <char> <fctr> <int> <int> <num> <num>
## 1: Pum1 Pum 14851783 7899658 5 19.41441
## 2: Qui1 Qui 15017682 8313533 5 18.90595
## 3: Qui2 Qui 14534059 7144141 4 21.60761
## 4: Qui3 Qui 14942618 8409699 5 20.51988
## 5: Qui4 Qui 14684326 7244626 4 19.29170
## 6: Rio1 Rio 15105950 8949697 5 20.31559
## 7: Rio2 Rio 14776205 7595578 5 20.09963
## 8: Rio3 Rio 14337226 6334602 4 21.13009
## 9: Rio4 Rio 14960910 8607331 5 20.54781
## methylated_cpg_pct CHG_meth_pct CHH_meth_pct
## <num> <num> <num>
## 1: 21.49834 0.6 0.6
## 2: 20.76845 0.8 0.8
## 3: 23.70096 0.7 0.8
## 4: 22.13247 0.7 0.8
## 5: 21.76495 0.7 0.7
## 6: 21.79080 0.6 0.7
## 7: 22.32245 0.6 0.7
## 8: 24.02625 0.8 0.8
## 9: 22.02631 0.7 0.8global %>%
pivot_longer(c(global_meth_pct, CHG_meth_pct, CHH_meth_pct),
names_to = "context", values_to = "pct") %>%
mutate(context = sub("_meth_pct", "", sub("global", "CpG", context))) %>%
ggplot(aes(population, pct, colour = population)) +
geom_jitter(width = 0.1, size = 2.5) +
facet_wrap(~context, scales = "free_y") +
labs(y = "% methylation (read-level)", title = "Methylation by context") +
theme_bw() + theme(legend.position = "none")
kruskal.test(global_meth_pct ~ population, data = global)##
## Kruskal-Wallis rank sum test
##
## data: global_meth_pct by population
## Kruskal-Wallis chi-squared = 1.2, df = 2, p-value = 0.5488Per-CpG distribution (pooled)
A mosaic methylome shows up as a bimodal distribution: most CpGs near 0%, a smaller fully methylated mode near 100%, and few intermediate sites.
ggplot(tibble(pct = cpg$pct), aes(pct)) +
geom_histogram(binwidth = 2, boundary = 0) +
labs(x = "% methylation per CpG (pooled)", y = "CpGs") + theme_bw()
class_summary <- as_tibble(table(class = cpg$class)) %>%
mutate(pct = 100 * n / sum(n))
write_tsv(class_summary, file.path(out_dir, "cpg_class_summary.tsv"))
class_summary## # A tibble: 3 x 3
## class n pct
## <chr> <int> <dbl>
## 1 unmethylated (<10%) 11196509 77.6
## 2 sparse (10-50%) 610968 4.23
## 3 methylated (>=50%) 2625633 18.2Genomic features
cpg_tab <- as_tibble(as.data.frame(mcols(cpg)[, c("feature", "pct", "class")]))
feature_summary <- composition %>%
left_join(as_tibble(table(feature = genome_cpg$feature)) %>% dplyr::rename(cpg_in_assembly = n),
by = "feature") %>%
left_join(cpg_tab %>%
group_by(feature) %>%
summarise(cpg_covered = n(),
mean_meth = mean(pct),
median_meth = median(pct),
pct_methylated = 100 * mean(class == meth_labels[3]),
pct_sparse = 100 * mean(class == meth_labels[2]),
pct_unmeth = 100 * mean(class == meth_labels[1])),
by = "feature") %>%
mutate(feature = factor(feature, feature_levels),
pct_cpg_covered = 100 * cpg_covered / cpg_in_assembly)
write_tsv(feature_summary, file.path(out_dir, "feature_summary.tsv"))
feature_summary## # A tibble: 5 x 15
## feature bp n_cpg gc_pct cpg_oe pct_genome cpg_per_kb cpg_in_assembly
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 exon 2.59e7 7.25e5 43.9 0.581 2.20 28.0 726550
## 2 intron 1.10e8 1.76e6 34.4 0.542 9.29 16.0 1758002
## 3 promoter 3.94e7 8.93e5 35.5 0.718 3.34 22.6 893067
## 4 downstream 3.74e7 7.45e5 35.1 0.647 3.17 19.9 745388
## 5 intergenic 9.68e8 2.20e7 35.9 0.706 82.0 22.7 21981174
## # i 7 more variables: cpg_covered <int>, mean_meth <dbl>, median_meth <dbl>,
## # pct_methylated <dbl>, pct_sparse <dbl>, pct_unmeth <dbl>,
## # pct_cpg_covered <dbl>cpg_tab %>%
ggplot(aes(feature, fill = class)) +
geom_bar(position = "fill") +
scale_fill_manual(values = c("grey80", "#9ecae1", "#08519c")) +
labs(y = "Proportion of covered CpGs", fill = NULL) + theme_bw()
cpg_tab %>%
ggplot(aes(feature, pct)) +
geom_violin(fill = "grey85", scale = "width") +
stat_summary(fun = mean, geom = "point", colour = "red") +
labs(y = "% methylation per CpG (pooled)") + theme_bw()
Where the methylated CpGs are
Share of the genome, of all assembly CpGs, of covered CpGs, and of methylated CpGs falling in each feature. Enrichment of methylated CpGs is tested per feature against all other features (Fisher), using covered CpGs as the background.
meth_cpg <- cpg_tab$class == meth_labels[3]
feature_share <- bind_rows(
composition %>% transmute(feature, set = "genome (bp)", n = bp),
as_tibble(table(feature = genome_cpg$feature)) %>% mutate(set = "assembly CpGs"),
as_tibble(table(feature = cpg_tab$feature)) %>% mutate(set = "covered CpGs"),
as_tibble(table(feature = cpg_tab$feature[meth_cpg])) %>% mutate(set = "methylated CpGs")
) %>%
group_by(set) %>% mutate(share = n / sum(n)) %>% ungroup() %>%
mutate(set = fct_inorder(set), feature = factor(feature, feature_levels))
ggplot(feature_share, aes(set, share, fill = feature)) +
geom_col() + scale_fill_brewer(palette = "Set2") +
labs(x = NULL, y = "Share") + theme_bw()
feature_enrichment <- map_dfr(feature_levels, function(f) {
in_f <- cpg_tab$feature == f
ft <- fisher.test(table(in_f, meth_cpg))
tibble(feature = f, odds_ratio = unname(ft$estimate),
ci_low = ft$conf.int[1], ci_high = ft$conf.int[2], p = ft$p.value)
}) %>% mutate(p_adj = p.adjust(p, "BH"))
write_tsv(feature_enrichment, file.path(out_dir, "feature_methylated_cpg_enrichment.tsv"))
feature_enrichment## # A tibble: 5 x 6
## feature odds_ratio ci_low ci_high p p_adj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 exon 2.12 2.11 2.14 0 0
## 2 intron 5.56 5.54 5.59 0 0
## 3 promoter 0.742 0.736 0.748 0 0
## 4 downstream 1.16 1.15 1.17 0 0
## 5 intergenic 0.323 0.322 0.324 0 0Consistency across oysters
Per-oyster coverage-weighted methylation in each feature. If the feature ranking holds in every animal, the pooled pattern isn’t driven by a few deep samples.
ps_feature %>%
mutate(pct = 100 * meth / cov) %>%
ggplot(aes(feature, pct, colour = population)) +
geom_point(position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.6)) +
labs(y = "% methylation (read-level, per oyster)") + theme_bw()
Repeats
# Not run until a repeat annotation for merged_out.fasta exists. Repeats overlap the
# five features, so they are reported as a second, orthogonal layer.
reps <- import(repeat_gff)
seqlevels(reps, pruning.mode = "coarse") <- seqlevels(si)
rep_class <- sub("/.*", "", coalesce(reps$Target, reps$type)) # adjust to the GFF's class field
cpg$repeat_class <- "non-repeat"
h <- findOverlaps(cpg, reps, ignore.strand = TRUE, select = "first")
cpg$repeat_class[!is.na(h)] <- rep_class[h[!is.na(h)]]
as_tibble(as.data.frame(mcols(cpg)[, c("feature", "repeat_class", "pct")])) %>%
group_by(feature, repeat_class) %>%
summarise(n = n(), mean_meth = mean(pct), .groups = "drop") %>%
arrange(feature, dplyr::desc(n))Chromosome-scale distribution
Methylation in 100 kb windows along each chromosome, against gene density. In mosaic methylomes methylation tracks gene density.
win <- tileGenome(seqlengths(si), tilewidth = window_size, cut.last.tile.in.chrom = TRUE)
win <- binnedAverage(win, coverage(gene_red), "gene_frac")
win$n_cpg_assembly <- countOverlaps(win, genome_cpg)
h <- findOverlaps(cpg, win)
win_meth <- data.table(w = subjectHits(h), meth = cpg$meth[queryHits(h)],
cov = cpg$cov[queryHits(h)], pct = cpg$pct[queryHits(h)]
)[, .(meth_pct = 100 * sum(meth) / sum(cov),
frac_methylated = mean(pct >= 50), n_cpg = .N), by = w]
windows <- as_tibble(as.data.frame(win)) %>%
mutate(w = row_number()) %>%
left_join(win_meth, by = "w") %>%
filter(n_cpg >= 20) %>%
mutate(cpg_per_kb = 1000 * n_cpg_assembly / width)
write_tsv(windows, file.path(out_dir, "window_methylation.tsv"))
windows %>%
transmute(seqnames, mb = (start + end) / 2e6,
`CpG methylation (%)` = meth_pct, `gene coverage (%)` = 100 * gene_frac) %>%
pivot_longer(-c(seqnames, mb)) %>%
ggplot(aes(mb, value)) +
geom_line(linewidth = 0.2) +
geom_smooth(method = "loess", span = 0.2, se = FALSE, colour = "red", linewidth = 0.5) +
facet_grid(name ~ seqnames, scales = "free", space = "free_x") +
labs(x = "Position (Mb)", y = NULL) + theme_bw() +
theme(strip.text.x = element_text(size = 7), axis.text.x = element_text(size = 6))
cor.test(windows$meth_pct, windows$gene_frac, method = "spearman")##
## Spearman's rank correlation rho
##
## data: windows$meth_pct and windows$gene_frac
## S = 1.3192e+11, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3651871cor.test(windows$meth_pct, windows$cpg_per_kb, method = "spearman")##
## Spearman's rank correlation rho
##
## data: windows$meth_pct and windows$cpg_per_kb
## S = 3.4353e+11, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.6531603ggplot(windows, aes(100 * gene_frac, meth_pct)) +
geom_bin2d(bins = 60) + geom_smooth(method = "loess", colour = "red") +
labs(x = "Gene coverage in window (%)", y = "Window CpG methylation (%)") + theme_bw()
Gene-body methylation
Per gene
Pooled coverage-weighted methylation over all covered CpGs in each gene, plus the same for its exons, introns, and 2 kb promoter. Genes need ≥ 5 covered CpGs. CpG O/E and GC are computed from sequence for every gene.
region_meth <- function(gr, ids, label) {
h <- findOverlaps(cpg, gr, ignore.strand = TRUE)
data.table(gene_id = ids[subjectHits(h)], meth = cpg$meth[queryHits(h)],
cov = cpg$cov[queryHits(h)], pct = cpg$pct[queryHits(h)]
)[, .(pct = 100 * sum(meth) / sum(cov), n = .N,
frac_meth_cpg = mean(pct >= 50)), by = gene_id] %>%
as_tibble() %>%
rename_with(~ paste0(label, "_", .x), -gene_id)
}
gene_promoters <- trim(promoters(genes, upstream = flank_bp, downstream = 0))
gene_introns <- psetdiff(granges(genes),
split(granges(exon_gr), factor(exon_gr$gene_id, names(genes))))
gene_introns <- unlist(gene_introns)
gene_seq <- getSeq(fa, genes)
gene_comp <- letterFrequency(gene_seq, c("C", "G")) * 1 # double: C * G overflows integer
n_exons <- table(factor(exon_gr$gene_id, names(genes)))
gene_meth <- tibble(gene_id = names(genes),
chr = as.character(seqnames(genes)),
length = width(genes),
n_exons = as.integer(n_exons[names(genes)]),
n_cpg = as.numeric(vcountPattern("CG", gene_seq)),
gc_pct = 100 * rowSums(gene_comp) / width(genes),
cpg_oe = n_cpg * width(genes) / (gene_comp[, "C"] * gene_comp[, "G"])) %>%
left_join(region_meth(genes, names(genes), "body"), by = "gene_id") %>%
left_join(region_meth(exon_gr, exon_gr$gene_id, "exon"), by = "gene_id") %>%
left_join(region_meth(gene_introns, names(gene_introns), "intron"), by = "gene_id") %>%
left_join(region_meth(gene_promoters, names(gene_promoters), "promoter"), by = "gene_id") %>%
mutate(body_n = replace_na(body_n, 0L),
body_class = cut(if_else(body_n >= min_gene_cpgs, body_pct, NA_real_),
meth_breaks, meth_labels, right = FALSE))
rm(gene_seq); invisible(gc())
gbm <- filter(gene_meth, !is.na(body_class))
message(nrow(gbm), " of ", nrow(gene_meth), " genes with >= ", min_gene_cpgs, " covered CpGs")
table(gbm$body_class)##
## unmethylated (<10%) sparse (10-50%) methylated (>=50%)
## 8103 702 4837Distribution
Bivalve gene-body methylation is usually bimodal: one group of genes heavily methylated, the other nearly unmethylated. The local minimum of the density between the modes is printed as a check on the fixed 10/50% class breaks.
dens <- density(gbm$body_pct, from = 0, to = 100)
mid <- dens$x > 10 & dens$x < 90
message("Density minimum between modes: ", round(dens$x[mid][which.min(dens$y[mid])], 1), "%")
ggplot(gbm, aes(body_pct)) +
geom_histogram(binwidth = 2, boundary = 0) +
geom_vline(xintercept = c(10, 50), linetype = 2, colour = "grey50") +
labs(x = "Gene-body CpG methylation (%)", y = "Genes") + theme_bw()
ggplot(gbm, aes(exon_pct, intron_pct)) +
geom_bin2d(bins = 60) + geom_abline(linetype = 2, colour = "red") +
labs(x = "Exon methylation (%)", y = "Intron methylation (%)") + theme_bw()
ggplot(gbm, aes(promoter_pct, body_pct)) +
geom_bin2d(bins = 60) +
labs(x = "Promoter (2 kb upstream) methylation (%)", y = "Gene-body methylation (%)") +
theme_bw()
CpG O/E and gene structure
Genes methylated in the germline lose CpGs over evolutionary time, so methylated genes are expected to have lower CpG O/E. A bimodal CpG O/E distribution across all genes (not only those with WGBS coverage) is the sequence-level signature of this.
ggplot(gene_meth, aes(cpg_oe)) +
geom_histogram(bins = 80) + xlim(0, 1.5) +
labs(x = "Gene CpG O/E (all genes)", y = "Genes") + theme_bw()
ggplot(gbm, aes(cpg_oe, fill = body_class)) +
geom_density(alpha = 0.5) + xlim(0, 1.5) +
scale_fill_manual(values = c("grey70", "#9ecae1", "#08519c")) +
labs(x = "Gene CpG O/E", fill = NULL) + theme_bw()
cor.test(gbm$body_pct, gbm$cpg_oe, method = "spearman")##
## Spearman's rank correlation rho
##
## data: gbm$body_pct and gbm$cpg_oe
## S = 6.7227e+11, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.5887704gbm %>%
transmute(body_class, `gene length (kb)` = length / 1000, `exons` = n_exons,
`GC (%)` = gc_pct) %>%
pivot_longer(-body_class) %>%
ggplot(aes(body_class, value, fill = body_class)) +
geom_boxplot(outlier.size = 0.3) + scale_y_log10() +
scale_fill_manual(values = c("grey70", "#9ecae1", "#08519c")) +
facet_wrap(~name, scales = "free_y") + labs(x = NULL) + theme_bw() +
theme(legend.position = "none", axis.text.x = element_text(angle = 20, hjust = 1))
write_tsv(gene_meth, file.path(out_dir, "gene_methylation.tsv"))Metagene profile
Methylation along a scaled gene body (TSS → TES, 20 bins) with 2 kb flanks (10 bins each), oriented by strand. Flanks can run into neighbouring genes; that is left in, as it is what a CpG near a TSS actually sees.
ext <- trim(resize(genes, width(genes) + 2 * flank_bp, fix = "center"))
h <- findOverlaps(cpg, ext, ignore.strand = TRUE)
flank_step <- flank_bp / n_flank_bins
mg <- data.table(gene_id = names(genes)[subjectHits(h)],
pos = start(cpg)[queryHits(h)],
gs = start(genes)[subjectHits(h)], ge = end(genes)[subjectHits(h)],
minus = as.character(strand(genes))[subjectHits(h)] == "-",
meth = cpg$meth[queryHits(h)], cov = cpg$cov[queryHits(h)],
pct = cpg$pct[queryHits(h)])
mg[, `:=`(rel = fifelse(minus, ge - pos, pos - gs), len = ge - gs + 1)] # 0 = TSS
mg[, bin := fcase(rel < 0, floor(rel / flank_step), # -10 .. -1
rel < len, floor(rel / len * n_body_bins), # 0 .. 19
default = n_body_bins + floor((rel - len) / flank_step))] # 20 .. 29
mg <- mg[bin >= -n_flank_bins & bin < n_body_bins + n_flank_bins]
mg <- merge(mg, as.data.table(gene_meth[, c("gene_id", "body_class")]), by = "gene_id")
metagene <- rbind(
mg[, .(body_class = "all genes", meth_pct = 100 * sum(meth) / sum(cov),
frac_methylated = mean(pct >= 50), n_cpg = .N), by = bin],
mg[!is.na(body_class), .(meth_pct = 100 * sum(meth) / sum(cov),
frac_methylated = mean(pct >= 50), n_cpg = .N),
by = .(bin, body_class = as.character(body_class))],
use.names = TRUE)
metagene[, body_class := factor(body_class, c("all genes", meth_labels))]
write_tsv(metagene, file.path(out_dir, "metagene_profile.tsv"))
ggplot(metagene, aes(bin + 0.5, meth_pct, colour = body_class)) +
geom_line() +
geom_vline(xintercept = c(0, n_body_bins), linetype = 2, colour = "grey50") +
scale_x_continuous(breaks = c(-n_flank_bins, 0, n_body_bins, n_body_bins + n_flank_bins),
labels = c(paste0("-", flank_bp / 1000, " kb"), "TSS", "TES",
paste0("+", flank_bp / 1000, " kb"))) +
scale_colour_manual(values = c("black", "grey60", "#9ecae1", "#08519c")) +
labs(x = NULL, y = "CpG methylation (%)", colour = NULL) + theme_bw()
rm(mg); invisible(gc())Exon and intron rank
Methylation by position of the exon/intron along the transcript (5’ → 3’), for multi-exon genes. In many invertebrates the first exon and 5’ end are less methylated than the rest of the body.
ex <- data.table(seqnames = as.character(seqnames(exon_gr)), start = start(exon_gr),
end = end(exon_gr), strand = as.character(strand(exon_gr)),
tx = exon_gr$tx, gene_id = exon_gr$gene_id)
setorder(ex, tx, start)
ex[, `:=`(n_ex = .N, idx = seq_len(.N)), by = tx]
ex[, rank := fifelse(strand == "-", n_ex - idx + 1L, idx)]
intr <- ex[n_ex > 1, .(seqnames = seqnames[1], start = end[-.N] + 1L, end = start[-1] - 1L,
strand = strand[1], n_ex = n_ex[1], idx = seq_len(.N - 1)), by = tx]
intr <- intr[end >= start]
intr[, rank := fifelse(strand == "-", n_ex - idx, idx)]
# Ranks 1-5, then "6+"; the final element of each transcript is pulled out as "last"
# (except a lone intron, which stays rank 1).
rank_meth <- function(dt, what) {
n_el <- if (what == "exon") dt$n_ex else dt$n_ex - 1L
label <- fifelse(dt$rank == n_el & dt$rank > 1, "last", as.character(pmin(dt$rank, 6)))
gr <- makeGRangesFromDataFrame(dt, seqinfo = si)
h <- findOverlaps(cpg, gr, ignore.strand = TRUE)
data.table(rank_label = label[subjectHits(h)],
meth = cpg$meth[queryHits(h)], cov = cpg$cov[queryHits(h)]
)[, .(meth_pct = 100 * sum(meth) / sum(cov), n_cpg = .N), by = rank_label
][, element := what][]
}
rank_tab <- rbind(rank_meth(ex[n_ex > 1], "exon"), rank_meth(intr, "intron")) %>%
mutate(rank_label = factor(recode(rank_label, `6` = "6+"),
c(as.character(1:5), "6+", "last")))
write_tsv(rank_tab, file.path(out_dir, "exon_intron_rank_methylation.tsv"))
ggplot(rank_tab, aes(rank_label, meth_pct, colour = element, group = element)) +
geom_point(size = 2) + geom_line() +
labs(x = "Rank from 5' end", y = "CpG methylation (%)", colour = NULL,
subtitle = "Multi-exon genes; 'last' = final exon / intron") + theme_bw()
Variability among oysters
Per-oyster gene-body methylation (reads summed over each gene) for genes with ≥ 30 calls in every oyster. In bivalves, heavily methylated genes tend to be the most stable among individuals.
gene_ind <- ps_gene[cov >= min_gene_calls & n_cpg >= min_gene_cpgs]
gene_ind <- gene_ind[, if (.N == nrow(samples)) .SD, by = gene_id]
gene_ind[, pct := 100 * meth / cov]
gene_var <- gene_ind[, .(mean_pct = mean(pct), sd_pct = sd(pct)), by = gene_id] %>%
as_tibble() %>%
left_join(select(gene_meth, gene_id, body_class), by = "gene_id")
write_tsv(gene_var, file.path(out_dir, "gene_methylation_variability.tsv"))
message(nrow(gene_var), " genes with per-oyster estimates in all ", nrow(samples), " oysters")
ggplot(gene_var, aes(mean_pct, sd_pct)) +
geom_bin2d(bins = 60) + geom_smooth(method = "loess", colour = "red") +
labs(x = "Mean gene-body methylation across oysters (%)",
y = "SD among oysters (% points)") + theme_bw()
# Between-oyster agreement in gene-body methylation (pairwise Spearman).
pm_gene <- dcast(gene_ind, gene_id ~ sample, value.var = "pct")
cc <- cor(as.matrix(pm_gene[, -1]), method = "spearman")
summary(cc[upper.tri(cc)])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.8497 0.8581 0.8620 0.8611 0.8651 0.8687Methylation and gene function
Functional annotation from 13.5 (DIAMOND blastp of GN
proteins vs Swiss-Prot, E ≤ 1e-8, GO from the best hit). GO terms are
the direct UniProt annotations, not propagated to parent terms, so broad
terms are under-counted; treat the term tests as descriptive.
annot <- read_tsv(annot_file, show_col_types = FALSE,
col_types = cols(.default = col_character()))
gene_fun <- gene_meth %>%
left_join(select(annot, gene_id, accession, protein_name, go_bp), by = "gene_id") %>%
mutate(annotation = case_when(is.na(accession) ~ "no Swiss-Prot hit",
is.na(go_bp) ~ "hit, no GO BP",
TRUE ~ "hit with GO BP"),
annotation = factor(annotation, c("hit with GO BP", "hit, no GO BP",
"no Swiss-Prot hit")))
write_tsv(gene_fun, file.path(out_dir, "gene_methylation_annotated.tsv"))Annotated vs unannotated genes
Genes with a Swiss-Prot hit are mostly conserved; unannotated genes are enriched for lineage-specific and mis-predicted genes. Conserved genes are expected to carry more gene-body methylation.
gf <- filter(gene_fun, !is.na(body_class))
gf %>% count(annotation, body_class) %>%
group_by(annotation) %>% mutate(pct = 100 * n / sum(n)) %>% ungroup()## # A tibble: 9 x 4
## annotation body_class n pct
## <fct> <fct> <int> <dbl>
## 1 hit with GO BP unmethylated (<10%) 2430 40.5
## 2 hit with GO BP sparse (10-50%) 306 5.1
## 3 hit with GO BP methylated (>=50%) 3264 54.4
## 4 hit, no GO BP unmethylated (<10%) 590 53.2
## 5 hit, no GO BP sparse (10-50%) 61 5.50
## 6 hit, no GO BP methylated (>=50%) 459 41.4
## 7 no Swiss-Prot hit unmethylated (<10%) 5083 77.8
## 8 no Swiss-Prot hit sparse (10-50%) 335 5.13
## 9 no Swiss-Prot hit methylated (>=50%) 1114 17.1chisq.test(table(gf$annotation, gf$body_class))##
## Pearson's Chi-squared test
##
## data: table(gf$annotation, gf$body_class)
## X-squared = 1983.4, df = 4, p-value < 2.2e-16ggplot(gf, aes(annotation, body_pct, fill = annotation)) +
geom_violin(scale = "width") + geom_boxplot(width = 0.1, outlier.shape = NA, fill = "white") +
labs(x = NULL, y = "Gene-body methylation (%)") + theme_bw() +
theme(legend.position = "none")
GO term tests
Background: genes with a gene-body estimate and ≥ 1 GO BP term. Two tests per term (10–500 genes):
- Wilcoxon rank-sum of gene-body methylation, term
genes vs the rest of the background (a continuous, direction-aware test;
delta_median> 0 = term genes more methylated). - Fisher for over-representation among methylated (≥ 50%) and among unmethylated (< 10%) genes.
go_long <- gf %>%
filter(!is.na(go_bp)) %>%
select(gene_id, body_pct, body_class, go_bp) %>%
separate_rows(go_bp, sep = "; ") %>%
mutate(go_id = str_extract(go_bp, "GO:\\d{7}"),
term = str_trim(str_remove(go_bp, "\\s*\\[GO:\\d{7}\\]$"))) %>%
filter(!is.na(go_id)) %>%
distinct(gene_id, go_id, .keep_all = TRUE)
bg <- distinct(go_long, gene_id, body_pct, body_class)
is_high <- set_names(bg$body_class == meth_labels[3], bg$gene_id)
is_low <- set_names(bg$body_class == meth_labels[1], bg$gene_id)
terms <- go_long %>% count(go_id, term) %>%
filter(n >= min_go_genes, n <= max_go_genes)
go_tests <- map_dfr(seq_len(nrow(terms)), function(k) {
ids <- go_long$gene_id[go_long$go_id == terms$go_id[k]]
in_t <- bg$gene_id %in% ids
wt <- wilcox.test(bg$body_pct[in_t], bg$body_pct[!in_t])
fisher_p <- function(flag) fisher.test(table(factor(in_t, c(TRUE, FALSE)),
factor(flag, c(TRUE, FALSE))),
alternative = "greater")$p.value
tibble(go_id = terms$go_id[k], term = terms$term[k], n_genes = sum(in_t),
median_meth = median(bg$body_pct[in_t]),
delta_median = median(bg$body_pct[in_t]) - median(bg$body_pct[!in_t]),
pct_methylated_genes = 100 * mean(is_high[in_t]),
wilcox_p = wt$p.value,
high_fisher_p = fisher_p(is_high),
low_fisher_p = fisher_p(is_low))
}) %>%
mutate(across(ends_with("_p"), ~ p.adjust(.x, "BH"), .names = "{.col}_adj")) %>%
arrange(wilcox_p)
write_tsv(go_tests, file.path(out_dir, "GO_BP_methylation_tests.tsv"))
message(nrow(bg), " background genes; ", nrow(terms), " GO BP terms tested; overall ",
round(100 * mean(is_high), 1), "% of background genes methylated")
go_tests %>% filter(wilcox_p_adj < 0.05) %>% count(direction = sign(delta_median))## # A tibble: 2 x 2
## direction n
## <dbl> <int>
## 1 -1 195
## 2 1 126top_go <- go_tests %>%
filter(wilcox_p_adj < 0.05) %>%
group_by(direction = if_else(delta_median > 0, "higher methylation", "lower methylation")) %>%
slice_min(wilcox_p, n = 15, with_ties = FALSE) %>% ungroup()
if (nrow(top_go) > 0) {
top_go %>%
mutate(term = fct_reorder(str_trunc(term, 60), median_meth)) %>%
ggplot(aes(median_meth, term, size = n_genes, colour = -log10(wilcox_p_adj))) +
geom_point() +
geom_vline(xintercept = median(bg$body_pct), linetype = 2) +
facet_grid(direction ~ ., scales = "free_y", space = "free_y") +
scale_colour_viridis_c() +
labs(x = "Median gene-body methylation of term genes (%)", y = NULL,
subtitle = "Dashed line: background median") + theme_bw()
}
Broad functional groups
A coarse view, grouping genes by keywords in their GO BP term names. The groups follow the housekeeping-vs-inducible contrast often reported for bivalve gene-body methylation (constitutively expressed, housekeeping functions methylated; environmentally responsive functions less so). Genes can fall into more than one group; edit the patterns to suit.
func_groups <- c(
"Translation / ribosome" = "translation|ribosom|rRNA processing",
"RNA processing / splicing" = "mRNA processing|RNA splicing|spliceosom",
"DNA replication / repair" = "DNA replication|DNA repair",
"Cell cycle" = "cell cycle|cell division|mitotic",
"Mitochondria / energy" = "mitochondri|electron transport|tricarboxylic|ATP synthesis",
"Protein degradation" = "ubiquitin|proteasom|protein catabolic",
"Signal transduction" = "signal transduction|signaling pathway|G protein-coupled",
"Immune / defense" = "immune|defense response|inflammatory|innate",
"Stress response" = "response to stress|oxidative stress|heat|hypoxia|unfolded protein",
"Development / differentiation" = "development|morphogenesis|differentiation")
group_meth <- imap_dfr(func_groups, function(pat, grp) {
go_long %>% filter(str_detect(term, regex(pat, ignore_case = TRUE))) %>%
distinct(gene_id, body_pct) %>% mutate(group = grp)
})
group_stats <- group_meth %>%
group_by(group) %>%
summarise(n_genes = n(), median_meth = median(body_pct),
pct_methylated_genes = 100 * mean(body_pct >= 50),
wilcox_p = wilcox.test(body_pct, bg$body_pct[!bg$gene_id %in% gene_id])$p.value) %>%
mutate(wilcox_p_adj = p.adjust(wilcox_p, "BH")) %>%
arrange(dplyr::desc(median_meth))
write_tsv(group_stats, file.path(out_dir, "functional_group_methylation.tsv"))
group_stats## # A tibble: 10 x 6
## group n_genes median_meth pct_methylated_genes wilcox_p wilcox_p_adj
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 RNA processin~ 172 85.9 83.1 3.04e-19 1.01e-18
## 2 Translation /~ 325 84.0 83.1 5.59e-33 5.59e-32
## 3 Cell cycle 438 83.3 69.4 1.80e-20 9.02e-20
## 4 Mitochondria ~ 256 82.0 80.5 4.26e-16 7.10e-16
## 5 DNA replicati~ 213 81.8 73.2 1.88e-11 2.35e-11
## 6 Protein degra~ 566 80.4 69.4 8.22e-18 2.06e-17
## 7 Stress respon~ 297 66.6 51.2 2.85e- 1 2.85e- 1
## 8 Development /~ 1468 27.2 46.6 3.14e- 9 3.48e- 9
## 9 Signal transd~ 1390 14.4 45.7 3.63e-16 7.10e-16
## 10 Immune / defe~ 467 1.22 37.7 2.21e-14 3.15e-14group_meth %>%
mutate(group = factor(group, group_stats$group)) %>%
ggplot(aes(body_pct, group)) +
geom_boxplot(outlier.size = 0.3) +
geom_vline(xintercept = median(bg$body_pct), linetype = 2) +
labs(x = "Gene-body methylation (%)", y = NULL,
subtitle = "Dashed line: median of all genes with GO BP") + theme_bw()
Summary
summary_tab <- tribble(
~metric, ~value,
"CpGs in assembly", length(genome_cpg),
"CpGs passing pooled filters", length(cpg),
"mean global CpG methylation, per oyster (%)", mean(global$global_meth_pct),
"methylated CpGs (>=50%, pooled) (%)", 100 * mean(cpg$pct >= 50),
"genes with gene-body estimate", nrow(gbm),
"methylated genes (>=50%) (%)", 100 * mean(gbm$body_class == meth_labels[3]),
"unmethylated genes (<10%) (%)", 100 * mean(gbm$body_class == meth_labels[1]),
"rho gene methylation vs CpG O/E", cor(gbm$body_pct, gbm$cpg_oe, method = "spearman",
use = "complete.obs"),
"rho window methylation vs gene coverage", cor(windows$meth_pct, windows$gene_frac,
method = "spearman"),
"GO BP terms tested", nrow(go_tests),
"GO BP terms, Wilcoxon FDR < 0.05", sum(go_tests$wilcox_p_adj < 0.05)
) %>% mutate(value = signif(value, 4))
write_tsv(summary_tab, file.path(out_dir, "landscape_summary.tsv"))
summary_tab## # A tibble: 11 x 2
## metric value
## <chr> <dbl>
## 1 CpGs in assembly 26100000
## 2 CpGs passing pooled filters 14430000
## 3 mean global CpG methylation, per oyster (%) 20.2
## 4 methylated CpGs (>=50%, pooled) (%) 18.2
## 5 genes with gene-body estimate 13640
## 6 methylated genes (>=50%) (%) 35.5
## 7 unmethylated genes (<10%) (%) 59.4
## 8 rho gene methylation vs CpG O/E -0.589
## 9 rho window methylation vs gene coverage 0.365
## 10 GO BP terms tested 839
## 11 GO BP terms, Wilcoxon FDR < 0.05 321Next steps
- Repeats / TEs. 69% of the assembly is repeat. Run
RepeatMasker (or EDTA) on
merged_out.fasta, setrepeat_gff, and the repeat chunk will split each feature into repeat vs non-repeat; intronic TEs are a common source of gene-body methylation. - Expression. The classic test of the gene-body model
is methylation vs expression level and breadth. Join
gene_methylation.tsvto per-gene expression from09(mean and CV of normalized counts across the 12 RNA-seq oysters) on GN gene ID. These are different animals from the WGBS oysters, so this is a gene-level, not individual-level, comparison. - GO propagation. With
GO.db(ortopGO), propagate annotations to ancestors before the term tests; that will mainly add power to the broad terms. - SNP filter. C/T SNPs at CpGs read as unmethylated.
For the landscape this mainly lowers the methylated fraction slightly;
it matters more for
13.4.
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] Rsamtools_2.22.0 Biostrings_2.74.1 XVector_0.46.0
## [4] rtracklayer_1.66.0 GenomicRanges_1.58.0 GenomeInfoDb_1.42.3
## [7] IRanges_2.40.1 S4Vectors_0.44.0 BiocGenerics_0.52.0
## [10] data.table_1.17.0 lubridate_1.9.4 forcats_1.0.0
## [13] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.4
## [16] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
## [19] ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2
## [3] farver_2.1.2 bitops_1.0-9
## [5] fastmap_1.2.0 RCurl_1.98-1.19
## [7] GenomicAlignments_1.42.0 XML_3.99-0.18
## [9] digest_0.6.37 timechange_0.3.0
## [11] lifecycle_1.0.4 magrittr_2.0.3
## [13] compiler_4.4.3 rlang_1.1.5
## [15] sass_0.4.9 tools_4.4.3
## [17] utf8_1.2.4 yaml_2.3.10
## [19] knitr_1.50 S4Arrays_1.6.0
## [21] labeling_0.4.3 bit_4.6.0
## [23] curl_6.2.2 DelayedArray_0.32.0
## [25] RColorBrewer_1.1-3 abind_1.4-8
## [27] BiocParallel_1.40.2 withr_3.0.2
## [29] grid_4.4.3 colorspace_2.1-1
## [31] scales_1.3.0 SummarizedExperiment_1.36.0
## [33] cli_3.6.4 rmarkdown_2.29
## [35] crayon_1.5.3 generics_0.1.3
## [37] httr_1.4.7 tzdb_0.5.0
## [39] rjson_0.2.23 cachem_1.1.0
## [41] zlibbioc_1.52.0 splines_4.4.3
## [43] parallel_4.4.3 restfulr_0.0.17
## [45] matrixStats_1.5.0 vctrs_0.6.5
## [47] Matrix_1.7-3 jsonlite_2.0.0
## [49] hms_1.1.3 bit64_4.6.0-1
## [51] jquerylib_0.1.4 glue_1.8.0
## [53] codetools_0.2-20 stringi_1.8.7
## [55] gtable_0.3.6 BiocIO_1.16.0
## [57] UCSC.utils_1.2.0 munsell_0.5.1
## [59] pillar_1.10.2 htmltools_0.5.8.1
## [61] GenomeInfoDbData_1.2.13 R6_2.6.1
## [63] vroom_1.6.5 evaluate_1.0.3
## [65] Biobase_2.66.0 lattice_0.22-7
## [67] bslib_0.9.0 SparseArray_1.6.2
## [69] nlme_3.1-168 mgcv_1.9-3
## [71] xfun_0.52 MatrixGenerics_1.18.1
## [73] pkgconfig_2.0.3