# Module 7: from gene lists to pathways.  Video: T9.
# Working directory: the "workshop" folder. Needs results/deseq_1D.txt from module 5.
suppressPackageStartupMessages({ library(piano); library(fgsea) })
set.seed(1)
dir.create("results", showWarnings = FALSE)

de <- read.delim("results/deseq_1D.txt", row.names = 1)
rownames(de) <- toupper(rownames(de))     # the gene sets use upper-case symbols
kegg <- loadGSC("data/KEGG.gmt")
sets <- kegg$gsc
length(sets)

# ---- A. over-representation by hand, for one pathway -----------------------------------------------
universe <- rownames(de)[!is.na(de$padj)]                       # genes that were actually tested
sig <- rownames(de)[!is.na(de$padj) & de$padj < 0.05]
ora <- function(pathway) {
  s <- intersect(sets[[pathway]], universe)
  k <- length(intersect(s, sig)); K <- length(s); n <- length(sig); N <- length(universe)
  c(K = K, observed = k, expected = K * n / N, p = phyper(k - 1, K, N - K, n, lower.tail = FALSE))
}
ora("Oxidative phosphorylation")
ora("HIF-1 signaling pathway")
# the same test as a 2 x 2 table:
s <- intersect(sets[["Oxidative phosphorylation"]], universe)
tab <- table(in_pathway = universe %in% s, significant = universe %in% sig)
tab
fisher.test(tab, alternative = "greater")$p.value
# EXERCISE A1: run ora() for every pathway, adjust with p.adjust(..., "BH") and list the top 10.
#              Then split `sig` into up and down genes and repeat. What changes?

# ---- B. piano (as in the original workshop) ---------------------------------------------------------------
lfc <- as.matrix(de$log2FoldChange); pval <- as.matrix(de$pvalue)
rownames(lfc) <- rownames(pval) <- rownames(de)
lfc[is.na(lfc)] <- 0; pval[is.na(pval)] <- 1
gsa <- runGSA(pval, lfc, gsc = kegg, geneSetStat = "reporter", signifMethod = "nullDist", nPerm = 1000)
res_piano <- GSAsummaryTable(gsa)
write.table(res_piano, "results/piano_kegg_1D.txt", sep = "\t", quote = FALSE, row.names = FALSE)
sum(res_piano$`p adj (dist.dir.up)` < 0.05, na.rm = TRUE)
sum(res_piano$`p adj (dist.dir.dn)` < 0.05, na.rm = TRUE)
pdf("results/piano_heatmap_1D.pdf", width = 8, height = 14)
GSAheatmap(gsa, adjusted = TRUE)
dev.off()
nw <- networkPlot2(gsa, class = "distinct", direction = "both", significance = 0.00005)
# a single self-contained HTML file needs pandoc (bundled with RStudio); otherwise a folder of files is written
visNetwork::visSave(nw, file = file.path(getwd(), "results/network_kegg_1D.html"), background = "white",
                    selfcontained = requireNamespace("rmarkdown", quietly = TRUE) && rmarkdown::pandoc_available())

# ---- C. fgsea: a ranked-list method -------------------------------------------------------------------------
stat <- de$stat; names(stat) <- rownames(de)
stat <- stat[!is.na(stat)]
fg <- fgsea(sets, stat, minSize = 10, maxSize = 500)
head(fg[order(fg$padj), c("pathway", "NES", "padj", "size")], 10)
plotEnrichment(sets[["Oxidative phosphorylation"]], stat) + ggplot2::ggtitle("Oxidative phosphorylation")

# ---- D. sanity check ---------------------------------------------------------------------------------------
# piano calls "Olfactory transduction" distinctly up in the heart. Look at the genes before believing it.
olf <- intersect(sets[["Olfactory transduction"]], rownames(de))
summary(de[olf, "baseMean"])
c(n = length(olf), filtered_or_untested = sum(is.na(de[olf, "padj"])))
# EXERCISE D1: does fgsea also call it? (Careful: the set has ~1,100 genes, more than maxSize = 500 above,
#              so fgsea silently skipped it. Rerun fgsea on sets["Olfactory transduction"] with maxSize = 2000.)
# EXERCISE D2: and ORA on the up genes only? If every method agrees, is the result therefore biology?
#              What else could make ~1,100 barely expressed genes look "up" in MI? (Hint: plotCounts on a few.)