---
title: "goatea vignette"
author: "Maurits Unkel"
package: goatea
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    toc_depth: 4
    number_sections: true
vignette: >
  %\VignetteIndexEntry{2. GOATEA R package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r options, include = FALSE, echo = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, 
  comment = "#>",
  
  fig.width = 12,
  fig.height = 6,
  out.width = "100%", 
  fig.align = "center"
)
```

# Introduction

Geneset Ordinal Association Test Enrichment Analysis (GOATEA) provides a Shiny interface with interactive visualizations and utility functions for performing and exploring automated gene set enrichment analysis using the 'goat' package. 

GOATEA is designed to support large-scale and user-friendly enrichment workflows across multiple gene lists and comparisons, with flexible plotting and output options. Visualizations pre-enrichment include interactive Volcano and UpSet (overlap) plots. Visualizations post-enrichment include interactive geneset split term dotplot, geneset hierarchical clustering treeplot, multi-genelist gene-effectsize heatmap, enrichment overview gene-geneset heatmap and bottom-up pathway-like STRING database of protein-protein-interactions network graph.

## Why use goatea? 

Some enrichment visualization packages are included in Bioconductor, for instance: enrichplot, vissE, CBNplot. 
goatea additionally provides interactive visualizations through an intuitive and customizable Shiny user interface and allows for multi-genelist comparisons. No bioinformatic expertise is needed. goatea includes commonly used visualizations and adds a novel enrichment overview heatmap and offers a bottom-up pathway-like protein-protein interaction graph with network statistics.

A number of enrichment analysis packages are already on Bioconductor: clusterProfiler, fgsea, DOSE and many more. 
In goatea, goat is implemented. goat is a novel method for efficient geneset enrichment analysis. It performs within seconds and results in more identified significant terms compared to other methods, see the reference. For researchers with coding expertise, goatea provides an automated analysis workflow for performing enrichment analysis with goat and obtaining overview tables and output figures to explore your data. 

## How goatea works 

Genelists are often obtained from transcriptomic and/or proteomic experiments. These genelist tables have to contain a 'gene' column with NCBI Entrez gene identifiers, a 'symbol' column with gene aliases (optional: can be mapped from their gene IDs), and 'effectsize' and 'pvalue' columns with measurement outcomes for a performed comparison. User defined thresholds for effect size and p-values are used to define gene significance. To get a first look and understanding of your genelist data, a volcano plot is provided. For multi-genelist data a UpSet (overlap) plot can additionally be generated to visualize significant genes between comparisons. 

Genesets can be obtained by downloading Gene Ontology organism specific genesets within goatea or from .gmt files, downloaded for instance from the Molecular Signatures Database. Genelist genes and genesets will be filtered and matched, then goat or an older version of gene set enrichment analysis can be performed. Extensive filtering options help you to explore only the terms and genes within your vision.

Post-enrichment visualizations help to identify genes and terms of interest or relevance for your experimental design. Identified and selected genes of interest are then used to specifically visualize your data by zooming in on specific biological or technical aspects. 

The steps above can be automated, an example workflow is described in the rest of this vignette.  

As intended for interactive and intuitive usage, run the goatea Shiny interface.

# Installation

goatea developer version is available through GitHub, goatea is also available on Bioconductor. 

```{r installation, eval = FALSE}
## GOATEA installation requires R (v4.5.0) and the lastest version of Rtools 
## Rtools is needed for package compilation, to download and install visit: 
# R: https://cran.r-project.org/mirrors.html
# Rtools: https://cran.r-project.org/bin/windows/Rtools/
## I recommend to run the installation and GUI startup code via Rstudio, to have a nice graphical environment: 
# Rstudio: https://posit.co/downloads/

## To install GOATEA from Bioconductor (v3.23) use:
if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager")
BiocManager::install("goatea")
# or via pak
if ( ! require("pak", quietly = TRUE)) install.packages('pak')
pak::pkg_install('goatea', dependencies = TRUE, upgrade = TRUE)

