GSVA.qmd

---
title: "GSVA"
author: "Harvard Chan Bioinformatics Core"
date: "`r Sys.Date()`"
format:
  html:
    code-fold: true
    code-tools: true
    code-overflow: wrap
    df-print: paged
    highlight-style: pygments
    number-sections: true
    self-contained: true
    theme: default
    toc: true
    toc-location: right
    toc-expand: false
    lightbox: true
params:
  # set column name and contrasts to be factors of interest
  column: "sample_type"
  contrasts:
    value:
      - ["sample_type", "tumor", "normal"]
  project_file: ../information.R
  params_file: ../00_params/params-example.R # example data
  functions_file: ../00_libs
  # select from gene set repository at https://github.com/bcbio/resources/tree/main/gene_sets/gene_sets
  # choose geneset, click "Raw', and copy url to work with mouse data
  geneset_fn: https://raw.githubusercontent.com/bcbio/resources/main/gene_sets/gene_sets/20240904/human/h.all.v2024.1.Hs.entrez.gmt
---

```{r}
#| message: FALSE
#| warning: FALSE
# This set up the working directory to this file so all files can be found
library(rstudioapi)
library(tidyverse)
# NOTE: This code will check version, this is our recommendation, it may work
# .      other versions
stopifnot(R.version$major >= 4) # requires R4
if (compareVersion(R.version$minor, "3.1") < 0) warning("We recommend >= R4.3.1")
stopifnot(compareVersion(as.character(BiocManager::version()), "3.18") >= 0)
```

