GSEA

From 太極
Revision as of 21:38, 18 April 2021 by Brb (talk | contribs)
Jump to navigation Jump to search

Basic

https://en.wikipedia.org/wiki/Gene_set_enrichment_analysis

Determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states

  • fgsea package and download stat 
    library(fgsea)
library(ggplot2)
data(examplePathways) # a list of length 1457 (pathways)
data(exampleRanks)    # a vector of length 12000 (genes)
set.seed(42)
plotEnrichment(examplePathways[["5991130_Programmed_Cell_Death"]],
               exampleRanks) + labs(title="Programmed Cell Death")

Two categories of GSEA procedures:

  • Competitive: compare genes in the test set relative to all other genes.
  • Self-contained: whether the gene-set is more DE than one were to expect under the null of no association between two phenotype conditions (without reference to other genes in the genome). For example the method by Jiang & Gentleman Bioinformatics 2007

See also BRB-ArrayTools -> GSEA.

Subramanian algorithm

In the plot, (x-axis) genes are sorted by their expression across all samples. Y-axis represents enrichment score. See HOW TO PERFORM GSEA - A tutorial on gene set enrichment analysis for RNA-seq. Bars represents genes being in the gene set. Genes on the LHS/RHS are more highly expressed in the experimental/control group. Small p means this gene set is enriched in this experimental sample.

ssGSEA

  • https://github.com/broadinstitute/ssGSEA2.0
  • ssGSEAProjection (v9.1.2). Each ssGSEA enrichment score represents the degree to which the genes in a particular gene set are coordinately up- or down-regulated within a sample.
  • single sample GSEA (ssGSEA) from http://baderlab.org/
  • How single sample GSEA works
  • Data formats gct (expression data format), gmt (gene set database format).
  • gsva() from the GSVA package has an option to compute ssGSEA. 单样本基因集富集分析 --- ssGSEA. The output of gsva() is a matrix of ES (# gene sets x # samples). It does not produce plots nor running the permutation tests. 
    library(GSVA)
library(heatmaply)
p <- 10 ## number of genes
n <- 30 ## number of samples
nGrp1 <- 15 ## number of samples in group 1
nGrp2 <- n - nGrp1 ## number of samples in group 2

## consider three disjoint gene sets
geneSets <- list(set1=paste("g", 1:3, sep=""),
                 set2=paste("g", 4:6, sep=""),
                 set3=paste("g", 7:10, sep=""))

## sample data from a normal distribution with mean 0 and st.dev. 1
y <- matrix(rnorm(n*p), nrow=p, ncol=n,
            dimnames=list(paste("g", 1:p, sep="") , paste("s", 1:n, sep="")))

## genes in set1 are expressed at higher levels in the last 'nGrp1+1' to 'n' samples
y[geneSets$set1, (nGrp1+1):n] <- y[geneSets$set1, (nGrp1+1):n] + 2

gsva_es <- gsva(y, geneSets, method="ssgsea")
dim(gsva_es) #  3 x 30
heatmaply(gsva_es)
heatmaply(y, Colv = F, Rowv= F, scale = "none")
  • A simple implementation of ssGSEA (single sample gene set enrichment analysis)
  • Use "ssgsea-gui.R". The first question is a folder containing input files GCT. The 2nd question is about gene set database in GMT format. This has to be very restrict. For example, "ptm.sig.db.all.uniprot.human.v1.9.0.gmt" and "ptm.sig.db.all.sitegrpid.human.v1.9.0.gmt" provided in github won't work with the example GCT file. 
    setwd("~/github/ssGSEA2.0/")
source("ssgsea-gui.R")
# select a folder containing gct files; e.g. PI3K_pert_logP_n2x23936.gct 
# select a gene set file; e.g. <ptm.sig.db.all.flanking.human.v1.8.1.gmt>

    A new folder (e.g. 2021-03-01) will be created under the same parent folder as the gct file folder. 

    tree -L 1 ~/github/ssGSEA2.0/example/gct/2021-03-20/                         

├── PI3K_pert_logP_n2x23936_ssGSEA-combined.gct
├── PI3K_pert_logP_n2x23936_ssGSEA-fdr-pvalues.gct
├── PI3K_pert_logP_n2x23936_ssGSEA-pvalues.gct
├── PI3K_pert_logP_n2x23936_ssGSEA-scores.gct
├── PI3K_pert_logP_n2x23936_ssGSEA.RData
├── parameters.txt
├── rank-plots
├── run.log
└── signature_gct

tree ~/github/ssGSEA2.0/example/gct/2021-03-20/rank-plots | head -3 
# 102 files. One file per matched gene set
├── DISEASE.PSP_Alzheime_2.pdf
├── DISEASE.PSP_breast_c_2.pdf

tree ~/github/ssGSEA2.0/example/gct/2021-03-20/signature_gct | head -3                    
# 102 files. One file per matched gene set
├── DISEASE.PSP_Alzheimer.s_disease_n2x23.gct
├── DISEASE.PSP_breast_cancer_n2x14.gct

    Since the GCT file contains 2 samples (the last 2 columns), ssGSEA produces one rank plot for each gene set (with adjusted p-value & the plot could contain multiple samples). The ES scores are saved in <PI3K_pert_logP_n2x23936_ssGSEA-scores> and adjust p-values are saved in <PI3K_pert_logP_n2x23936_ssGSEA-fdr-pvalues>.

    SsGSEA.png

  • Some discussions from biostars.org. Find -> "ssgsea"
  • Some papers.
  • 【生信分析 3】教你看懂GSEA和ssGSEA分析结果. No groups/classes in the data (6:33). Output is a heatmap. Each value is computed sample by sample. Rows = gene set. Columns = (sorted by the 1st gene set) samples.

Selected papers

Popularity and performance of bioinformatics software: the case of gene set analysis Xie 2021

Retrieved from "https://wiki.taichimd.us/index.php?title=GSEA&oldid=28868"