Gene Set Analysis
BIBS Workshop
Training materials were derived or adapted from excellent online tutorials. Don’t hesitate to explore these resources to complete this training:
1 Preamble : Rolling your own statistics
1.1 Statistics by hand : “Is the red color enriched?”
Hint
Do you keep the same answer if I add the bag where the marbles have been drawn ?
So, is the red color enriched?
It depends on the background!
Without a background set, you cannot determine if a category is over-represented in your gene list. The same number of red balls could be:
- Enriched — if the background has only a few red balls
- Not enriched — if the background has many red balls
- Depleted — if the background has even more red balls than expected
The statistical test (hypergeometric test) compares the proportion of genes in your category against the proportion of genes in that category in the full genome.
For example: - Your gene list has 10 genes, 4 of which are in the “red category” - The background (all genes) has 13,000 genes, 35 of which are in the “red category” - The question: is 4/10 = 40% significantly more than 35/13,000 ≈ 0.27%?
That is what Over-Representation Analysis (ORA) answers — it compares the proportion of a category in your list to the proportion of that same category in the background (universe).
1.2 Start again: “Is the red color enriched?”
Hint
What do you intuitively use as information to respond ? What is the difference with the previous cases ?
So, is the red color enriched?
In this case, you use the size (a ranking metric) to order the marbles from biggest to smallest. By sorting all marbles by size, you can see whether the red marbles (the gene set of interest) are clustered at one end (the big end) — that is enrichment — or spread evenly across all sizes — not enriched.
- A2 (Enriched) : The red marbles are mostly on the left (biggest side), meaning the gene set members are concentrated at the top of the ranked list.
- B2 (Not enriched) : The red marbles are randomly scattered across all sizes, meaning the gene set members are not coordinated with the ranking.
This is the Gene Set Enrichment Analysis (GSEA) principle : you do not need a threshold (p-value / significant vs non-significant). Instead, you use all genes ranked by a metric (e.g., log fold change), and test whether a gene set’s members are preferentially at the top or bottom of that ranked list. This approach is more powerful when you have only a small number of significantly DE genes, because it can capture weak but coordinated signals across a whole pathway — while ORA would miss them entirely since the individual genes never reached the significance cutoff.
1.3 Summary
First take-home message
- ORA (Over-Representation Analysis): You need a gene list (e.g., significant DE genes) + a background (universe). It tests if a category is over-represented — i.e., whether the proportion in your list is higher than expected based on the background.
- FCS/GSEA (Ranking-based): You need all genes ranked by a metric (e.g., fold change). It tests if a gene set’s members are clustered at one end of the ranked list — no arbitrary threshold needed.
For a robust and reproducible analysis, you always need to know: - What is your background? (What set of marbles/genes could you possibly draw?) - What database was used? (Which categories were tested?) - Without this information, enrichment results cannot be interpreted — and therefore cannot be reproduced.
2 Overview of Gene set analysis methods
Generally for any differential expression analysis, it is useful to interpret the resulting genes lists using freely available web- and R-based tools. To gain greater biological insight on the differentially expressed genes (DEGs) you may want to determine whether there is enrichment of known biological functions, interactions, or pathways. While tools for functional analysis span a wide variety of techniques, they can loosely be categorized into three main types: over-representation analysis, functional class scoring, and pathway topology [1].
In this training, we will only focus on the first two approaches. (If you are interested in pathway topology analysis, you can find a review here and a review of topology-based methods here.)
Depending on the method you wish to use, you will need the corresponding input data.
- For over-representation based methods, you need a gene list (DE, up-regulated genes, down-regulated genes, …) of interest and a background list (the list of genes used as control - “the universe”).
- For ranking-based methods, you need all studied genes ranked according to a score (for example, the log fold change).
Note that all tools described below are great tools to validate experimental results and to make hypotheses. These tools suggest genes/pathways that may be involved with your condition of interest; however, you should NOT use these tools to make conclusions about the pathways involved in your experimental process. You will need to perform experimental validation of any suggested pathways.
2.1 Over-Representation Analysis (ORA)
There are a plethora of functional enrichment tools that perform some type of “over-representation” analysis by querying databases containing information about gene function and interactions.
These databases typically categorize genes into groups (gene sets) based on shared function, or involvement in a pathway, or presence in a specific cellular location, or other categorizations, e.g. functional pathways, etc. Essentially, known genes are binned into categories that have been consistently named (controlled vocabulary) based on how the gene has been annotated functionally. These categories are independent of any organism, however each organism has distinct categorizations available.
To determine whether any categories are over-represented, you can determine the probability of having the observed proportion of genes associated with a specific category in your gene list based on the proportion of genes associated with the same category in the background set (gene categorizations for the appropriate organism).
The statistical test that will determine whether something is actually over-represented is the Hypergeometric test. Using the example of the first functional category above, hypergeometric distribution is a probability distribution that describes the probability of 25 genes (k) being associated with “Functional category 1”, for all genes in our gene list (n=1000), from a population of all of the genes in entire genome (N=13,000) which contains 35 genes (K) associated with “Functional category 1”.
The calculation of probability of k successes follows the formula:
This test will result in an adjusted p-value (after multiple test correction) for each category tested.
2.2 Ranking-based method/Functional Class Scoring (FCS) (ex: GSEA)
In addition to over-representation analysis, other types of analyses can be equally important or informative for obtaining biological insight from your results. The hypothesis behind functional class scoring (FCS) methods is that although large changes in individual genes can have significant effects on pathways (and will be detected via ORA methods), weaker but coordinated changes in sets of functionally related genes (i.e., pathways) can also have significant effects.
Thus, rather than setting an arbitrary threshold to identify ‘significant genes’, all genes are considered in the analysis. The gene-level statistics from the dataset are aggregated to generate a single pathway-level statistic, and the statistical significance of each pathway is reported. This type of analysis can be particularly helpful if the differential expression analysis only outputs a small list of significant DE genes.
2.2.1 Gene set enrichment analysis using clusterProfiler and Pathview
One commonly used tool which is classified under Functional class scoring (FCS), is GSEA. Gene set enrichment analysis utilizes gene-level statistics or log2 fold changes for all genes to determine whether gene sets for particular biological pathways are enriched among the large positive or negative fold changes.
Gene sets are pre-defined groups of genes, which are functionally related. Commonly used gene sets include those derived from KEGG pathways, Gene Ontology terms, MSigDB, Reactome, or gene groups that share some other functional annotations, etc.
This image describes the theory of GSEA, with the ‘gene set S’ showing the metric used (in our case, ranked log2 fold changes) to determine enrichment of genes in the gene set. The left-most image is representing this metric used for the GSEA analysis. The log2 fold changes for each gene in the ‘gene set S’ is shown as a line in the middle image.
The large positive log2 fold changes are at the top of the gene set image, while the largest negative log2 fold changes are at the bottom of the gene set image. In the right-most image, the gene set is turned horizontally, underneath which is an image depicting the calculations involved in determining enrichment, as described below.
Step 1: Calculation of enrichment score:
An enrichment score for a particular gene set is calculated by walking down the list of log2 fold changes and increasing the running-sum statistic every time a gene in the gene set is encountered and decreasing it when genes are not part of the gene set. The size of the increase/decrease is determined by magnitude of the log2 fold change.
Larger (positive or negative) log2 fold changes will result in larger increases or decreases. The final enrichment score is where the running-sum statistic is the largest deviation from zero.
Step 2: Estimation of significance:
The significance of the enrichment score is determined using permutation testing, which performs rearrangements of the data points to determine the likelihood of generating an enrichment score as large as the enrichment score calculated from the observed data. Essentially, for this step, the first permutation would reorder the log2 fold changes and randomly assign them to different genes, reorder the gene ranks based on these new log2 fold changes, and recalculate the enrichment score. The second permutation would reorder the log2 fold changes again and recalculate the enrichment score, and this would continue for the total number of permutations run. Therefore, the number of permutations run will increase the confidence in the significance estimates.
Step 3: Adjust for multiple test correction
After all gene sets are tested, the enrichment scores are normalized for the size of the gene set, then the p-values are corrected for multiple testing.
The Normalized Enrichment Score (NES) is the enrichment score after normalization for the size of each gene set: larger gene sets tend to have higher raw enrichment scores simply because they contain more genes, so NES corrects for this bias and makes scores comparable across different gene set sizes. NES is the primary metric shown in GSEA dotplots (x-axis).
The GSEA output will yield the core genes in the gene sets that most highly contribute to the enrichment score. The genes output are generally the genes at or before the running sum reaches its maximum value (e.g., the most influential genes driving the differences between conditions for that gene set).
2.3 Gene sets databases
A gene set is an unordered collection of genes that are functionally related. A pathway can be interpreted as a gene set by ignoring functional relationships among genes.
2.3.1 Gene Ontology (GO)
Gene Ontology defines concepts/classes used to describe gene function, and relationships between these concepts. It classifies functions along three aspects:
- MF: Molecular Function molecular activities of gene products
- CC: Cellular Component where gene products are active
- BP: Biological Process pathways and larger processes made up of the activities of multiple gene products
GO terms are organized in a directed acyclic graph, where edges between terms represent parent-child relationship.
GO term hierarchy
Some gene products are well-researched, with vast quantities of data available regarding their biological processes and functions. However, other gene products have very little data available about their roles in the cell.
For example, the protein “p53” would contain a wealth of information about its roles in the cell, whereas another protein might only be known as a “membrane-bound protein” with no other information available.
The GO ontologies were developed to describe and query biological knowledge with differing levels of information available. To do this, GO ontologies are loosely hierarchical, ranging from general, ‘parent’, terms to more specific, ‘child’ terms. The GO ontologies are “loosely” hierarchical since ‘child’ terms can have multiple ‘parent’ terms.
Some genes with less information may only be associated with general ‘parent’ terms or no terms at all, while other genes with a lot of information be associated with many terms.
2.3.2 Kyoto Encyclopedia of Genes and Genomes (KEGG)
KEGG is a collection of manually drawn pathway maps representing molecular interaction and reaction networks. These pathways cover a wide range of biochemical processes that can be divided into 7 broad categories:
- Metabolism
- Genetic information processing
- Environmental information processing
- Cellular processes
- Organismal systems
- Human diseases
- Drug development.
Available pathways are listed: here.
2.3.3 Other gene sets
GO and KEGG are the most frequently used for functional analysis. They are typically the first choice because of their long-standing curation and availability for a wide range of species.
Other gene sets include but are not limited to Disease Ontology (DO), Disease Gene Network (DisGeNET), wikiPathways, Molecular Signatures Database (MSigDb).
2.3.4 User defined gene sets (GMT files)
The GMT file format is a tab delimited file format that describes gene sets. In the GMT format, each row represents a gene set; in the GMX format, each column represents a gene set.
3 Practicals
3.1 Get “toy” data
In this training, we will use the results of a differential expression analysis described in Kenny PJ et al, Cell Rep 2014.
The RNA-Seq was performed on HEK293F cells that were either transfected with a MOV10 transgene, or treated with siRNA to knock down Mov10 expression, or with non-specific (irrelevant) siRNA. This resulted in 3 conditions: Mov10 oe (overexpression), Mov10 kd (knockdown) and Irrelevant kd, respectively. The number of replicates is as shown below.
Using these data, the authors evaluate transcriptional patterns associated with perturbation of MOV10 expression. Please note that the irrelevant siRNA will be treated as our control condition.
What is the purpose of these datasets? What does Mov10 do?
The authors are investigating interactions between various genes involved in Fragile X syndrome, a disease in which there is aberrant production of the FMRP protein.
FMRP is “most commonly found in the brain, is essential for normal cognitive development and female reproductive function. Mutations of this gene can lead to fragile X syndrome, mental retardation, premature ovarian failure, autism, Parkinson’s disease, developmental delays and other cognitive deficits.” - from wikipedia
MOV10, is a putative RNA helicase that is also associated with FMRP in the context of the microRNA pathway.
The hypothesis tested in the paper is that FMRP and MOV10 associate and regulate the translation of a subset of RNAs.
3.2 Result files
Up regulated genes between MOV10 over expressed samples and control samples: here
Down regulated genes between MOV10 over expressed samples and control samples: here
Differentially (up + down) expressed genes between MOV10 over expressed samples and control samples: here
All the studied genes = background genes list here
The complete result of the differential expression analysis, with gene identifiers, logfoldchange values and pvalue adjusted: here
Download all these files to your computer.
3.3 Web tools
In the first part, we will conduct an over-representation analysis (ORA) using a web app.
There are several solutions, including gprofiler, DAVID, …
Here, we will use gprofiler: https://biit.cs.ut.ee/gprofiler/gost
3.3.1 gprofiler
- Go to the gprofiler website and make an analysis on the DEG list
- Upload the “DEG_list.txt”
- On the right, in Data sources select only GO databases and KEGG
- Click on Run query
- In a new page, do the same thing but add a background list
- Upload the “DEG_list.txt”
- On the right, in Data sources select only GO databases and KEGG
- Click on Advanced options
- Choose Custom in Statistical domain scope and upload the background_list.txt
- Click on Run query
- Compare the results
Depending on the background list, ORA results can vary considerably!
So, for a reproducible science, it is important to provide:
- the used tools/app and its version
- the code or the instructions
- the used databases and its versions
- the genes list of interest and the background
3.4 Using a more reproducible way (command lines)
Now we will perform the same analysis in a more reproducible way. For this, we will use the clusterProfiler R package, which combines:
- Over-Representation Analysis
- Ranking-based method
- Most used databases
- and visualization functions!
- Create a working directory (GSA_BIBS_training)
- Download the code template file available here and put it in “GSA_BIBS_training”
- Create a “data” directory and a “results” directory
- In the data directory, place the complete result “DE_result_table.tsv” from the differential expression analysis available here
Your working directory should look like this:
GSA_BIBS_training/
├── data/
│ └──DE_complete_results.tsv
├── GSA_template.Rmd
└── results/Depending on the parameters used and/or the methods used, Gene Set Analysis results can vary considerably!
-> keeping track of all the commands/instructions that generate your results is very important (the use of a script is strongly recommended)
All tools described today are great tools to validate experimental results and to make hypotheses. These tools suggest genes/pathways that may be involved with your condition of interest; however, you should NOT use these tools to make conclusions about the pathways involved in your experimental process. You will need to perform experimental validation of any suggested pathways.
3.5 Time to work on your data
4 Ressources
- Tips for working with GO terms: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003343
- Urgent need for consistent standards in functional enrichment analysis: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1009935
- Gene Set Analysis: Challenges, Opportunities, and Future Research: https://www.frontiersin.org/articles/10.3389/fgene.2020.00654/full
- Interpreting omics data with pathway enrichment analysis: https://doi-org.insb.bib.cnrs.fr/10.1016/j.tig.2023.01.003
These materials have been developed by members of the BIBS team at CIRI (https://ciri.ens-lyon.fr/). These are open access materials distributed under the terms of the Creative Commons Attribution license (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
- Some materials used in these lessons were derived or adapted from work made available by the Harvard Chan Bioinformatics Core (HBC) (https://github.com/hbctraining) under the Creative Commons Attribution license (CC BY 4.0).