params: functional_title: "Functional Analysis" title: "r params$functional_title" date: "r format(Sys.time(), '%d %B %Y %H:%M:%S', 'Europe/Vienna', T)" output: rmarkdown::html_document: self_contained: yes theme: null highlight: null mathjax: null


library(DOSE)
library(ggplot2)
library("clusterProfiler")
library("deseq2analysis")
options(knitr.table.format = "html") 
theme_set(theme_bw())
outtabledir <- paste(outbasefun,"_functional_analysis_tables",sep="")
dir.create(outtabledir)


deseq2parameterfile <- paste(outbasefun,"_deseq2parameters.RDS",sep="")
deseq2parameters <- NULL
if(file.exists(deseq2parameterfile)){
  deseq2parameters <- readRDS(deseq2parameterfile)
}

goup <- readRDS(paste(outbasefun,"_up.RDS",sep=""))
godn <- readRDS(paste(outbasefun,"_dn.RDS",sep=""))
gode <- readRDS(paste(outbasefun,"_de.RDS",sep=""))

showKegg <- !is.null(goup$kegg)
showReactome <- !is.null(goup$react)

GSEA_PATH <- paste(outbasefun,"_gse.RDS",sep="")
showGSEA <- file.exists(GSEA_PATH) 
gseaL <- NULL
if(showGSEA){
  gseaL <- readRDS(showGSEA)
}

All functional analyses were done with three sets of genes: upregulated genes, downregulated genes and all deregulated genes together. r goparameterdescription(deseq2parameters) The GeneOntology analysis test for overrepresentation of groups of genes within the GO categories. Because the gene-ontology is a directed acyclic graph - picture a branching tree - there is some redundancy in the terms and genes.

Upregulated Genes

GO

GO molecular function

   emptyTable(goup$mf, datatabGO, "UpGoMf", outtabledir, outbasefun)

GO biological process

   emptyTable(goup$bp, datatabGO, "UpGoBp", outtabledir, outbasefun)

GO cellular component

   emptyTable(goup$cc, datatabGO, "UpGoCc", outtabledir, outbasefun)

KEGG

   emptyTable(goup$kegg, datatabKEGG, "UpKEGG", outtabledir, outbasefun)

Reactome

   emptyTable(as.data.frame(goup$react), datatabReact, "UpReact", outtabledir, outbasefun)
   plotEmpty(goup$react,  enrichplot::dotplot(goup$react, showCategory=15))
   plotEmpty(goup$react,  enrichplot::emapplot(goup$react))

Downregulated Genes

GO

GO molecular function

   emptyTable(godn$mf, datatabGO, "DnGoMf", outtabledir, outbasefun)

GO biological process

   emptyTable(godn$bp, datatabGO, "DnGoBb", outtabledir, outbasefun)

GO cellular component

   emptyTable(godn$cc, datatabGO, "DnGoCc", outtabledir, outbasefun)

KEGG

   emptyTable(godn$kegg, datatabKEGG, "DnKEGG", outtabledir, outbasefun)

Reactome

   datatabReact(as.data.frame(godn$react), "DnReact", outtabledir, outbasefun)
   plotEmpty(godn$react,  enrichplot::dotplot(godn$react, showCategory=15))
   plotEmpty(godn$react,  enrichplot::emapplot(godn$react))

Deregulated Genes

GO

GO molecular function

   emptyTable(gode$mf, datatabGO, "DrGoMf", outtabledir, outbasefun)

GO biological process

   emptyTable(gode$bp, datatabGO, "DrGoBp", outtabledir, outbasefun)

GO cellular component

   emptyTable(gode$cc, datatabGO, "DrGoCc", outtabledir, outbasefun)

KEGG

   emptyTable(gode$kegg, datatabKEGG, "DrKEGG", outtabledir, outbasefun)

Reactome

   emptyTable(as.data.frame(gode$react), datatabReact, "DrReact", outtabledir, outbasefun)
   plotEmpty(gode$react,  enrichplot::dotplot(gode$react, showCategory=15))
   plotEmpty(gode$react,  enrichplot::emapplot(gode$react))

Gene Set Enrichment Analysis

For Gene Set Enrichment Analysis (GSEA) no arbitraty cut-offs are used. The whole list, ranked by log-fold change is compared to gene sets. Multiple Gene Set categories are extracted e.g. the Hallmark gene sets and for each gene sets the top 10 up and downregulated pathways are shown.

Hallmark gene sets

Hallmark gene sets are coherently expressed signatures derived by aggregating many MSigDB gene sets to represent well-defined biological states or processes.

   topPathways <- extractGSEAOverviewCollapsedTop(gseaL, 10, "H", gseaL$pvalGo) 
   fgsea::plotGseaTable(gseaL$pathways[topPathways], gseaL$stats, gseaL$gsea, gseaParam = 0.5)
   datatabGSEA(gseaL$gsea, topPathways, "GSEAH", outtabledir, outbasefun)       

Curated gene sets

Curated gene sets from online pathway databases, publications in PubMed, and knowledge of domain experts.

   topPathways <- extractGSEAOverviewCollapsedTop(gseaL, 10, gseaL$pvalGo, "C2")
   fgsea::plotGseaTable(gseaL$pathways[topPathways], gseaL$stats, gseaL$gsea, gseaParam = 0.5)
   datatabGSEA(gseaL$gsea, topPathways, "GSEAC2", outtabledir, outbasefun)

Go gene sets

GO gene sets consist of genes annotated by the same GO terms.

   topPathways <- extractGSEAOverviewCollapsedTop(gseaL, 10, gseaL$pvalGo, "C5")
   fgsea::plotGseaTable(gseaL$pathways[topPathways], gseaL$stats, gseaL$gsea, gseaParam = 0.5)
   datatabGSEA(gseaL$gsea, topPathways, "GSEAC5", outtabledir, outbasefun)

Motif gene sets

Motif gene sets based on conserved cis-regulatory motifs from a comparative analysis of the human, mouse, rat, and dog genomes.

   topPathways <- extractGSEAOverviewCollapsedTop(gseaL, 10, gseaL$pvalGo, "C3")
   fgsea::plotGseaTable(gseaL$pathways[topPathways], gseaL$stats, gseaL$gsea, gseaParam = 0.5)
   datatabGSEA(gseaL$gsea, topPathways, "GSEAC3", outtabledir, outbasefun)
   zipfunctabs(outtabledir, KNITDIR)

session info

   sessioninfo::session_info()


idot/deseq2analysis documentation built on Aug. 14, 2021, 5:12 a.m.