--- title: "07 - Genomic annotation (exploratory)" output: html_document --- Place CpG methylation, and the DML/DMR lists from step 6, in genomic features: exons, introns, promoters, intergenic regions and transposable elements. **This step is exploratory.** Step 6 found no evidence of an oil effect: the real oil vs. control split ranked 9 of 10 among all possible splits, for both DMLs and DMRs (see `plan.md`). So: - the main output here is a description of the methylation landscape itself, which does not depend on treatment; - the DML/DMR lists are always compared against the lists from the 9 random splits, not just against background. If the real list looks like the random ones, it says nothing about oil. - Input: `output/06-differential-methylation/` (CpG methylation, tested CpGs/tiles, DML/DMR lists, permutation lists), `data/genome/` (FASTA, GFF, RepeatMasker, GO annotations from step 3), `output/03-genome-prep/chrom-sizes.tsv` - Output feature tracks: `data/annotation/` (gitignored) - Output summaries: `output/07-annotation/` Uses bedtools from the existing `myflow` env, and data.table, GenomicRanges, Biostrings and ggplot2 from the lab R container. ```{r setup, include=FALSE} # Paths are relative to code/, so chunks work the same run one at a time or knitted knitr::opts_chunk$set(echo = TRUE, eval = TRUE) library(data.table) library(GenomicRanges) library(ggplot2) ``` # Feature tracks Built from NCBI's RefSeq annotation (GFF) and RepeatMasker output for C_virginica-3.0. - **genes**: every `gene` and `pseudogene` feature (protein-coding, lncRNA, tRNA, ...). - **exons**: all exons of all transcripts, merged. - **introns**: gene bodies minus exons. - **promoters**: 1 kb upstream of each gene's start, strand-aware. - **intergenic**: everything that is not a gene body or promoter. - **TEs**: RepeatMasker elements other than simple repeats, low-complexity sequence, satellites and RNA genes. NCBI's RepeatMasker run for this assembly classifies almost all of these as `DNA?` (unclassified DNA transposons), so this track is incomplete and coarse. CpGs get one exclusive category, by priority exon > intron > promoter > intergenic, plus a separate TE flag. ```{bash features} set -euo pipefail export PATH=/mmfs1/gscratch/srlab/sr320/miniforge3/envs/myflow/bin:$PATH gff=../data/genome/GCF_002022765.2_C_virginica-3.0_genomic.gff.gz rm_out=../data/genome/GCF_002022765.2_C_virginica-3.0_rm.out.gz out=../data/annotation mkdir -p "${out}" # bedtools needs the genome file and BEDs in the same chromosome order sort -k1,1 ../output/03-genome-prep/chrom-sizes.tsv > "${out}/genome.txt" # Genes (BED6, name = LOC id) and a gene table (LOC id, NCBI GeneID, biotype) zcat "${gff}" | awk -F'\t' 'BEGIN { OFS = "\t" } $3 == "gene" || $3 == "pseudogene" { id = "NA"; gid = "NA"; bt = "NA" if (match($9, /ID=gene-[^;]+/)) id = substr($9, RSTART + 8, RLENGTH - 8) if (match($9, /GeneID:[0-9]+/)) gid = substr($9, RSTART + 7, RLENGTH - 7) if (match($9, /gene_biotype=[^;]+/)) bt = substr($9, RSTART + 13, RLENGTH - 13) print $1, $4 - 1, $5, id, 0, $7 > "/dev/stdout" print id, gid, bt > "/dev/stderr" }' 2> "${out}/genes-info.tsv" | sort -k1,1 -k2,2n > "${out}/genes.bed" # Product description of each gene, from its first transcript zcat "${gff}" | awk -F'\t' 'BEGIN { OFS = "\t" } $3 ~ /^(mRNA|lnc_RNA|transcript|tRNA|rRNA)$/ && match($9, /gene=[^;]+/) { g = substr($9, RSTART + 5, RLENGTH - 5) if (!(g in seen) && match($9, /product=[^;]+/)) { seen[g] = 1; print g, substr($9, RSTART + 8, RLENGTH - 8) } }' > "${out}/genes-product.tsv" zcat "${gff}" | awk -F'\t' 'BEGIN { OFS = "\t" } $3 == "exon" { print $1, $4 - 1, $5 }' \ | sort -k1,1 -k2,2n | bedtools merge > "${out}/exons.bed" bedtools merge -i "${out}/genes.bed" > "${out}/gene-bodies.bed" bedtools subtract -a "${out}/gene-bodies.bed" -b "${out}/exons.bed" > "${out}/introns.bed" bedtools flank -i "${out}/genes.bed" -g "${out}/genome.txt" -l 1000 -r 0 -s \ | sort -k1,1 -k2,2n > "${out}/promoters.bed" cat "${out}/gene-bodies.bed" <(cut -f1-3 "${out}/promoters.bed") | sort -k1,1 -k2,2n | bedtools merge \ | bedtools complement -i - -g "${out}/genome.txt" > "${out}/intergenic.bed" # RepeatMasker .out: 3 header lines; columns 5-7 = sequence, start, end (1-based); 9 = strand (+ or C); 10-11 = repeat, class/family zcat "${rm_out}" | awk 'BEGIN { OFS = "\t" } NR > 3 && $11 !~ /^(Simple_repeat|Low_complexity|Satellite|rRNA|tRNA|snRNA|scRNA|srpRNA)/ { print $5, $6 - 1, $7, $10, $11, ($9 == "C" ? "-" : "+") }' | sort -k1,1 -k2,2n > "${out}/tes.bed" for f in genes exons introns promoters intergenic tes; do printf "%-11s %8d features %8.1f Mb\n" "${f}" "$(wc -l < "${out}/${f}.bed")" \ "$(awk '{ s += $3 - $2 } END { print s / 1e6 }' "${out}/${f}.bed")" done ``` # Helpers ```{r helpers} out_dir <- "../output/07-annotation" in_dir <- "../output/06-differential-methylation" ann_dir <- "../data/annotation" dir.create(out_dir, showWarnings = FALSE) categories <- c("exon", "intron", "promoter", "intergenic") read_bed_gr <- function(file) { b <- fread(file, header = FALSE, select = 1:3, col.names = c("chr", "start", "end")) GRanges(b$chr, IRanges(b$start + 1, b$end)) # BED is 0-based, half-open } feature_files <- c(exon = "exons", intron = "introns", promoter = "promoters", intergenic = "intergenic", tes = "tes", `gene-bodies` = "gene-bodies") features <- lapply(feature_files, function(f) read_bed_gr(file.path(ann_dir, paste0(f, ".bed")))) # One exclusive category per range (exon > intron > promoter > intergenic) plus a TE flag annotate <- function(gr) { cat <- rep(NA_character_, length(gr)) for (f in categories) cat[is.na(cat) & overlapsAny(gr, features[[f]])] <- f data.table(category = factor(cat, levels = categories), in_TE = overlapsAny(gr, features$tes)) } # Hit lists written by step 6 (BED) read_hits <- function(file) { if (file.size(file) == 0) return(GRanges()) read_bed_gr(file) } samples <- c("NB3", "NB6", "NB11", "HB16", "HB2") real_split <- "HB16+HB2" perm_files <- function(prefix) { files <- list.files(file.path(in_dir, "permutations"), paste0("^", prefix, "_.*\\.bed$"), full.names = TRUE) setNames(files, sub(paste0("^", prefix, "_(.*)\\.bed$"), "\\1", basename(files))) } ``` # Methylation landscape Where do CpGs in the genome fall, where do the CpGs that MBD-BS covers fall, and how methylated is each feature? MBD enrichment pulls down methylated DNA, so the covered CpGs should be concentrated in gene bodies, which are the main methylated compartment in oyster and most invertebrates. ```{r genome-cpgs} genome <- Biostrings::readDNAStringSet("../data/genome/bismark/C_virginica-3.0.fa") names(genome) <- sub(" .*", "", names(genome)) cg <- Biostrings::vmatchPattern("CG", genome) starts <- Biostrings::startIndex(cg) genome_cpgs <- GRanges(rep(names(genome), lengths(starts)), IRanges(unlist(starts, use.names = FALSE), width = 1)) rm(genome, cg, starts); invisible(gc()) cat("CpGs in the genome:", format(length(genome_cpgs), big.mark = ","), "\n") genome_ann <- annotate(genome_cpgs) ``` ```{r landscape} meth <- fread(file.path(in_dir, "cpg-methylation-5x.csv.gz")) meth[, mean_meth := rowMeans(.SD), .SDcols = samples] meth_gr <- GRanges(meth$chr, IRanges(meth$start, width = 1)) meth <- cbind(meth, annotate(meth_gr)) landscape <- merge( genome_ann[, .(genome_CpGs = .N), by = category][, genome_pct := round(100 * genome_CpGs / sum(genome_CpGs), 1)], meth[, .(tested_CpGs = .N, mean_meth = round(mean(mean_meth), 1), pct_over_50 = round(100 * mean(mean_meth > 50), 1)), by = category][ , tested_pct := round(100 * tested_CpGs / sum(tested_CpGs), 1)], by = "category") landscape[, enrichment := round(tested_pct / genome_pct, 2)] setcolorder(landscape, c("category", "genome_CpGs", "genome_pct", "tested_CpGs", "tested_pct", "enrichment", "mean_meth", "pct_over_50")) te <- data.table( set = c("genome", "tested (5x)"), pct_in_TE = round(100 * c(mean(genome_ann$in_TE), mean(meth$in_TE)), 1), mean_meth_in_TE = c(NA, round(meth[in_TE == TRUE, mean(mean_meth)], 1)), mean_meth_outside_TE = c(NA, round(meth[in_TE == FALSE, mean(mean_meth)], 1)) ) fwrite(landscape, file.path(out_dir, "landscape-by-feature.csv")) fwrite(te, file.path(out_dir, "landscape-TEs.csv")) landscape te ``` `enrichment` > 1 means MBD-BS covers that feature more than its share of genome CpGs. ```{r landscape-plots, fig.width = 10, fig.height = 4} share <- melt(landscape[, .(category, genome = genome_pct, `tested (5x)` = tested_pct)], id.vars = "category", variable.name = "set", value.name = "pct") p1 <- ggplot(share, aes(set, pct, fill = category)) + geom_col() + labs(x = NULL, y = "% of CpGs", fill = "Feature", title = "Where CpGs fall") + theme_bw() p2 <- ggplot(meth, aes(category, mean_meth)) + geom_violin(fill = "grey80", scale = "width") + geom_boxplot(width = 0.12, outlier.shape = NA) + labs(x = NULL, y = "Mean % methylation (5 samples)", title = "Methylation by feature (CpGs >=5x)") + theme_bw() ggsave(file.path(out_dir, "landscape-share.png"), p1, width = 5, height = 4, dpi = 150) ggsave(file.path(out_dir, "landscape-methylation.png"), p2, width = 5, height = 4, dpi = 150) p1 p2 ``` ## Gene-body methylation Invertebrate genes tend to be either heavily methylated or unmethylated across the gene body, giving a bimodal distribution. MBD enrichment pulls down methylated DNA, so only genes with enough coverage appear here, and unmethylated genes are probably under-represented relative to the genome. ```{r gene-body, fig.width = 6, fig.height = 4} genes_bed <- fread(file.path(ann_dir, "genes.bed"), header = FALSE, col.names = c("chr", "start", "end", "gene", "score", "strand")) genes_info <- fread(file.path(ann_dir, "genes-info.tsv"), header = FALSE, col.names = c("gene", "GeneID", "biotype")) genes_product <- fread(file.path(ann_dir, "genes-product.tsv"), header = FALSE, col.names = c("gene", "product"), quote = "", sep = "\t") genes_product[, product := utils::URLdecode(product), by = gene] # GFF escapes commas etc. as %2C genes <- Reduce(function(x, y) merge(x, y, by = "gene", all.x = TRUE), list(genes_bed[, .(gene, chr, start = start + 1, end, strand)], genes_info, genes_product)) genes_gr <- GRanges(genes$chr, IRanges(genes$start, genes$end), gene = genes$gene) hits <- findOverlaps(meth_gr, genes_gr) gene_meth <- data.table(gene = genes_gr$gene[subjectHits(hits)], mean_meth = meth$mean_meth[queryHits(hits)])[ , .(CpGs = .N, mean_meth = round(mean(mean_meth), 1)), by = gene] gene_meth <- merge(genes[, .(gene, GeneID, biotype, product)], gene_meth, by = "gene") fwrite(gene_meth, file.path(out_dir, "gene-body-methylation.csv")) gm <- gene_meth[biotype == "protein_coding" & CpGs >= 5] cat("Protein-coding genes with >=5 CpGs at >=5x:", nrow(gm), "of", genes[biotype == "protein_coding", .N], "\n") p <- ggplot(gm, aes(mean_meth)) + geom_histogram(binwidth = 2.5) + labs(x = "Mean gene-body methylation (%)", y = "Genes", title = "Gene-body methylation", subtitle = "Protein-coding genes with >=5 CpGs at >=5x") + theme_bw() ggsave(file.path(out_dir, "gene-body-methylation.png"), p, width = 6, height = 4, dpi = 150) p ``` # Where do the DMLs and DMRs fall? For each list, the share of hits in each feature, next to the background (everything tested). The key comparison is real vs. the 9 random splits: if the real split's shares fall inside the range of the random ones, its genomic distribution is no different from noise. ```{r hit-distribution} feature_shares <- function(gr) { a <- annotate(gr) s <- a[, .N, by = category][CJ(category = factor(categories, levels = categories)), on = "category"] s[is.na(N), N := 0L] rbind(s[, .(feature = as.character(category), pct = 100 * N / max(sum(N), 1))], data.table(feature = "TE", pct = 100 * mean(a$in_TE))) } # The real split's list is one of the 10 permutation files, made with identical settings compare_lists <- function(background_bed, prefix) { bg <- feature_shares(read_bed_gr(background_bed))[, .(feature, background_pct = pct)] lists <- perm_files(prefix) shares <- rbindlist(lapply(names(lists), function(split) { feature_shares(read_hits(lists[[split]]))[, split := split] })) real <- shares[split == real_split, .(feature, real_pct = pct)] rand <- shares[split != real_split, .(random_min = min(pct), random_max = max(pct), random_mean = mean(pct)), by = feature] out <- Reduce(function(x, y) merge(x, y, by = "feature"), list(bg, real, rand)) out[, real_outside_random := real_pct < random_min | real_pct > random_max] out[, (2:6) := lapply(.SD, round, 1), .SDcols = 2:6] out[match(c(categories, "TE"), feature)] } dml_dist <- compare_lists(file.path(in_dir, "tested-cpgs-5x.bed"), "dml-5x") dmr_dist <- compare_lists(file.path(in_dir, "tested-tiles-1kb.bed"), "dmr-1kb") fwrite(dml_dist, file.path(out_dir, "dml-5x-feature-distribution.csv")) fwrite(dmr_dist, file.path(out_dir, "dmr-1kb-feature-distribution.csv")) dml_dist dmr_dist ``` Tiles (1 kb) often span more than one feature; each is assigned its highest-priority overlap. # Genes with DMLs or DMRs (real split) Genes whose body or promoter overlaps a DML or DMR from the real split. This is a list for browsing, **not a set of oil-responsive genes**: the random splits produce lists of the same size. ```{r hit-genes} promoters_bed <- fread(file.path(ann_dir, "promoters.bed"), header = FALSE, col.names = c("chr", "start", "end", "gene", "score", "strand")) regions <- rbind(genes[, .(gene, chr, start, end, part = "gene body")], promoters_bed[, .(gene, chr, start = start + 1, end, part = "promoter")]) regions_gr <- GRanges(regions$chr, IRanges(regions$start, regions$end)) hits_per_gene <- function(bed, label) { b <- fread(bed, header = FALSE, col.names = c("chr", "start", "end", "name")) if (nrow(b) == 0) return(data.table()) gr <- GRanges(b$chr, IRanges(b$start + 1, b$end)) diff <- as.numeric(sub("meth.diff=([-0-9.]+);.*", "\\1", b$name)) ov <- findOverlaps(gr, regions_gr) data.table(gene = regions$gene[subjectHits(ov)], part = regions$part[subjectHits(ov)], meth_diff = diff[queryHits(ov)], type = label) } hit_genes <- rbind(hits_per_gene(file.path(in_dir, "dml-5x.bed"), "DML"), hits_per_gene(file.path(in_dir, "dmr-1kb.bed"), "DMR")) gene_hits <- hit_genes[, .(DMLs = sum(type == "DML"), DMRs = sum(type == "DMR"), parts = paste(sort(unique(part)), collapse = ", "), mean_meth_diff = round(mean(meth_diff), 1)), by = gene] gene_hits <- merge(genes[, .(gene, GeneID, biotype, product, chr, start, end)], gene_hits, by = "gene") setorder(gene_hits, -DMRs, -DMLs) fwrite(gene_hits, file.path(out_dir, "genes-with-dml-dmr.csv")) cat("Genes overlapping a real-split DML or DMR:", nrow(gene_hits), "\n") head(gene_hits[, .(gene, product, DMLs, DMRs, parts, mean_meth_diff)], 20) ``` # GO sanity check A simple over-representation test: for each GO term, are genes containing a DML (gene body or promoter) more often annotated with it than genes with any tested CpG? Hypergeometric test, BH-adjusted, terms with 10–500 background genes. Run for the real split and the 9 random splits. If the real split gives no more significant terms than random ones, there is no functional signal to report. ```{r go-check} gaf <- fread(cmd = "zcat ../data/genome/GCF_002022765.2_C_virginica-3.0_gene_ontology.gaf.gz | grep -v '^!'", header = FALSE, sep = "\t", select = c(2, 5, 9), col.names = c("GeneID", "GO", "aspect"), quote = "") gaf <- unique(gaf) gaf[, GeneID := as.character(GeneID)] genes[, GeneID := as.character(GeneID)] genes_with <- function(bed) { gr <- read_hits(bed) if (length(gr) == 0) return(character()) unique(regions$gene[subjectHits(findOverlaps(gr, regions_gr))]) } background_genes <- genes_with(file.path(in_dir, "tested-cpgs-5x.bed")) bg <- merge(gaf, genes[gene %in% background_genes, .(gene, GeneID)], by = "GeneID") term_sizes <- bg[, .(bg_genes = uniqueN(gene)), by = .(GO, aspect)][bg_genes >= 10 & bg_genes <= 500] bg <- bg[GO %in% term_sizes$GO] n_bg <- uniqueN(bg$gene) go_test <- function(hit_genes) { hits <- intersect(hit_genes, bg$gene) counts <- bg[gene %in% hits, .(hit_genes = uniqueN(gene)), by = GO] res <- merge(term_sizes, counts, by = "GO", all.x = TRUE) res[is.na(hit_genes), hit_genes := 0L] res[, p := phyper(hit_genes - 1, bg_genes, n_bg - bg_genes, length(hits), lower.tail = FALSE)] res[, fdr := p.adjust(p, "BH")] res[order(p)] } go_by_split <- lapply(perm_files("dml-5x"), function(bed) go_test(genes_with(bed))) go_summary <- rbindlist(lapply(names(go_by_split), function(split) { r <- go_by_split[[split]] data.table(split = split, real_split = split == real_split, terms_tested = nrow(r), terms_fdr_05 = sum(r$fdr < 0.05), min_fdr = signif(min(r$fdr), 3)) })) setorder(go_summary, -terms_fdr_05, min_fdr) fwrite(go_summary, file.path(out_dir, "go-check-by-split.csv")) fwrite(go_by_split[[real_split]], file.path(out_dir, "go-check-real-split.csv")) cat("Background genes with GO terms:", n_bg, "; terms tested:", nrow(term_sizes), "\n") go_summary ``` # Summary ```{r session} sessionInfo() ``` What to take from this step: - **Landscape** (`landscape-by-feature.csv`, `gene-body-methylation.png`): where MBD-BS coverage and methylation sit. This holds regardless of the oil result and is the most reportable output. - **DML/DMR distribution:** only meaningful if `real_outside_random` is TRUE for some feature. - **Gene list and GO check:** only meaningful if the real split stands out from the random splits in `go-check-by-split.csv`.