13.4 - WGBS DNA methylation across three populations
2026-10-08
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 RPaths 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.
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.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 <- 8Samples
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.
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)## # A tibble: 8 x 3
## sample population site
## <chr> <fct> <chr>
## 1 Qui1 Qui Quilhua
## 2 Qui2 Qui Quilhua
## 3 Qui3 Qui Quilhua
## 4 Qui4 Qui Quilhua
## 5 Rio1 Rio Rio_Pudeto
## 6 Rio2 Rio Rio_Pudeto
## 7 Rio3 Rio Rio_Pudeto
## 8 Rio4 Rio Rio_PudetoAlignment 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%.
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## # A tibble: 9 x 8
## sample pairs_analysed mapping_pct dup_pct CpG_meth_pct CHG_meth_pct
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Pum1 38376265 65.7 7.01 19.8 0.6
## 2 Qui1 38720614 67.2 7.01 19.2 0.8
## 3 Qui2 36627727 66.9 6.4 22 0.7
## 4 Qui3 41446330 67.4 7.12 20.9 0.7
## 5 Qui4 34443887 67.3 5.3 19.6 0.7
## 6 Rio1 41519830 68.5 6.4 20.6 0.6
## 7 Rio2 35781535 67.5 6.18 20.4 0.7
## 8 Rio3 32773197 67.9 4.81 21.5 0.8
## 9 Rio4 39972937 68.2 4.99 20.8 0.7
## # i 2 more variables: CHH_meth_pct <dbl>, population <chr>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.
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.
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.")
}## # A tibble: 8 x 4
## sample cpgs masked pct_masked
## <chr> <int> <int> <dbl>
## 1 Qui1 15436617 47599 0.308
## 2 Qui2 15005656 47599 0.317
## 3 Qui3 15294522 47603 0.311
## 4 Qui4 15156835 47597 0.314
## 5 Rio1 15472656 47606 0.308
## 6 Rio2 15203272 47601 0.313
## 7 Rio3 14877535 47599 0.32
## 8 Rio4 15404322 47600 0.309# 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## # A tibble: 8 x 6
## sample cpgs_1x cpgs_5x cpgs_10x median_cov global_meth_pct
## <chr> <int> <int> <int> <dbl> <dbl>
## 1 Qui1 15389018 9677048 2590550 6 19.3
## 2 Qui2 14958057 8272322 2566791 5 22.1
## 3 Qui3 15246919 9679776 3329617 6 20.9
## 4 Qui4 15109238 8476831 2034932 5 19.7
## 5 Rio1 15425050 10256709 3856697 6 20.7
## 6 Rio2 15155671 8825006 2548738 5 20.5
## 7 Rio3 14829936 7468570 1930883 5 21.5
## 8 Rio4 15356722 9946247 3450767 6 20.9ggplot(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
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 filteredDistribution
Invertebrate methylomes are typically mosaic/bimodal: most CpGs near 0% or near 100%, with methylation concentrated in gene bodies.
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.
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.
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)##
## exon intron promoter intergenic
## 136346 692944 280255 5126410Differential 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.
pairs <- combn(levels(samples$population), 2, simplify = FALSE)
contrasts <- set_names(pairs, map_chr(pairs, paste, collapse = "_vs_"))
names(contrasts)## [1] "Qui_vs_Rio"# 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)
dml <- imap(contrasts, function(g, name) {
d <- run_pair(meth_unite, g)
write_hits(d, name, "DML")
})
map_int(dml, length)## Qui_vs_Rio
## 14551 kb tiles (DMRs)
Tiling pools reads across neighbouring CpGs, which is the better-powered test at this depth.
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## # A tibble: 1 x 8
## contrast n_dmr hyper hypo exon intron promoter intergenic
## <chr> <int> <int> <int> <int> <int> <int> <int>
## 1 Qui_vs_Rio 379 200 179 31 31 9 308Gene-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.
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()
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)## Qui_vs_Rio
## 2sig <- 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_*.tsvto the09limma tables ongene_id(GN gene IDs; the09StringTie 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), usingSPannot-GO.tsvvia the03join.
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 pheatmap_1.0.12 limma_3.62.2
## [13] rtracklayer_1.66.0 methylKit_1.32.1 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] bitops_1.0-9 rlang_1.1.5
## [3] magrittr_2.0.3 matrixStats_1.5.0
## [5] compiler_4.4.3 mgcv_1.9-3
## [7] vctrs_0.6.5 reshape2_1.4.4
## [9] pkgconfig_2.0.3 crayon_1.5.3
## [11] fastmap_1.2.0 XVector_0.46.0
## [13] labeling_0.4.3 utf8_1.2.4
## [15] Rsamtools_2.22.0 rmarkdown_2.29
## [17] tzdb_0.5.0 UCSC.utils_1.2.0
## [19] bit_4.6.0 xfun_0.52
## [21] zlibbioc_1.52.0 cachem_1.1.0
## [23] jsonlite_2.0.0 DelayedArray_0.32.0
## [25] BiocParallel_1.40.2 parallel_4.4.3
## [27] R6_2.6.1 bslib_0.9.0
## [29] stringi_1.8.7 RColorBrewer_1.1-3
## [31] jquerylib_0.1.4 numDeriv_2016.8-1.1
## [33] Rcpp_1.0.14 SummarizedExperiment_1.36.0
## [35] knitr_1.50 R.utils_2.13.0
## [37] Matrix_1.7-3 splines_4.4.3
## [39] timechange_0.3.0 tidyselect_1.2.1
## [41] qvalue_2.38.0 abind_1.4-8
## [43] yaml_2.3.10 codetools_0.2-20
## [45] curl_6.2.2 lattice_0.22-7
## [47] plyr_1.8.9 Biobase_2.66.0
## [49] withr_3.0.2 coda_0.19-4.1
## [51] evaluate_1.0.3 mclust_6.1.3
## [53] Biostrings_2.74.1 pillar_1.10.2
## [55] MatrixGenerics_1.18.1 generics_0.1.3
## [57] vroom_1.6.5 RCurl_1.98-1.19
## [59] emdbook_1.3.14 hms_1.1.3
## [61] munsell_0.5.1 scales_1.3.0
## [63] gtools_3.9.5 glue_1.8.0
## [65] tools_4.4.3 BiocIO_1.16.0
## [67] data.table_1.17.0 GenomicAlignments_1.42.0
## [69] mvtnorm_1.4-1 XML_3.99-0.18
## [71] grid_4.4.3 bbmle_1.0.25.1
## [73] bdsmatrix_1.3-7 colorspace_2.1-1
## [75] nlme_3.1-168 GenomeInfoDbData_1.2.13
## [77] restfulr_0.0.17 cli_3.6.4
## [79] fastseg_1.52.0 S4Arrays_1.6.0
## [81] gtable_0.3.6 R.methodsS3_1.8.2
## [83] sass_0.4.9 digest_0.6.37
## [85] SparseArray_1.6.2 farver_2.1.2
## [87] rjson_0.2.23 htmltools_0.5.8.1
## [89] R.oo_1.27.1 lifecycle_1.0.4
## [91] httr_1.4.7 statmod_1.5.2
## [93] bit64_4.6.0-1 MASS_7.3-65