This code is in this ![](https://img.shields.io/badge/status-draft-grey) revision.

```{r libraries}
#| message: FALSE
#| warning: FALSE
## load libraries
library(GSVA)
library(GSEABase)
library(reshape2)
library(ChIPpeakAnno)
library(org.Hs.eg.db)
# library(org.Mm.eg.db)
library(AnnotationDbi)
library(DESeq2)
library(limma)
library(gridExtra)
library(bcbioR)
library(ggprism)
library(knitr)
library(clusterProfiler)

colors <- cb_friendly_cols(1:15)
ggplot2::theme_set(theme_prism(base_size = 14))
opts_chunk[["set"]](
  cache = F,
  cache.lazy = FALSE,
  dev = c("png", "pdf"),
  error = TRUE,
  highlight = TRUE,
  message = FALSE,
  prompt = FALSE,
  tidy = FALSE,
  warning = FALSE,
  echo = T,
  fig.height = 4)

# set seed for reproducibility
set.seed(1234567890L)
```


```{r load_params}
#| cache: FALSE
#| message: FALSE
#| warning: FALSE
source(params$project_file)
source(params$params_file)
purr::map(list.files(params$functions_file, pattern = "*.R$", full.names = T), source) %>% invisible()
column <- params$column
contrasts <- params$contrasts
subset_column <- params$subset_column
subset_value <- params$subset_value
```

```{r sanitize_datatable}
sanitize_datatable <- function(df, ...) {
  # remove dashes which cause wrapping
  DT::datatable(df, ...,
    rownames = gsub("-", "_", rownames(df)),
    colnames = gsub("-", "_", colnames(df))
  )
}
```

# Overview

-   Project: `r project`
-   PI: `r PI`
-   Analyst: `r analyst`
-   Experiment: `r experiment`
-   Aim: `r aim`

```{r load_data}
coldata <- load_coldata(
  coldata_fn, column,
  subset_column, subset_value
)
coldata[[contrasts[[1]][1]]] <- relevel(as.factor(coldata[[contrasts[[1]][1]]]), contrasts[[1]][3])
coldata$sample <- row.names(coldata)

counts <- load_counts(counts_fn)
counts <- counts[, coldata$sample]
```

# Method

Gene Set Variation Analysis (GSVA) is a non-parametric, unsupervised method for estimating variation of gene set enrichment through the samples of a expression data set. GSVA performs a change in coordinate systems, transforming the data from a gene by sample matrix to a gene-set by sample matrix, thereby allowing the evaluation of pathway enrichment for each sample. This new matrix of GSVA enrichment scores facilitates applying standard analytical methods like functional enrichment, survival analysis, clustering, CNV-pathway analysis or cross-tissue pathway analysis, in a pathway-centric manner. More info in the vignette [here](https://bioconductor.org/packages/release/bioc/vignettes/GSVA/inst/doc/GSVA.html).

Hänzelmann S, Castelo R, Guinney J (2013). “GSVA: gene set variation analysis for microarray and RNA-Seq data.” BMC Bioinformatics, 14, 7. doi:10.1186/1471-2105-14-7, [https://doi.org/10.1186/1471-2105-14-7](https://doi.org/10.1186/1471-2105-14-7)

# Data

```{r show_coldata}
coldata %>% sanitize_datatable()
```

```{r normalize_data}
dds <- DESeqDataSetFromMatrix(counts,
  colData = coldata,
  design = ~1
)

dds <- DESeq(dds)
norm_counts <- counts(dds, normalized = TRUE)
```


```{r ensembl_to_entrez}
## convert ensembl to entrez

entrezIDs_all <- convert2EntrezID(
  IDs = rownames(norm_counts), orgAnn = "org.Hs.eg.db",
  ID_type = "ensembl_gene_id"
)

entrezid <- mapIds(org.Hs.eg.db, keys = rownames(norm_counts), keytype = "ENSEMBL", column = "ENTREZID")

counts_entrez <- norm_counts
stopifnot(nrow(counts_entrez) == length(entrezid))
rownames(counts_entrez) <- entrezid
library(matrixStats)
if (any(duplicated(rownames(counts_entrez)))) {
  counts_entrez <- rowsum(as.matrix(counts_entrez), group = rownames(counts_entrez))
}
```


# Prep and run GSVA

```{r load_genesets}
gene_sets <- read.gmt(params$geneset_fn)

genes_by_pathway <- split(gene_sets$gene, gene_sets$term)
```


```{r GSVA}
#| message: FALSE
#| warning: FALSE
gsvaPar <- GSVA::gsvaParam(counts_entrez, genes_by_pathway, kcdf = "Poisson")

gsva.es <- gsva(gsvaPar, verbose = F)
```

## Test for Significance

```{r limma}
mod <- model.matrix(~ factor(coldata[[column]]))
fit <- lmFit(gsva.es, mod)
fit <- eBayes(fit)
res <- topTable(fit, coef = paste0("factor(coldata[[column]])", contrasts[[1]][2]), number = Inf, sort.by = "P")

res %>% sanitize_datatable()
```

## Graph top 5 pathways

```{r graph_pathways}
#| results: 'asis'

scores <- t(gsva.es)

sig <- subset(res, res$adj.P.Val < 0.1)

if (nrow(sig) >= 5) {
  pathways <- rownames(sig)[1:5]
} else if (nrow(sig) == 0) {
  pathways <- c()
} else {
  pathways <- rownames(sig)
}

if (length(pathways) > 0) {
  to_graph <- data.frame(scores[, pathways]) %>%
    rownames_to_column("sample") %>%
    pivot_longer(!sample, names_to = "pathway", values_to = "enrichment_score")
  to_graph <- left_join(to_graph, coldata)

  for (single_pathway in pathways) {
    cat("### ", single_pathway, "\n")

    to_graph_single_pathway <- to_graph %>% filter(pathway == single_pathway)
    p <- ggplot(to_graph_single_pathway, aes(x = .data[[column]], y = enrichment_score)) +
      geom_boxplot() +
      geom_point(alpha = 0.5) +
      ggtitle(single_pathway)
    print(p)

    cat("\n\n")
  }
} else {
  cat("No pathways were detected as significantly enriched")
}
```



# R session

List and version of tools used for the QC report generation.

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