## goatea requires at least one of the following available organism genome wide annotation packages:
### format: organism (taxid): org.Xx.eg.dg
# Human (9606)--------: org.Hs.eg.db
# Mouse (10090)-------: org.Mm.eg.db
# Fruit Fly (7227)----: org.Dm.eg.db
# Rhesus monkey (9544): org.Mmu.eg.db
# Rat (10116)---------: org.Rn.eg.db
# Worm (6239)---------: org.Ce.eg.db
# Chimpanzee (9598)---: org.Pt.eg.db
# Zebrafish (7955)----: org.Dr.eg.db
if ( ! require("pak", quietly = TRUE)) install.packages('pak')
pak::pkg_install(c(
  "org.Hs.eg.db", 
  "org.Mm.eg.db", 
  "org.Dm.eg.db", 
  "org.Mmu.eg.db", 
  "org.Rn.eg.db", 
  "org.Ce.eg.db", 
  "org.Pt.eg.db", 
  "org.Dr.eg.db"
))
```

# Running GOATEA: Shiny application

For an application demonstration, see the [GUI vignette](goatea_GUI.html)

Easiest: run GOATEA in your browser via HuggingFace: https://huggingface.co/spaces/Mausaya/GOATEA
Note: this may be somewhat slower, as 16GB RAM and 2CPU are shared across all users. 




For local usage, simply install the package and its dependencies, then run the script below.

Note that the GOATEA color scheme is easily customizable, have fun creating your own theme! 

Click the '?' buttons and hover over UI elements for explanation in tooltips. 

```{r goatea_Shiny, eval=FALSE, message=FALSE, warning=FALSE}
library(goatea)

## customizable coloring
colors <- list(
  main_bg = "#222222",
  darker_bg = "#111111",
  focus = "#32CD32", 
  hover = "#228B22",
  border = "#555555",
  text = "#FFFFFF"
)

## set max file size to ...MB for uploads (default: 100)
options(shiny.maxRequestSize = 1024^2 * 100)

## run the goatea Shiny application
shiny::shinyApp(
  ui = goatea:::goatea_ui,
  server = function(input, output, session) {
    goatea:::goatea_server(
      input, output, session, 
      css_colors = colors)
  }
)
```

# Running goatea: automated analyses

## Initialization

Note that all plots in the vignette are clickable to full screen, in case the images load too small. 

```{r library, message=FALSE, warning=FALSE}
# load goatea package library
library(goatea)
```

Set parameter values used throughout goatea. 

```{r initialization}
# set chosen organism taxid (for genome wide annotation package: org.Hs.eg.db)
taxid <- 9606 # Human (Homo Sapiens)

# set and create output folder
outdir <- tempdir()
dir.create(outdir, recursive = TRUE)

# set significance thresholds
p_value_threshold <- 0.05
absolute_effectsize_threshold <- 1
```

## Loading genelists

Transcriptomic and/or proteomic experiments often end up in comparing groups through differential gene expression or something similar. There, lists of genes with measurement weights and statistical probability values are obtained. Genelists have to be in .csv/.tsv/.xlsx format. Genelist tables contain a 'gene' column with NCBI Entrez gene identifiers, a 'symbol' column with gene aliases (optional: can be mapped from their gene IDs), and 'effectsize' and 'pvalue' columns with measurement outcomes for a performed comparison. Here, goat manuscript example data will be used, and expected use through uploading your own file is shown. 

To upload multiple files at once, ensure they are located in the same folder. You can then select them simultaneously by holding the Ctrl (Windows) or Command (Mac) key while clicking individual files, or using the Shift key to select a continuous range.

```{r loading_genelists}
# GOAT paper example data (example_Colameo_data):
load(system.file("extdata", "example_Colameo.rda", package = "goatea"))
## to manually download (to a specific folder, tempdir() by default) use: 
# data <- goat::download_goat_manuscript_data(output_dir = outdir)

# We select datasets from Colameo 2021 for our automated workflow.
# This is RNA-seq with matched mass-spectrometry data. 
# Rat data are mapped to Human NCBI Entrez gene identifiers within the goat package.
genelists <- list(
  rna = example_Colameo_data$`Colameo 2021:RNA-seq:PMID34396684`,
  ms = example_Colameo_data$`Colameo 2021:mass-spec:PMID34396684`
)

