In protViz/prora: ORA, sigORA and GSEA functionalities for proteomic data

knitr::opts_chunk$set(
  echo = FALSE,
  message = FALSE,
  warning = FALSE
)

if(!exists("progress")){
  progress <- function(howmuch, detail){
    invisible(NULL)
  }
}

GSEA <- eval(params$GSEA)

library(tidyverse)
library(knitr)

Introduction

The following analysis compares the enrichment of particular gene/protein set members towards the upper and lower end of the provided ranked protein list (e.g. ranked by fold changes, P-values, henceforth denoted generally as "score"). This analysis is commonly referred to as Gene Set Enrichment Analysis and a more detailed description of the method can be found in @Subramanian15545. In principle, the protein list is ranked by the provided scores and an enrichment score is calculated based on the relative positions the members of a particular gene set take in the whole list. To calculate a P-value and a corresponding FDR, adjusted for multiplicity [@BH1995], a permutation test approach is used. The default number of permutations WebGestaltR uses is $n_\text{perm}=1000$ \footnote{Note that due to run time issues the number of permutations is 10 per default in the fgczgseaora run scripts. If more accurate results are required nperm has to be adjusted accordingly in the run script.}.

Parameters

Organism: r GSEA$organism
Target Database: r GSEA$target
Contrast: r GSEA$fpath
Number of permutations: r GSEA$nperm
The rankTable.rnk table used.{target="_blank"}

`WebGestaltR` Results

repfile <- list.files(pattern = "Report_")

Since WebGestaltR provides its own summary report, please refer to WebGestalt result file{target="_blank"}.

GSEA Results {.tabset .tabset-pills}

Enriched Pathways

genesets <- GSEA$merged_data %>% 
  distinct(geneSet, .keep_all = TRUE) %>% 
  dplyr::mutate(set = paste0("<a href='",link,"'>",geneSet,"</a>")) %>% 
  dplyr::select(set, description, enrichmentScore, normalizedEnrichmentScore, pValue, FDR)

genesets <- genesets %>% arrange(normalizedEnrichmentScore)
DT::datatable(genesets,
              escape = FALSE,
              colnames = c("Pathway", "Description", "enrichmentScore", "normalizedEnrichmentScore", "P-value","Adj. P-value"),
              style = "bootstrap") %>% 
  DT::formatRound(digits = 3, columns = c("enrichmentScore", "normalizedEnrichmentScore", "pValue", "FDR"))

Input Data

colnames(GSEA$input_data) <- c("ID", "Score")
GSEA$input_data %>%
  dplyr::mutate(`Score` = round(`Score`, 2)) %>% 
  dplyr::arrange(Score) %>% 
  DT::datatable(colnames = c("ID", "Score"), width = 700, style = "bootstrap")

Visualisation

gettop <- function(data, n=10){
  data %>% 
    distinct(geneSet, .keep_all = TRUE) %>%
    top_n(n = n, wt = abs(normalizedEnrichmentScore)) %>% 
    dplyr::select(geneSet) %>% 
    inner_join(., data, by = "geneSet") -> out
  return(out)
}

out <- gettop(GSEA$merged_data)

out %>% 
  dplyr::select(userId, geneSet, score) %>% 
  spread(key = geneSet, value = score) -> out2
fields::image.plot(t(as.matrix(out2[,-1])), axes = FALSE)
axis(side = 1, at = seq(0, 1, length.out = ncol(out2)-1), labels = colnames(out2)[-1], lwd = 0, lwd.ticks = 0, las = 2, cex.axis = 0.8)
axis(side = 2, at = seq(0, 1, length.out = nrow(out2)), labels = out2$userId, lwd = 0, lwd.ticks = 0, las = 1, cex.axis = 0.5)

toplot <- GSEA$merged_data %>%
  dplyr::select(geneSet, geneSymbol) %>%
  rowid_to_column("IDD")
toplot <- toplot %>%
  tidyr::spread(geneSet, geneSymbol) %>%
  dplyr::select(-IDD) 

toplot <- toplot %>%
  as.list() %>%
  lapply(na.omit)

if(length(toplot) <= 1) {
  message("UpSetR plot cannot be displayed. Only one pathway enriched.")
} else {
  UpSetR::fromList(toplot) %>%
    UpSetR::upset(mb.ratio = c(0.5, 0.5), order.by = "freq", nsets = 20)
}