---
title: "08 - Genomic-feature and GO annotation"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(eval = FALSE, message = FALSE, warning = FALSE)
library(rtracklayer); library(GenomicRanges); library(data.table)
```

# Reference annotation

Annotate against the NCBI RefSeq GFF for the assembly. Confirm the GFF and cov
files use the same seqid namespace (both RefSeq accessions, e.g. NC_047559.1).

```{r gff}
dir.create("data", showWarnings = FALSE)
gff_url <- paste0("https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/902/806/645/",
  "GCF_902806645.1_cgigas_uk_roslin_v1/",
  "GCF_902806645.1_cgigas_uk_roslin_v1_genomic.gff.gz")
if (!file.exists("data/cgigas_roslin.gff.gz"))
  download.file(gff_url, "data/cgigas_roslin.gff.gz", mode = "wb")

gff <- import("data/cgigas_roslin.gff.gz")
genes <- gff[gff$type == "gene"]
exons <- gff[gff$type == "exon"]
cds   <- gff[gff$type == "CDS"]
introns <- setdiff(reduce(genes), reduce(exons))
proms <- promoters(genes, upstream = 1000, downstream = 0)
features <- list(cds = cds, exons = exons, introns = introns, proms = proms)
```

# Hierarchical feature assignment

Priority CDS > exon > intron > promoter > intergenic.

```{r assign}
assign_feature <- function(gr, features) {
  feat <- rep("intergenic", length(gr))
  feat[overlapsAny(gr, features$proms)]            <- "promoter"
  feat[overlapsAny(gr, reduce(features$introns))]  <- "intron"
  feat[overlapsAny(gr, reduce(features$exons))]    <- "exon"
  feat[overlapsAny(gr, reduce(features$cds))]      <- "CDS"
  factor(feat, levels = c("CDS","exon","intron","promoter","intergenic"))
}
tocsv_gr <- function(dt) GRanges(dt$chr, IRanges(dt$start, dt$end))

bg  <- fread("diffmeth_interaction_all.csv.gz")           # all tested CpGs = background
sets <- list(
  background   = bg,
  ploidy_main  = fread("dml_ploidy_main.csv"),
  pH_within_2N = fread("dml_pH_within_2N.csv"),
  pH_within_3N = fread("dml_pH_within_3N.csv"),
  interaction  = fread("dml_interaction.csv"))
featdist <- rbindlist(lapply(names(sets), function(nm){
  f <- assign_feature(tocsv_gr(sets[[nm]]), features)
  data.table(contrast = nm, feature = names(table(f)),
             pct = 100 * as.numeric(table(f)) / length(f))
}))
fwrite(featdist, "dml_feature_distribution.csv")
```

Expected: DMLs strongly depleted in intergenic (~6-7% vs ~25% background) and
enriched in gene bodies (intron ~50-56%, CDS ~28-35%) - canonical invertebrate
gene-body methylation.

# DML -> gene -> GO, hypergeometric enrichment

Gene-to-GO mapping uses the project's `geneid2go` table (numeric NCBI GeneID
parsed from the GFF `Dbxref`, joined to transcripts then GO terms). The GO
universe is restricted to genes overlapping any tested CpG; terms need >= 5
genes; p-values BH-adjusted.

```{r go}
go_enrich <- function(sel_genes, bg_genes, gid2go, min_term = 5) {
  g <- gid2go[geneid %in% bg_genes]
  univ <- unique(g$geneid); Nuniv <- length(univ)
  term_sz <- g[, .(N = uniqueN(geneid)), by = go]
  sel <- intersect(sel_genes, univ); ng <- length(sel)
  hit <- g[geneid %in% sel, .(k = uniqueN(geneid)), by = go]
  tab <- merge(hit, term_sz, by = "go")[N >= min_term]
  tab[, pval := phyper(k - 1, N, Nuniv - N, ng, lower.tail = FALSE)]
  tab[, fold := k / (N * ng / Nuniv)]
  tab[, qval := p.adjust(pval, "BH")]
  setorder(tab, pval)[]
}
# gid2go: data.table(geneid, go) built from Dbxref + geneid2go.tab; then e.g.
# go_res <- go_enrich(genes_in_interaction_DMLs, background_genes, gid2go)
# fwrite(go_res, "go_enrichment_all.csv")
```

Interaction-DML genes are enriched for cytoskeleton/dynein/microtubule, TGF-beta
and Ras signalling, and transcriptional regulation.