head(genelists[['rna']])
```

To load your own data, you can follow this alternative example.

```{r loading_genelists_user_data}
# set the filepath to your own data files 
file <- system.file('extdata','example_genelist.csv', package = 'goatea')
# Check file formatting, see function documentation for additional parameters 
genelist <- goatea::read_validate_genelist(file = file)
# genelists <- list(
  # csv = genelist
# )
head(genelist)
```

## Loading genesets

symbolsets are terms of biological aspects associated with a set of symbols. The Gene Ontology (GO) and Molecular Signatures databases (MSigDB) are commonly used repositories containing organism specific symbolsets. goat provides GO symbolsets for download, also to supply a .gmt file containing symbolsets downloaded from MSigDB for instance. 

```{r loading_genesets}
genesets <- goat::load_genesets_go_bioconductor(taxid = taxid)

## optionally: save and load from .rda to speed up after downloading
# save(genesets, file = file.path(outdir, 'genesets.rda'))
# load(file.path(outdir, 'genesets.rda')) # loading .rda to variable: genesets

## alternatively: load via .gmt file (downloaded from e.g.: MSigDB)
# then use: goat::load_genesets_gmtfile()

head(genesets)
```

# GOAT Enrichment Analysis 

Gene significance is set per genelist by the user-defined thresholds. The genelist is then ordered by effect size. For this example, if there are more genes than the amount of precomputed nulldistributions goat uses for its enrichment testing the amount of genes is limited to that amount. If you do want to use all genes, set the method to 'goat_bootstrap'. 

Genelist genes and genesets will be filtered and matched, then goat or an older version of gene set enrichment analysis can be performed. Too small or too large genesets can be irrelevant or mess up testing statistics, leaving parameter values to defaults is a good starting point. 

Geneset enrichment with goat is by default performed on 'effectsize'. The user can also go for 'pvalue' or a specific direction by gene effectsize sign or absolute values. Besides the goat method, older methods are also implemented to enable easy comparison.

```{r enrichment_analysis}
filtered_genesets_list <- list()
enrichment_results <- list()
for (name in names(genelists)) {
  genelist <- genelists[[name]]
  ## set signif column by given p-value- and effectsize significance thresholds
  genelist$signif <- genelist$pvalue <= p_value_threshold & 
                     abs(genelist$effectsize) >= absolute_effectsize_threshold
  # order genelist by effectsize 
  genelist <- genelist[order(abs(genelist$effectsize), decreasing = TRUE),]
  ## keep max n genes for regular goat method
  # optionally: remove this line and use goat_bootstrap method 
  if (nrow(genelist) > max(goat::goat_nulldistributions$N)) {
    genelist <- genelist[1:max(goat::goat_nulldistributions$N),]
  }
  # update genelist
  genelists[[name]] <- genelist
  
  ## filter genesets 
  filtered_genesets <- goat::filter_genesets(
    genesets = genesets,
    genelist = genelist,
    min_overlap = 10L,
    max_overlap = NA,
    max_overlap_fraction = 0.5,
    min_signif = NA,
    max_size = NA,
    dedupe = FALSE
  )
  filtered_genesets_list[[name]] <- filtered_genesets

  ## run enrichment 
  enrichment_results[[name]] <- goatea::run_geneset_enrichment(
    genesets = filtered_genesets,
    genelist = genelist,
    method = "goat",
    score_type = "effectsize",
    padj_method = "BH",
    padj_sources = TRUE,
    padj_cutoff = p_value_threshold,
    padj_min_signifgenes = 0L
  )
}
head(enrichment_results[[1]])
```

## Genes overview table

Here, the genelists are formatted into a gene-based overview table. If there are multiple genelists, for all genelists, each gene and their p-values, effect sizes, effect size percentages, gene set ratios and overlapping significant genes will be shown. 

```{r genes_overview}
genes_overview <- goatea::run_genelists_overlap(genelists)
genes_overview <- goatea::calculate_geneSetRatio(
  enrichment_results, 
  genes_overview
)
head(genes_overview)
```

## Enrichment searching and filtering

You can filter enrichment by querying the term names and by thresholding the numerical features. In this way you can explore relevant biological aspects and/or filter out technicalities. After filtering, use the filtered enrichment as input with the visualizations for zooming in on specific terms and their genes and interactions. 

```{r searching_and_filtering}
filtered_enrichment <- goatea::filter_enrichment(
  df = enrichment_results[[1]],
  terms_query = c("synapse", "synaptic"), 
  terms_antiquery = c("asymmetric"),
  min_ngenes = 100, 
  min_ngenes_signif = 10
)
head(filtered_enrichment$name)
```

## Gene selection

In the goatea app, the visualizations are used to select specific genes to pass to next visualizations in order to zoom in on specific terms or groups of genes. Here we simply take the significant genes of the first filtered enrichment to show some visualizations that use gene selections. 

```{r selecing_genes}
selected_genes <- filtered_enrichment$genes_signif[[1]]
```

# Plotting

## pre-enrichment: Volcano plot

A Volcano plot is commonly used to have a visual representation of the genelist. It shows an overview of genes as points with effectsize vs p-value and significance thresholds.  
```{r genelist_volcano}
goatea::plot_EnhancedVolcano(
  genelist = genelists[[1]],
  effectsize_threshold = absolute_effectsize_threshold,
  pvalue_threshold = p_value_threshold,
  background_color = 'white',
  foreground_color = 'black'
)
```

## pre-enrichment: UpSet overlap plot

If multiple genelists are loaded, their overlap of significant genes is shown through an UpSet plot. These show the total number of genes in a genelist and the number of overlapping significant genes for each intersection.

UpSet plot offers the intersect, distinct and union modes for overlapping and then selecting genes, which are taken along for downstream heatmap visualizations and PPI network plotting. Intersect are genes present in all selected contrasts, to visualize a shared core. Distinct genes are unique to specific intersections, to visualize exclusive signatures. Union genes are all genes found in any of the selected contrasts, to visualize a comprehensive overview. 

```{r genelist_overlap, fig.width = 8, fig.height = 4}
if (length(genelists) > 1) {
  goatea::plot_genelists_overlap_upsetjs(
    genelists = genelists, 
    mode = 'distinct'
  )
}
``` 

## Post-enrichment: termtree plot 
The tree plot shows a clustered hierarchical split between terms based on their term similarity. Similarity is based on semantics of the term. Terms are clustered together by semantics, meaning words that look or mean alike. A summary of the clustered semantics is given to the right of the vertical lines for each cluster. 

NOTE: As the 'enrichplot' package used internally no longer allows for GOAT enrichment to be converted to the class needed for plotting their 'treeplot', I rerun a local over-representation analysis on all genes of the selected genelist versus the selected (filtered) genesets of the current enrichment view.

```{r tree_plot}
goatea::plot_termtree(
  genelist = genelists[["rna"]], 
  genesets = filtered_genesets_list[["rna"]], 
  map_organism = taxid, 
  effectsize_threshold = absolute_effectsize_threshold, 
  Nterms = 20, 
  Nwords = 5, 
  Nclusters = 3
)
```

## Post-enrichment: split-dot plot

The dotplot shows terms split by up- and down regulation. The regulation is defined by the sign of the effect sizes of the genes in the genelists. The gene ratio shows the percentage of genes of the term that were matched with the genes of the genelist. 

```{r term_split_dotplot, fig.width = 16, fig.height = 9}
# filter to keep a single geneset source
single_source_enrichment <- enrichment_results[[1]] %>% 
  dplyr::filter(source == enrichment_results[[1]]$source[1])
