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. Chosen with the ref render
parameter. merged (default) uses
merged_out.fasta + GN.gene.gff3 and
output/13-wgbs/. hapa uses the
WGBS_REF=hapa alignment to
Och_HapA_assembly.fa (output/13-wgbs-hapa/),
with the GN gene models lifted onto HapA by 13.11 (same
gene IDs, so 13.5 applies) and a HapA SNP mask from
13.7 run with WGBS_REF=hapa; results go to
output/13-wgbs-hapa/13.6-landscape:
rmarkdown::render("13.6-wgbs-methylation-landscape.Rmd", params = list(ref = "hapa")).
Either way 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)ref <- match.arg(params$ref, c("merged", "hapa"))
ref_cfg <- list(
merged = list(out_root = "../output/13-wgbs", genome_dir = "../data/13-wgbs-genome",
fasta = "merged_out.fasta", gff = "GN.gene.gff3"),
hapa = list(out_root = "../output/13-wgbs-hapa", genome_dir = "../data/13-wgbs-genome-hapa",
fasta = "Och_HapA_assembly.fa", gff = "GN.gene.hapa.gff3") # gff from 13.11
)[[ref]]
out_root <- ref_cfg$out_root
out_dir <- file.path(out_root, "13.6-landscape")
meth_dir <- file.path(out_root, "methylation")
genome_dir <- ref_cfg$genome_dir
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 <- file.path(out_root, "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, ref_cfg$fasta))
gff <- import(file.path(genome_dir, ref_cfg$gff))
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 21302904 589293 43.8 0.576 2.70 27.7
## 2 intron 90860377 1451415 34.4 0.541 11.5 16.0
## 3 promoter 32576339 735676 35.5 0.716 4.13 22.6
## 4 downstream 30660238 609162 35.1 0.645 3.88 19.9
## 5 intergenic 614230198 13819898 35.8 0.703 77.8 22.5Read 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 17208692 100
## 2 covered in >= 1 oyster 16516489 96
## 3 pass pooled filters 15288213 88.8# 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(file.path(out_root, "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 15116179 9139093 6 19.87521
## 2: Qui1 Qui 15389019 9677048 6 19.28576
## 3: Qui2 Qui 14958058 8272322 5 22.07488
## 4: Qui3 Qui 15246920 9679776 6 20.94474
## 5: Qui4 Qui 15109239 8476832 5 19.72137
## 6: Rio1 Rio 15425051 10256709 6 20.69646
## 7: Rio2 Rio 15155672 8825006 5 20.46411
## 8: Rio3 Rio 14829937 7468570 5 21.54410
## 9: Rio4 Rio 15356723 9946247 6 20.88673
## methylated_cpg_pct CHG_meth_pct CHH_meth_pct
## <num> <num> <num>
## 1: 21.82387 0.6 0.6
## 2: 21.04088 0.8 0.8
## 3: 23.91153 0.7 0.8
## 4: 22.43464 0.7 0.8
## 5: 22.04242 0.7 0.7
## 6: 22.12881 0.6 0.7
## 7: 22.56393 0.7 0.7
## 8: 24.17103 0.8 0.8
## 9: 22.25073 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 = 0.86667, df = 2, p-value = 0.6483Per-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%) 11642750 76.2
## 2 sparse (10-50%) 754584 4.94
## 3 methylated (>=50%) 2890879 18.9Genomic 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.13e7 5.89e5 43.8 0.576 2.70 27.7 590542
## 2 intron 9.09e7 1.45e6 34.4 0.541 11.5 16.0 1452754
## 3 promoter 3.26e7 7.36e5 35.5 0.716 4.13 22.6 735851
## 4 downstream 3.07e7 6.09e5 35.1 0.645 3.88 19.9 609316
## 5 intergenic 6.14e8 1.38e7 35.8 0.703 77.8 22.5 13820229
## # 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.04 2.03 2.05 0 0
## 2 intron 5.24 5.22 5.25 0 0
## 3 promoter 0.743 0.738 0.748 0 0
## 4 downstream 1.12 1.12 1.13 2.48e-262 2.48e-262
## 5 intergenic 0.344 0.343 0.345 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 = 4.9532e+10, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3969998cor.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 = 1.4476e+11, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.7622579ggplot(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%)
## 9638 881 5644Distribution
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 = 1.1029e+12, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.5671771gbm %>%
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.8562 0.8620 0.8663 0.8649 0.8681 0.8720Methylation 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%) 2681 40.7
## 2 hit with GO BP sparse (10-50%) 267 4.05
## 3 hit with GO BP methylated (>=50%) 3639 55.2
## 4 hit, no GO BP unmethylated (<10%) 681 53.6
## 5 hit, no GO BP sparse (10-50%) 65 5.11
## 6 hit, no GO BP methylated (>=50%) 525 41.3
## 7 no Swiss-Prot hit unmethylated (<10%) 6276 75.6
## 8 no Swiss-Prot hit sparse (10-50%) 549 6.61
## 9 no Swiss-Prot hit methylated (>=50%) 1480 17.8chisq.test(table(gf$annotation, gf$body_class))##
## Pearson's Chi-squared test
##
## data: table(gf$annotation, gf$body_class)
## X-squared = 2291.3, 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 240
## 2 1 139top_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~ 182 84.7 90.1 1.48e-24 7.41e-24
## 2 Translation /~ 327 84.4 88.4 2.04e-41 2.04e-40
## 3 Cell cycle 474 81.5 70.9 5.98e-21 9.97e-21
## 4 Mitochondria ~ 292 80.9 85.3 4.10e-22 1.10e-21
## 5 DNA replicati~ 247 80.1 69.2 4.78e- 8 5.31e- 8
## 6 Protein degra~ 620 78.6 69.4 1.47e-15 2.10e-15
## 7 Stress respon~ 321 63.2 51.4 1.16e- 1 1.16e- 1
## 8 Development /~ 1607 31.8 47.9 2.87e-11 3.59e-11
## 9 Signal transd~ 1526 4.23 45.7 4.42e-22 1.10e-21
## 10 Immune / defe~ 523 1.14 36.1 5.93e-22 1.19e-21group_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 17210000
## 2 CpGs passing pooled filters 15290000
## 3 mean global CpG methylation, per oyster (%) 20.6
## 4 methylated CpGs (>=50%, pooled) (%) 18.9
## 5 genes with gene-body estimate 16160
## 6 methylated genes (>=50%) (%) 34.9
## 7 unmethylated genes (<10%) (%) 59.6
## 8 rho gene methylation vs CpG O/E -0.567
## 9 rho window methylation vs gene coverage 0.397
## 10 GO BP terms tested 960
## 11 GO BP terms, Wilcoxon FDR < 0.05 379Next 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