GSEA: Difference between revisions
Line 564: | Line 564: | ||
* [https://peerj.com/articles/4262/ Interpretation] | * [https://peerj.com/articles/4262/ Interpretation] | ||
** If the '''coefficient is positive''', miRNA inhibition on target mRNAs might be lifted/reduced, thereby leading to '''less suppression''' on the gene set of interest in the '''experimental group'''. | ** If the '''coefficient is positive''', miRNA inhibition on target mRNAs might be lifted/reduced, thereby leading to '''less suppression''' on the gene set of interest in the '''experimental group'''. | ||
*** When the statement mentions a ‘positive coefficient’, it’s referring to a situation where an increase in miRNA levels is associated with an increase in the levels of target mRNAs. This is unusual because '''we’d normally expect more miRNA to lead to less mRNA, due to the inhibition process''' I mentioned earlier. | *** When the statement mentions a ‘positive coefficient’, it’s referring to a situation where an increase in miRNA levels is associated with an increase in the levels of target mRNAs. ''This is unusual'' because '''we’d normally expect more miRNA to lead to less mRNA, due to the inhibition process''' I mentioned earlier. | ||
*** So, if the coefficient is positive, it might mean that the miRNA’s inhibitory effect on the mRNA is being lifted or reduced. As a result, there’s less suppression of the mRNA, which could lead to more of its associated protein being produced. | *** So, if the coefficient is positive, it might mean that the miRNA’s inhibitory effect on the mRNA is being lifted or reduced. As a result, there’s less suppression of the mRNA, which could lead to more of its associated protein being produced. | ||
*** If these genes are the ones producing the target mRNAs, then less suppression by the miRNAs would lead to '''increased activity of these genes in the experimental group.''' | *** If these genes are the ones producing the target mRNAs, then less suppression by the miRNAs would lead to '''increased activity of these genes in the experimental group.''' |
Revision as of 16:19, 20 October 2023
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
- https://www.bioconductor.org/help/course-materials/2015/SeattleApr2015/E_GeneSetEnrichment.html, https://www.bioconductor.org/help/course-materials/2009/SeattleApr09/gsea/HyperG_Lecture.pdf. For hypergeometric test, it is
- finding Biocarta or KEGG pathways significantly enriched in the user's gene list
- Are DE genes in the set more common than DE genes not in the set?.
- Are selected genes more often in the GO category than expected by chance? (one-tailed)
We can draw a 2x2 table or a venn diagram to see the diagram.
In gene set Yes No DE Yes No
- Gene List Enrichment Analysis dhyper(), binom.test(), fisher.test().
- Pathway Commons
- Gene set analysis methods: statistical models and methodological differences 2014
- http://software.broadinstitute.org/gsea/index.jsp, Subramanian, et al 2005 paper
- HOW TO PERFORM GSEA - A tutorial on gene set enrichment analysis for RNA-seq (video)
- Algorithm. GSEA walks down the ranked list of genes, increasing a running-sum statistic when a gene belongs to the set and decreasing it when the gene does not.
- Interpretation of 3 enrichment plots. The 1st plot (ES on y-axis) tells you how over or under expressed is your gene respect to the ranked list. The 2nd part of the graph (barcode-like) shows where the members of the gene set appear in the ranked list of genes. The 3rd graph (y=ranked list metric, x=rank) shows how your metric is distributed along the list.
- slides with formulas.
- What does a negative enrichment score mean? A negative NES will indicate that the genes in the set S will be mostly at the bottom of your list L.
- GSEA R Implementation from GSEA-MSigDB.
- READMe.
- Using GSEA.1.0.R
- Statistical power of gene-set enrichment analysis is a function of gene set correlation structure by SWANSON 2017
- clusterProfiler package and the online book.
- Gene-set Enrichment with Regularized Regression Fang 2019
- msigdbr package. MSigDB Gene Sets for Multiple Organisms in a Tidy Data Format.
- Best method/package for Gene Set Enrichment Analysis in R? and the gage package
- ES could be negative; see Genome 559: Introduction to Statistical and Computational Genomics
- multiGSEA: a GSEA-based pathway enrichment analysis for multi-omics data
- How to use DAVID for functional annotation of genes, Using DAVID for Functional Enrichment Analysis in a Set of Genes (Part 1), (Part 2) (video)
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.
Interpretation
- Gene Set Enrichment Analysis (GSEA) takes an alternative approach : it focuses on cumulative changes in expression of multiple genes as a group (belonging to a same gene-set/pahtway), which shifts the focus from individual genes to groups of genes. See this.
- XXX class/subtype is associated with the OOO gene set(s) (GSVA vignette)
- XXX subtype (of samples) is characterized by the expression of OOO markers, thus we expect it to correlate with the OOO2 gene set (GSVA vignette)
- The XXX subtype (of samples) shows high expression of OOO genes, thus the OOO gene set is highly enriched for this group (GSVA vignette). OR if we find a gene set is enriched in XXX subtype (e.g. sensitive models), it means genes in that gene set are more highly expressed in the XXX subtype samples compared to the YYY subtype samples.
- Negative ES Interpretation.
- A positive enrichment score should always reflect enrichment on the positive side of the zero cross (although not necessarily all genes on the positive side) and be enrichment in whichever phenotype was selected to be first in the comparison in the Phenotype selection dialogue. And vice versa, with a negative enrichment score reflecting enrichment in the genes on the negative side of the ranked list.
- GSEA is known to have some issues with highly skewed gene distributions, but that shouldn't affect the raw enrichment scores, just NES and significance calculations when GSEA runs in to a condition where its only sampled from one side of the distribution.
- Enrichment Score Interpretation
- If, for example, you provide a gene list ranked by a combination of fold change and p-value (e.g., sign(FC) * log10(pvalue)), then the positive scores are associated with upregulated genes and negative scores are associated with downregulated genes.
- fgsea paper
- The more positive is the value of sr(p) the more enriched the gene set is in the positively-regulated genes (with Si > 0). Accordingly, negative sr(p) corresponds to enrichment in the negatively regulated genes.
- The opposite of "enriched" is "depleted". See Wikipedia.
Calculation
- Subramanian paper
- https://youtu.be/bT00oJh2x_4
- *pathwaycommons.org
- Compute cumulative sum over ranked genes. MD Anderson lecture
- Increase sum when gene in set, decrease it otherwise. That is, +1/-1 weights in cumulative sum is used to represent whether genes are in the interested gene set.
- Magnitude of increment depends on correlation of gene with phenotype.
- Record the maximum deviation from zero as the enrichment score
FDR cutoff
Why does GSEA use a false discovery rate (FDR) of 0.25 rather than the more classic 0.05?
Given the lack of coherence in most expression datasets and the relatively small number of gene sets being analyzed, using a more stringent FDR cutoff may lead you to overlook potentially significant results.
piano
- piano package. Platform for integrative analysis of omics data.
- Available GSA methods.
- The main function is runGSA(). It contains a wide range of GSEA methods like maxmean, fisher, gsea, fgsea, etc.
- Conceptual overview of the runGSA workflow
- It was used in Xeva paper.
edgeR
Over-representation analysis. ?goana and ?kegga. See UserGuide 2.14 Gene ontology (GO) and pathway analysis.
fgsea
fgsea package and download stat
- ?exampleRanks gives a hint about how it was obtained. It is actually the t-statistic, not logFC.
... fit2 <- eBayes(fit2) de <- data.table(topTable(fit2, adjust.method="BH", number=12000, sort.by = "B"), keep.rownames = TRUE) colnames(de) plot(de$t, exampleRanks[match(de$rn, names(exampleRanks))]) # 45 degrees straight line range(de$t) # [1] -62.22877 52.59930 range(exampleRanks) # [1] -63.33703 53.28400 range(de$logFC) # [1] -11.56067 12.57562
- Paper is still on biorxiv only. Scholar.google.com. Computer tech lab in Russia.
- Documentation
- ?collapsePathways. Collapse list of enriched pathways to independent ones.
- DESeq results to pathways in 60 Seconds with the fgsea package
- Performing GSEA using MSigDB gene sets in R
- Using the fast preranked gene set enrichment analysis (fgsea) package
Are fgsea and Broad Institute GSEA equivalent
Are fgsea and Broad Institute GSEA equivalent?
Examples
- vignette,
- plotEnrichment() for a single pathway (including the output plot). No need to run fgsea(). Source code. Enrichment score (ES) on the plot is calculated by calcGseaStat()$res. The ES value is the same as the one shown in plotGseaTable() though plotEnrichment does not return it.
- plotGseaTable() for a bunch of selected pathways (including the output plot).
- Basic usage
library(fgsea) library(ggplot2) data(examplePathways) # a list of length 1457 (pathways) data(exampleRanks) # a vector of length 12000 (genes) set.seed(42) fgseaRes <- fgsea(pathways = examplePathways, stats = exampleRanks, minSize = 4, # default minSize=1, maxSize=Inf maxSize =10) # I used a very small maxSize in order to see details later # So the results here can't be compared with the default class(fgseaRes) # [1] "data.table" "data.frame" dim(fgseaRes) # [1] 386 8 fgseaRes[1:2, ] # pathway pval padj # 1: 1368092_Rora_activates_gene_expression 0.8630705 0.9386018 # 2: 1368110_Bmal1:Clock,Npas2_activates_circadian_gene_expression 0.4312268 0.8159487 # log2err ES NES size leadingEdge # 1: 0.05412006 -0.3087414 -0.6545067 5 11865,12753,328572,20787 # 2: 0.08312913 0.4209054 1.0360288 9 20893,59027,19883 sapply(examplePathways, length) |> range() # [1] 1 2366 # Question: minSize, maxSize mean 'matched' genes? Can't verify sum(sapply(examplePathways, length) < 3) # [1] 104 sum(sapply(examplePathways, length) > 10) # [1] 939 1457 - 104 - 939 # [1] 414 examplePathways2 <- lapply(examplePathways, function(x) x[x %in% names(exampleRanks)]) sapply(examplePathways2, length) |> range() # [1] 0 968 sum(sapply(examplePathways2, length) < 3) # [1] 232 sum(sapply(examplePathways2, length) > 10) # [1] 730 1457 - 232 - 730 # [1] 495 range(fgseaRes$ES) # [1] -0.8442416 0.9488163 range(fgseaRes$NES) # [1] -2.020695 2.075729 order(fgseaRes$ES)[1:5] # [1] 75 289 320 249 312 order(fgseaRes$NES)[1:5] # [1] 289 75 102 320 200 # choose the top gene set in order to zoom in & see the detail head(fgseaRes[order(pval), ], 1) # pathway pval padj log2err ES # 1: 5991601_Phosphorylation_of_Emi1 2.461461e-07 9.501241e-05 0.6749629 0.9472236 # NES size leadingEdge # 1: 2.082967 6 107995,12534,18817,67141,268697,56371 plotEnrichment(examplePathways[["5991601_Phosphorylation_of_Emi1"]], exampleRanks) # 12000 total genes though only 6 genes are matched in this top gene set. debug(plotEnrichment) # Browse[2]> toPlot # x y # 1 0 0.000000000 # 2 47 -0.003918626 # 3 48 0.308356695 # 4 322 0.285511939 # 5 323 0.452415897 # 6 407 0.445412396 # 7 408 0.593677191 # 8 447 0.590425565 # 9 448 0.730277162 # 10 617 0.716186784 # 11 618 0.833696389 # 12 638 0.832028888 # 13 639 0.947223612 # 14 12001 0.000000000 examplePathways[["5991601_Phosphorylation_of_Emi1"]] # [1] "12534" "18817" "56371" "67141" "107995" "268697" # [7] "434175" "102643276" match(examplePathways[["5991601_Phosphorylation_of_Emi1"]], names(exampleRanks)) # [1] 11678 11593 11362 11553 11953 11383 NA NA # Since there are 12000 total genes and exampleRanks have been sorted, # these matched genes have very high ranks. # See the plot below. length(exampleRanks) # [1] 12000 exampleRanks[1:10] # 170942 109711 18124 12775 72148 16010 11931 # -63.33703 -49.74779 -43.63878 -41.51889 -33.26039 -32.77626 -29.42328 # 13849 241230 665113 # -28.83475 -28.65111 -27.81583 sum(exampleRanks < 0) # [1] 6469 sum(exampleRanks > 0) # [1] 5531 plot(exampleRanks)
There are 6*2 points (excluding the 0 value at the starting and end) in the figure.
- The range of NES is not [-1, 1] but ES seems to be in [-1, 1]. The example sort pathways by NES but it can be adapted to sort by pval. Source code. Also ES and NES are not in the same order (see below) though we cannot tell it from the small plot. According to the manual,
- ES – enrichment score, same as in Broad GSEA implementation;
- NES – enrichment score normalized to mean enrichment of random samples of the same size (the number of genes in the gene set).
- Hπow to interpret NES(normalized enrichment score)?
- Negative normalized Enrichment Score (NES) in GSEA analysis
- question about GSEA about the sign of ES/NES. "Does this means in positive NES upregulated genes are enriched while in negative NES DN genes are enriched?"π Step 1. rank all genes in your list (L) based on "how well they divide the conditions". Step 2, you want to see whether the genes present in a gene set (S) are at the top or at the bottom of your list...or if they are just spread around randomly. A positive NES will indicate that genes in set S will be mostly represented at the top of your list L. a negative NES will indicate that the genes in the set S will be mostly at the bottom of your list L.
with(fgseaRes, plot(abs(ES), pval)) # relationship is not consistent for a lot of pathways with(fgseaRes, plot(abs(NES), pval)) # relationship is more consistent with(fgseaRes, plot(ES, NES)) head(order(fgseaRes$pval)) # [1] 121 112 289 68 187 188 head(rev(order(abs(fgseaRes$NES)))) # [1] 121 289 68 188 187 22
- The computation speed is FAST!
system.time(invisible(fgsea(examplePathways, exampleRanks))) # user system elapsed # 31.162 38.427 15.800
plotEnrichment()
- Understand plotEnrichment() function. Plot shape (concave or convex) requires the stats of genes from both in- and out- of the gene set. This affects the interpretation of the plot. The plot always starts with 0 and end with 0 in Y (enrichment score). plotEnrichment() does not return enrichment scores even it calculates them for the plot.
# Reduce the total number of genes i <- match(examplePathways[["5991601_Phosphorylation_of_Emi1"]], names(exampleRanks)) i <- na.omit(i) # Exclude 2 genes in the pathways but not in the data if (FALSE) { # Total genes = pathway genes will not work? plotEnrichment(examplePathways[["5991601_Phosphorylation_of_Emi1"]], exampleRanks[i]) # Error: GSEA statistic is not defined when all genes are selected } # PS values in 'exampleRanks' object are sorted from the smallest to the largest # Case 1: include 1 down-regulated gene for our total genes plotEnrichment(examplePathways[["5991601_Phosphorylation_of_Emi1"]], exampleRanks[c(1,i)]) # saved as fgseaPlotSmall.png exampleRanks[c(1,i)] # 170942 12534 18817 56371 67141 107995 268697 # -63.33703 15.16265 13.46935 10.46504 12.70504 28.36914 10.67534 # # Rank and add +1/-1, # 28.36 15.16 13.46 12.70 10.67 10.46 -63.33 # +1 +1 +1 +1 +1 +1 -1 debug(plotEnrichment) # Browse[2]> toPlot # x y # 1 0 0.0000000 # 2 0 0.0000000 # 3 1 0.3122753 # 4 1 0.3122753 # 5 2 0.4791793 # 6 2 0.4791793 # 7 3 0.6274441 # 8 3 0.6274441 # 9 4 0.7672957 # 10 4 0.7672957 # 11 5 0.8848053 # 12 5 0.8848053 # 13 6 1.0000000 # 14 8 0.0000000 x <-c(28.36, 15.16, 13.46, 12.70, 10.67, 10.46, -63.33) plot(0:7, c(0, cumsum(x)), type = "b") # Case 2: include another one up-regulated gene (so all are up-reg) for our total genes # However, the largest gene (added) is not in the gene set. # It will make the plot starting from a negative value (ignore 0 from the very left) plotEnrichment(examplePathways[["5991601_Phosphorylation_of_Emi1"]], exampleRanks[c(i, 12000)]) # saved as fgseaPlotSmall2.png exampleRanks[c(i, 12000)] # 12534 18817 56371 67141 107995 268697 80876 # 15.16265 13.46935 10.46504 12.70504 28.36914 10.67534 53.28400 # # Rank and add +1/-1 # 53.28 28.36 15.16 13.46 12.70 10.67 10.46 # -1 +1 +1 +1 +1 +1 +1
For figure 1 (case 1), it has 6 points. Figure 2 is the one I generate using plot() and exampleRanks. For figure 3 (case 2), it has 7 points (excluding the 0 value at the starting and end). The curve goes down first because the gene has the largest rank (statistic) and is not in the pathway. So the x-axis location is determined by the rank/statistics and going up or down on the y-axis is determined by whether the gene is in the pathway or not.
- Shifting exampleRanks values will change the enrichment scores because the directions of genes' statistics change?
plotEnrichment(examplePathways[[1]], exampleRanks) + labs(title="Programmed Cell Death") # like sin() with one period range(exampleRanks) # [1] -63.33703 53.28400 plotEnrichment(examplePathways[[1]], exampleRanks+64) + labs(title="Programmed Cell Death") # like a time series plot plotEnrichment(examplePathways[[1]], exampleRanks-54) + labs(title="Programmed Cell Death") # same problem if we shift exampleRanks to all negative
including NE, NES, pval, etc. So it is not purely rank-based.
set.seed(42) fgseaRes2 <- fgsea(pathways = examplePathways, stats = exampleRanks+2, # OR 2*(exampleRanks+2) minSize = 4, maxSize =10) identical(fgseaRes, fgseaRes2) # [1] FALSE R> fgseaRes[1:3, c("ES", "NES", "pval")] ES NES pval 1: -0.3087414 -0.6545067 0.8630705 2: 0.4209054 1.0360288 0.4312268 3: -0.4065043 -0.8617561 0.6431535 R> fgseaRes2[1:3, c("ES", "NES", "pval")] ES NES pval 1: 0.4554877 0.8284722 0.7201767 2: 0.5520914 1.1604366 0.2833333 3: -0.2908641 -0.5971951 0.9442724
However, scaling won't change anything (because scaling does not change the directions?).
set.seed(42) fgseaRes3 <- fgsea(pathways = examplePathways, stats = 2*exampleRanks, minSize = 4, maxSize =10) identical(fgseaRes, fgseaRes3) # [1] TRUE
plotGseaTable()
- Obtain enrichment scores by calling fgsea(); i.e. fgseaMultilevel().
examplePathways[["5991601_Phosphorylation_of_Emi1"]] |> length() # [1] 8 set.seed(1) fgsea2 <- fgsea(pathways = examplePathways["5991601_Phosphorylation_of_Emi1"], stats = exampleRanks) tibble::as_tibble(fgsea2) # A tibble: 1 × 8 # pathway pval padj log2err ES NES size leadi…¹ # <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <list> # 1 5991601_Phosphor… 9.39e-7 9.39e-7 0.659 0.947 2.05 6 <chr> # … with abbreviated variable name ¹leadingEdge plotGseaTable(examplePathways["5991601_Phosphorylation_of_Emi1"], exampleRanks, fgsea2)
The signs of NES/ES can be confirmed from the gene rank plot. Below is a plot from 10 enriched pathways
and a plot from 10 non-enriched pathways
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.
Kolmogorov Smirnov test statistic
- In Gene Set Enrichment Analysis (GSEA), the Kolmogorov-Smirnov test is used to calculate an enrichment score for each gene-set (pathway).
- This score reflects how often members (genes) included in that gene-set (pathway) occur at the top or bottom of the ranked data set.
- The ranking is based on a measure of differential gene expression.
- The enrichment score is calculated using the Kolmogorov-Smirnov statistic, which quantifies the distance between the empirical distribution function of a sample and the cumulative distribution function of a reference distribution, or between the empirical distribution functions of two samples.
- In GSEA, the empirical distribution function represents the ranked list of genes, and the reference distribution represents the gene-set (pathway) being tested for enrichment.
- The enrichment score is then used to determine whether the gene-set (pathway) is significantly enriched in the ranked list of genes.
- An example
- Let’s say we have a ranked list of 10 genes based on their differential expression between two conditions, with the most highly expressed genes at the top and the least expressed genes at the bottom. The list might look like this: GeneA, GeneB, GeneC, GeneD, GeneE, GeneF, GeneG, GeneH, GeneI, GeneJ.
- Now let’s say we have a gene-set (pathway) that contains three genes: GeneB, GeneD, and GeneH. We want to determine if this gene-set is significantly enriched in our ranked list of genes.
- To do this, we can use the Kolmogorov-Smirnov test to calculate an enrichment score for the gene-set. The test will calculate the maximum distance between the empirical distribution function of our ranked list of genes and the cumulative distribution function of our gene-set.
- The empirical distribution function of our ranked list of genes is a step function that increases by 1/10 at each gene in the list. The cumulative distribution function of our gene-set is also a step function that increases by 1/3 at each gene in the set.
- The Kolmogorov-Smirnov statistic will then calculate the maximum vertical distance between these two functions. If this distance is large, it suggests that the members of our gene-set (GeneB, GeneD, and GeneH) are occurring more often at the top or bottom of our ranked list of genes than would be expected by chance. This would indicate that our gene-set is significantly enriched in our ranked list of genes.
- To determine if the distance is statistically significant, we can compare the calculated KS statistic to a critical value obtained from a reference distribution. If the KS statistic is greater than the critical value, we can conclude that the gene-set is significantly enriched in the ranked list of genes.
- What is the reference distribution of KS statistic? In the context of Gene Set Enrichment Analysis (GSEA), the reference distribution for the Kolmogorov-Smirnov (KS) statistic is obtained by permuting the sample labels (resampling genes actually) and recalculating the enrichment score for the gene-set many times. This creates a null distribution of enrichment scores that represents what we would expect to see if the gene-set was not enriched in the ranked list of genes.
- The critical value for the KS statistic is then obtained from this null distribution, typically by taking a high percentile (e.g., the 95th percentile) of the distribution. If the observed KS statistic is greater than this critical value, we can conclude that the gene-set is significantly enriched in the ranked list of genes.
- Why do we need to rank genes first? The ranking of genes is also used in the calculation of the enrichment score for a gene-set. The enrichment score reflects how often members (genes) included in that gene-set (pathway) occur at the top or bottom of the ranked list of genes. If the members of a gene-set occur more often at the top or bottom of the ranked list of genes than would be expected by chance, it suggests that the gene-set is significantly enriched in the ranked list of genes.
- R
# Define a ranked list of genes based on their differential expression # between two conditions, with the most highly expressed genes at the top # and the least expressed genes at the bottom ranked_genes <- c("GeneA", "GeneB", "GeneC", "GeneD", "GeneE", "GeneF", "GeneG", "GeneH", "GeneI", "GeneJ") # Define a gene-set (pathway) that contains three genes gene_set <- c("GeneB", "GeneD", "GeneH") # Calculate the enrichment score for the gene-set using the KS statistic enrichment_score <- max(ks.test(ranked_genes, gene_set)$statistic) # Permute the sample labels and recalculate the enrichment score many times # to create a null distribution of enrichment scores set.seed(123) n_permutations <- 1000 null_distribution <- replicate(n_permutations, { permuted_genes <- sample(ranked_genes) ks.test(permuted_genes, gene_set)$statistic }) # Calculate the critical value for the KS statistic from the null distribution critical_value <- quantile(null_distribution, 0.95) # Determine if the gene-set is significantly enriched in the ranked list of genes if (enrichment_score > critical_value) { print("The gene-set is significantly enriched in the ranked list of genes.") } else { print("The gene-set is not significantly enriched in the ranked list of genes.") }
Rafael A Irizarry
Gene set enrichment analysis made simple
microRNA/miRNA
- https://en.wikipedia.org/wiki/MicroRNA 小分子核糖核酸
- MicroRNAs (miRNAs) are a class of small, non-coding RNA molecules (無法進一步轉譯成蛋白質的RNA) that play a role in the regulation of gene expression. They are typically 21-23 nucleotides in length and are found in plants, animals, and some viruses. miRNAs function by binding to complementary sequences in messenger RNA (mRNA) molecules, leading to the silencing of the target gene through the degradation of the mRNA or the inhibition of its translation into protein.
- A target gene is a gene whose expression is regulated by a specific regulatory molecule, such as a transcription factor or a microRNA. The regulatory molecule binds to specific sequences in the DNA or RNA of the target gene, leading to changes in the expression of that gene. For example, a transcription factor may bind to the promoter region of a target gene, leading to an increase or decrease in the transcription of that gene into messenger RNA (mRNA). Similarly, a microRNA may bind to the mRNA of a target gene, leading to the degradation of the mRNA or the inhibition of its translation into protein.
- MiRNAs are involved in a wide range of cellular processes, including cell growth, differentiation, development, and apoptosis, and have been implicated in many diseases
RBiomirGS
- RBiomirGS R package
require("RBiomirGS") # two input files # liver.csv # kegg.v5.2.entrez.gmt, 186 pathways raw <- read.csv(file = "liver.csv", header = TRUE, stringsAsFactors = FALSE) dim(raw) # [1] 85 3 raw[1:3,] # miRNA FC pvalue # 1 mmu-miR-let-7f-5p 0.58 0.034 # 2 mmu-miR-1a-5p 0.32 0.037 # 3 mmu-miR-1b-5p 0.50 0.001 # Target mRNA mapping # creating 95 csv files; e.g. "mmu-miR-1a-5p_mRNA.csv" rbiomirgs_mrnascan(objTitle = "mmu_liver_predicted", mir = raw$miRNA, sp = "mmu", addhsaEntrez = TRUE, queryType = "predicted", parallelComputing = TRUE, clusterType = "PSOCK") # > head -3 mmu-miR-10b-5p_mRNA.csv # "database","mature_mirna_acc","mature_mirna_id","target_symbol","target_entrez","target_ensembl","score" # "diana_microt","MIMAT0000208","mmu-miR-10b-5p","Bdnf","12064","ENSMUSG00000048482","1" # "diana_microt","MIMAT0000208","mmu-miR-10b-5p","Bcl2l11","12125","ENSMUSG00000027381","1" # GSEA # generating 3 files: # mirnascore.csv # mrnascore.csv # mirna_mrna_iwls_GS.csv rbiomirgs_logistic(objTitle = "mirna_mrna_iwls", mirna_DE = raw, var_mirnaName = "miRNA", var_mirnaFC = "FC", var_mirnaP = "pvalue", mrnalist = mmu_liver_predicted_mrna_hsa_entrez_list, mrna_Weight = NULL, gs_file = "kegg.v5.2.entrez.gmt", optim_method = "IWLS", p.adj ="fdr", parallelComputing = FALSE, clusterType = "PSOCK") > wc -l mirnascore.csv 86 mirnascore.csv > head -3 mirnascore.csv 13:46:47 "miRNA","S_mirna" "mmu-miR-let-7f-5p",1.46852108295774 "mmu-miR-1a-5p",1.43179827593301 > wc -l mrnascore.csv 15297 mrnascore.csv > head -3 mrnascore.csv 13:46:53 "EntrezID","S_mrna" "1",-3.20388726317687 "10",-11.7197772358718 > wc -l mirna_mrna_iwls_GS.csv 187 mirna_mrna_iwls_GS.csv > head -3 mirna_mrna_iwls_GS.csv # "GS","converged","loss","gene.tested","coef","std.err","t.value","p.value","adj.p.val" # "KEGG_GLYCOLYSIS_GLUCONEOGENESIS","Y",0.0231448486023595,54,0.103936359331543,...
# Histogram of all gene sets # y-axis = model coefficients # x-axis = gene set rbiomirgs_bar(gsadfm = mirna_mrna_iwls_GS, n = "all", y.rightside = TRUE, yTxtSize = 8, plotWidth = 300, plotHeight = 200, xLabel = "Gene set", yLabel = "model coefficient") # Histogram of top 50 ranked KEGG pathways based on absolute value of the coefficient # (with significant adjusted p values < 0.05) # y-axis = gene set # x-axis = log odds ratio rbiomirgs_bar(gsadfm = mirna_mrna_iwls_GS, gs.name = TRUE, n = 50, y.rightside = FALSE, yTxtSize = 8, plotWidth = 200, plotHeight = 300, xLabel = "Log odss ratio", signif_only = TRUE) # Volcano plot of all gene sets # y-axis = -log10(p-value) # x-axis = model coefficient rbiomirgs_volcano(gsadfm = mirna_mrna_iwls_GS, topgsLabel = TRUE, n =15, gsLabelSize = 3, sigColour = "blue", plotWidth = 250, plotHeight = 220, xLabel = "model coefficient")
, ,
- miRNAs are small molecules that can bind to mRNAs and prevent them from being translated into proteins. This process is known as ‘inhibition’.
- Interpretation
- If the coefficient is positive, miRNA inhibition on target mRNAs might be lifted/reduced, thereby leading to less suppression on the gene set of interest in the experimental group.
- When the statement mentions a ‘positive coefficient’, it’s referring to a situation where an increase in miRNA levels is associated with an increase in the levels of target mRNAs. This is unusual because we’d normally expect more miRNA to lead to less mRNA, due to the inhibition process I mentioned earlier.
- So, if the coefficient is positive, it might mean that the miRNA’s inhibitory effect on the mRNA is being lifted or reduced. As a result, there’s less suppression of the mRNA, which could lead to more of its associated protein being produced.
- If these genes are the ones producing the target mRNAs, then less suppression by the miRNAs would lead to increased activity of these genes in the experimental group.
- Furthermore, with a positive coefficient, a unit increase in Smrna results in an increased odds ratio of a gene belonging to the gene set of interest.
- It needs to be clarified that a positive model coefficient for a gene set means that the gene set of interest might be under more miRNA-dependent inhibition in the control group, as opposed to being activated under the experimental condition.
- If the coefficient is positive, miRNA inhibition on target mRNAs might be lifted/reduced, thereby leading to less suppression on the gene set of interest in the experimental group.
- miRBase: the microRNA database
- One microRNA/miR is one set. For human, its name (miRBase ID) looks like hsa-miR-181a, miRBase Accession Number is MIMAT0000256.
- MicroRNA and microRNA target database
- Correlation of expression profiles between microRNAs and mRNA targets using NCI-60 data
MIENTURNET
- http://userver.bio.uniroma1.it/apps/mienturnet/
- MIENTURNET: an interactive web tool for microRNA-target enrichment and network-based analysis 2019
multiGSEA
multiGSEA - Combining GSEA-based pathway enrichment with multi omics data integration.
Plot
- Chapter 12 Visualization of Functional Enrichment Result, enrichplot::gseaplot() or enrichplot::gseaplot2() from Bioconductor
- GSEA 从数据到美图 where step4output.Rdata is available on github. Some quirks:
# geneset$term = str_remove(geneset$term,"HALLMARK_")
- fgsea::plotGseaTable() from Bioc
- limma:barcodeplot(), It is of course not the same as the Broad Institute's GSEA plot.
- 15 Visualization of functional enrichment result. Especially bar plots.
- Plot rank densities from singscore::plotRankDensity() from Bioc
Single sample
singscore
- singscore
- Paper Single sample scoring of molecular phenotypes. Our approach does not depend upon background samples and scores are thus stable regardless of the composition and number of samples being scored. In contrast, scores obtained by GSVA, z-score, PLAGE and ssGSEA can be unstable when less data are available (NS < 25). Simulation was conducted.
- SingscoreAMLMutations
- TotalScore = UpScore + DownScore
- centerScore and knownDirection(=TRUE by default) parameters used in generateNull() and simpleScore() functions.
- On the paper, epithelial and mesenchymal gene sets are up-regulated and TGFβ-EMT signature is bidirectional.
- The pvalue calculation seems wrong. For example, the first sample D_Ctrl_R1 returns p=0.996 but it should be 2*(1-.996) according to the null distribution plot. However, the 1st sample is a control sample so we don't expect it has a large score; its p-value should be large?
- Gene-set enrichment analysis workshop. Also use ExperimentHub, SummarizedExperiment, emtdata, GSEABase, msigdb, edgeR, limma, vissE, igraph & patchwork packages.
- Functions in the GSEABase package help with reading, parsing and processing the signatures.
ssGSEA & GSVA package
- 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).
- gmt file can be read by GSA.read.gmt()
- gmt file can be read by clusterProfiler::read.gmt(). See R做GSEA富集分析
- ssGSEA produced by Broad Gene Pattern is different from one produced by gsva package.
- Actually they are not compatible. See a plot in singscore paper.
- Gene Pattern also compute p-values testing the Spearman correlation
- GSVA and SSGSEA for RNA-Seq TPM data.
- Generally, negative enrichment values imply down-regulation of a signature / pathway; whereas, positive values imply up-regulation.
- The idea is to then conduct a differential signature / pathway analysis (using, for example, limma) so that you can have, in addition to differentially expressed genes, differentially expressed signatures / pathways.
- GSVA vignette. In GSEA Subramanian et al. (2005) it is also observed that the empirical null distribution obtained by permuting phenotypes is bimodal and, for this reason, significance is determined independently using the positive and negative sides of the null distribution.
- Tips:
- Pre-processing RNA-seq data (normalization and transformation) for GSVA
- pre-filtering gene sets for GSVA/ssGSEA
- kcdf parameter it is only relevant when method="gsva". See Method and kcdf arguments in gsva package and debugging source code.
- mx.diff: Offers two approaches to calculate the enrichment statistic (ES) from the KS random walk statistic. mx.diff=FALSE: ES is calculated as the maximum distance of the random walk from 0. mx.diff=TRUE (default): ES is calculated as the magnitude difference between the largest positive and negative random walk deviations.
- min.sz=1, max.sz=Inf
- 纯R代码实现ssGSEA算法评估肿瘤免疫浸润程度. The pheatmap package was used to draw the heatmap. Original paper.
- fpkm expression data was sorted by gene name & median values. Duplicated genes are removed. The data is transformed by log2(x+1) before sent to gsva()
- NES scores are scaled and truncated to [-2, 2]. The scores are further scaled to have the range [0,1] before sending it to the heatmap function
- gsva() from the GSVA package has an option to compute ssGSEA. 单样本基因集富集分析 --- ssGSEA (it includes the formula from Barbie's paper). The output of gsva() is a matrix of ES (# gene sets x # samples). It does not produce plots nor running the permutation tests. Notice the option ssgsea.norm. Note ssgsea.norm = TRUE (default) option will scale ES by the absolute difference of the max and min ES; Unscaled ES / (max(unscaled ES) - min(unscaled ES)). So the scaled ES values depends on the included samples. But it seems the impact of the included samples is small from some real data. See the source code on github and on how to debug an S4 function.
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 set.seed(1234) 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 hist(gsva_es) # bi-modal range(gsva_es) # [1] -0.4390651 0.5609349 ## build design matrix design <- cbind(sampleGroup1=1, sampleGroup2vs1=c(rep(0, nGrp1), rep(1, nGrp2))) ## fit the same linear model now to the GSVA enrichment scores fit <- lmFit(gsva_es, design) ## estimate moderated t-statistics fit <- eBayes(fit) ## set1 is differentially expressed topTable(fit, coef="sampleGroup2vs1") # logFC AveExpr t P.Value adj.P.Val B # set1 0.5045008 0.272674410 8.065096 4.493289e-12 1.347987e-11 17.067380 # set2 -0.1474396 0.029578749 -2.315576 2.301957e-02 3.452935e-02 -4.461466 # set3 -0.1266808 0.001380826 -2.060323 4.246231e-02 4.246231e-02 -4.992049 heatmaply(gsva_es) # easy to see a pattern # samples' clusters are not perfect heatmaply(gsva_es, scale = "none") # 'scale' is not working? heatmaply(y, Colv = F, Rowv= F, scale = "none") # not easy to see a pattern
- If we set ssgsea.norm=FALSE, do we get the same results when we compute ssgsea using a subset of samples?
gsva1 <- gsva(y, geneSets, method="ssgsea", ssgsea.norm = FALSE) gsva2 <- gsva(y[, 1:2], geneSets, method="ssgsea", ssgsea.norm = FALSE) range(abs(gsva1[, 1:2] - gsva2)) # [1] 0 0
- Verify ssgsea.norm = TRUE option.
gsva_es <- gsva(y, geneSets, method="ssgsea") gsva1 <- gsva(y, geneSets, method="ssgsea", ssgsea.norm = FALSE) gsva2 <- gsva1 / abs(max(gsva1) - min(gsva1)) # abs(max(gsva1) - min(gsva1)) 8.9 range(gsva_es) # [1] -0.4390651 0.5609349 range(abs(gsva_es - gsva2)) # [1] 0 0
- Only ranks matters! If I replace the sample 1 gene expression values with the ranks, the ssgsea scores are not changed at all.
y2 <- y y2[, 1] <- rank(y2[, 1]) gsva2 <- gsva(y2, geneSets, method="ssgsea") gsva_es[, 1] # set1 set2 set3 # 0.1927056 0.1699799 -0.1782934 gsva2[, 1] # set1 set2 set3 # 0.1927056 0.1699799 -0.1782934 all.equal(gsva_es, gsva2) # [1] TRUE # How about the reverse ranking? No that will change everything. # The one with the smallest value was assigned one according to # the definition of 'rank' y3[, 1] <- 11 - y2[, 1] gsva3 <- gsva(y3, geneSets, method="ssgsea") gsva3[, 1] # set1 set2 set3 #-0.054991319 -0.009603802 0.191733634 cbind(y[, 1], y2[, 1], y3[, 1]) # [,1] [,2] [,3] # g1 -1.2070657 2 9 # g2 0.2774292 7 4 # g3 1.0844412 10 1 # g4 -2.3456977 1 10 # g5 0.4291247 8 3 # g6 0.5060559 9 2 # g7 -0.5747400 4 7 # g8 -0.5466319 6 5 # g9 -0.5644520 5 6 # g10 -0.8900378 3 8
- Is mx.diff useful? No. GSVA with and without absolute scaling method. The mx.diff and abs.ranking only apply when method="gsva". The only parameters that apply when method="ssgsea" are tau and ssgsea.norm.
gsva4 <- gsva(y, geneSets, method="ssgsea", mx.diff = FALSE) all.equal(gsva_es, gsva4) # [1] TRUE
- 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>.
- Some possible problems using ssGSEA to replace GSEA for DE analysis? See a toy example on Single Sample GSEA (ssGSEA) and dynamic range of expression. PS: the rank values table is wrong; they should be reversed.
- In general terms, PLAGE and z-score are parametric and should perform well with close-to-Gaussian expression profiles, and ssGSEA and GSVA are non-parametric and more robust to departures of Gaussianity in gene expression data. See Method and kcdf arguments in gsva package
- 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.
- Gene set analysis: GSVA, Z-Score, and fold-changes
escape
escape - Easy single cell analysis platform for enrichment. Github.
sTAM
- Paper sTAM: An Online Tool for the Discovery of miRNA-Set Level Disease Biomarkers 2020
- http://mir.rnanut.net/stam
Benchmark
Toward a gold standard for benchmarking gene set enrichment analysis Geistlinger, 2021. GSEABenchmarkeR package in Bioconductor.
phantasus
- phantasus Visual and interactive gene expression analysis. read.gct(), write.gct(), getGSE(), loadGEO(), gseaPlot().
- Using Phantasus application
PhenoExam
PhenoExam: gene set analyses through integration of different phenotype databa
Selected papers
- On the influence of several factors on pathway enrichment analysis Mubeen 2022
- Popularity and performance of bioinformatics software: the case of gene set analysis Xie 2021
- Gene set analysis approaches for RNA-seq data: performance evaluation and application guideline Rahmatallah 2015. Simulated RNA-Seq data was considered too.
- Gene set analysis methods: a systematic comparison Mathur, 2018
- A Comparison of Gene Set Analysis Methods in Terms of Sensitivity, Prioritization and Specificity Tarca 2013
Malacards
- MalaCards - human disease database. Paper 2018
- GSEABenchmarkeR package used it.
- (Not Malacards) KEGGandMetacoreDzPathwaysGEO Disease Datasets from GEO. This is a collection of 18 data sets for which the phenotype is a disease with a corresponding pathway in either KEGG or metacore database.This collection of datasets were used as gold standard in comparing gene set analysis methods.
- 生信数据库之疾病数据库
GSDA
Gene-set distance analysis (GSDA): a powerful tool for gene-set association analysis
Genekitr
https://github.com/GangLiLab/genekitr. twitter
MSigDB
https://www.gsea-msigdb.org/gsea/msigdb/
- C1: Hallmark
- C2: Curated gene sets including BioCarta, KEGG, Reactome
- C5: Oncology
Hallmark
- HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION/EMT from MsigDB
- Epithelial and mesenchymal gene signatures. See here. I save the gene lists in Github.
- 'Table S1A. Generic EMT signature for tumour ' is really more interesting since Epi is more tumor related. The first 145 are epi, the rest (315-145=170) are mes.
- Epithelial (170 genes), mesenchymal gene signatures 218-170=48 genes). Tan et al. 2014 and can be found in the 'Table S1B. Generic EMT signature for cell line ’.