goatea::plot_splitdot(single_source_enrichment, topN = 20)
```

## Post-enrichment: gene-effectsize heatmap

This heatmap shows the selected genes for each genelist to compare their effect sizes and p-values as a simple overview. 

```{r gene_genelist_heatmap, fig.width = 12, fig.height = 8}
selected_symbols <- filtered_enrichment$symbol[[1]][
  filtered_enrichment$genes[[1]] %in% filtered_enrichment$genes_signif[[1]]
][1:15]
p <- goatea::plot_gene_effectsize_ComplexHeatmap(
  genes = selected_symbols, 
  genes_overview = genes_overview
)
```

## Post-enrichment: gene-genesets heatmap

This heatmap serves as an overview of the enrichment result, combined with all genelists information. Use the search functionality and pass the resulting filtered enrichment to zoom in on relevant biological aspects and focus on filtered aspects of the enrichment for exploration. In the goatea app this heatmap is also used to select genes/terms from identified aspects of interest. 

```{r gene_term_heatmap}
p <- goatea::plot_ComplexHeatmap(
  enrichment_result = enrichment_results[[1]],
  genelist = genelists[[1]],
  n_cluster = 5,
  n_top_terms = 20,
  n_top_genes = 50,
  genelist_overlap = genes_overview
)
```

## Post-enrichment: protein-protein interaction network graph

Gene sets miss information of gene-gene interactions, this visualization can offer bottom-up relationships between genes to gain insight. The STRING database is used to obtain protein-protein interaction data by efficient streaming. For the selected genes of interest, an igraph network is built in the background in order to obtain network statistics. The igraph is passed to visNetwork for (in Shiny app) interactive visualization and exploration of the network. 

First we select symbols to download ppi data for our chosen organism. 

```{r ppi_network_select_symbols}
# selecting 10 symbols of the searchad and significant first enriched term
selected_symbols <- filtered_enrichment$symbol[[1]][
  filtered_enrichment$genes[[1]] %in% filtered_enrichment$genes_signif[[1]]
][1:10]
```

```{r ppi_network_example_default}
assign("ppi_data", get(load(system.file("extdata", "example_ppi_data.rda", package = "goatea"))))
```

Normally, we would download STRING database PPI data, yet for this vignette, we will use example data by default. 
Try out this code for yourself to explore with downloaded data. 

```{r ppi_network_example_download, eval = FALSE}
## this will download the STRING database protein-protein interaction data 
## to your selected folder 
# downloading should only take a couple of seconds
# downlaod happens only if the files are not found in the specified folder
# if the download does not work, double check if the STRING database is online
ppi_data <- goatea::get_string_ppi(
  aliases = selected_symbols, 
  organism = taxid, 
  folder = file.path(outdir)
)
```

I advise to regard this output as a basic and toy visualization, the goatea in Shiny app version contains node and edge hover information with all network statistics. Edges can be clicked to view their interaction scores and information on the STRING website. Nodes are clustered by STRING interaction scores. The nodes and edges can be colored by each statistical feature. Nodes can be filtered by each feature and made into a subnetwork for more detailed exploration. 

```{r ppi_network}
# creating igraph in order to do network statistics
ppi_graph <- goatea::get_ppigraph(ppi_data = ppi_data)
# converting igraph to visnetwork for (Shiny app) interactive visualization
vis_network <- goatea::get_visNetwork(
  ppigraph = ppi_graph, 
  genes_overview = genes_overview
)
## showing a basic version of the visNetwork
## the goatea Shiny app version has full interactive capability 
visNetwork::visNetwork(
  vis_network$nodes, 
  vis_network$edges, 
  width = "100%", 
  height = "100%") %>% 
  visNetwork::visPhysics(stabilization = FALSE)
```

# Information

A thank you for your time and effort in using goatea, I hope it may aid you in exploring your data! 

## Contact 

goatea was developed by Maurits Unkel at the Erasmus Medical Center in the department of Psychiatry.
Source code is available through GitHub: https://github.com/mauritsunkel/goatea
Questions regarding goatea can be sent to: mauritsunkel@gmail.com

## License

goatea is licensed under Apache 2.0.

## Issues and contributions

For issues and contributions, please find the goatea Github page: https://github.com/mauritsunkel/goatea/issues

## GOAT reference

Koopmans, F. GOAT: efficient and robust identification of gene set enrichment. Commun Biol 7, 744 (2024). https://doi.org/10.1038/s42003-024-06454-5

## Session info

```{r session_info}
sessionInfo()
```