R/BatchGSEA.R

In Rvirgenslane/Hotgenes: For Visualizing Differential Expression Analysis

#' Performs GSEA for all Hotgenes contrasts
#' @param HotgenesObj R object generated by DEseq2_export function
#' @param nTop number of pathways to return
#' @param m_df GeneSet generated using msigdbr function
#' @param cut_padj Adjusted pvalue cut off. Default is 0.1.
#' @inheritParams fgsea::fgsea
#' @importFrom stats reorder
#' @return HotgenesObj appended with GSEA for each contrast
#' @export
#' @examples Example_Hotgenes_dir<-system.file("extdata",
#' "Example_Hotgenes.Rdata",
#' package = "Hotgenes", mustWork = TRUE)
#' load(Example_Hotgenes_dir)
#' library(msigdbr)
#' 
#' # GO annotations
#' m_df = msigdbr(species = "Homo sapiens", category = "C5", subcategory = "BP")
#' 
#' # Reactome annotations
#' # m_df = msigdbr(species = "Homo sapiens", category = "C2", 
#' # subcategory = "CP:REACTOME")  
#' 
#' qbat<-BatchGSEA(HotgenesObj=Example_Hotgenes,
#' m_df= m_df)
#' 
#' # View enriched pathways
#' lapply(qbat$OuputGSEA, function(x) head(x$fgRes$pathway))
#'
#' # plot of enriched pathways
#' qbat$OuputGSEA$shEWS$g
#' 
#' # Prep for GSEA plot
#' m_list = qbat$m_list
#' 
#' # top pathways
#' qbat$OuputGSEA$shEWS$top$pathway
#' pyid<-qbat$OuputGSEA$shEWS$top$pathway[[2]]
#' pyid
#' 
#' Gene_Ranks <- qbat$OuputGSEA$shEWS$Gene_Ranks
#' fgsea::plotEnrichment(m_list[[pyid]], Gene_Ranks,
#' gseaParam = 1, ticksSize = 0.2) +
#' ggplot2::labs(title=stringr::str_wrap(pyid, 20))


BatchGSEA<-function(m_df=NULL, HotgenesObj=NULL,
nTop=10, cut_padj=0.1,
minSize=15,
maxSize=1000,
nperm=100000,
nproc=1){

m_df$gs_name<-gsub("_", " ", m_df$gs_name) # fix names
m_list = split(x = m_df$gene_symbol, f = m_df$gs_name)

contrast_name<-names(HotgenesObj$Output_DE)
contrast_name_lists <- setNames(vector(length(contrast_name),
mode="list"), contrast_name)


for(DE_sel in contrast_name){

clean_table<-HotgenesObj$Output_DE[[DE_sel]]
ranks<-clean_table[order(clean_table$stat, decreasing = TRUE),]
Gene_Ranks <- ranks$stat
names(Gene_Ranks) <- rownames(ranks)

# fgsea of pathways -------------------------------------------------------
fgseaRes <- fgsea::fgsea(pathways = m_list, 
stats = Gene_Ranks,
minSize=minSize,
maxSize=maxSize,
nperm=nperm,
nproc=nproc)

fgRes<-fgseaRes[fgseaRes$padj<cut_padj,]

# Filtering
Up <- head(fgRes[order(fgRes$NES, decreasing = TRUE)], n=nTop)
Down <- head(fgRes[order(fgRes$NES, decreasing = FALSE)], n = nTop)

top<-rbind(Up,Down)

top$Enrichment = ifelse(top$NES > 0, "Up-regulated", "Down-regulated")
top<-unique(data.frame(top, stringsAsFactors = FALSE))

filtRes = top
filtRes$pathway<-stringr::str_wrap(filtRes$pathway, 50)

# plot
g <- ggplot(filtRes, aes(reorder(.data$pathway, .data$NES), .data$NES)) +
geom_segment(aes(reorder(.data$pathway, .data$NES),
xend=.data$pathway, y=0, yend=.data$NES)) +
geom_point(size=3, aes(fill = .data$Enrichment),
shape=21, stroke=0.6) +
scale_fill_manual(values = c("Down-regulated" = "dodgerblue",
"Up-regulated" = "firebrick") ) +
coord_flip() +
labs(title= DE_sel,x="Pathway", y="Normalized Enrichment Score") +
theme_minimal(base_size = 8) +
theme(axis.text  = element_text(size=8))


GSEA_OUT<-list(fgRes=fgRes,
top=top,
g=g,
Gene_Ranks=Gene_Ranks)

contrast_name_lists[[DE_sel]]<-GSEA_OUT

}

HotgenesObj$OuputGSEA<-contrast_name_lists
HotgenesObj$m_list<-m_list
return(HotgenesObj)
}