Functional Enrichment Analysis Study Guide

RNA-Seq Data Analysis Framework

Pathway enrichment analysis is the final stage of a multi-step RNA-seq data analysis pipeline. This pipeline involves the following sequence:

Raw FASTQ files: Initial sequencing data undergo quality control measures.
Trimmed and Filtered FASTQ files: Processed sequences are aligned to the reference genome.
Genome Alignment: Aligned reads are used for gene expression quantification.
Count Matrix: A matrix comparing samples (e.g., Sample A, B, C, D) against genes (e.g., Gene 1, Gene 2, Gene 3, Gene 4).
Normalized Counts: Adjustment of counts to allow for meaningful gene expression comparisons.
Differentially Expressed Genes (DEGs): Identification of genes with significant changes in expression between conditions, often visualized using volcano plots showing the relationship between $\log_{2}(\text{Fold Change})$ and $-log_{10}(p\text{-value})$ .
Pathway Enrichment Analysis: Interpreting the biological meaning behind the list of differentially expressed genes.

Understanding Differentially Expressed Genes (DEGs)

A standard output from a differential expression study includes a table of genes characterized by various statistical metrics. For example, a result table might include:

ENTREZID & SYMBOL: Unique identifiers for genes (e.g., 55753, OGDHL).
GENENAME: Full descriptive name (e.g., oxoglutarate dehydrogenase L).
logFC (Log Fold Change): Magnitude of change (e.g., OGDHL: $3.294304$ , NES: $-4.649815$ ).
P.Value: Statistical significance (e.g., OGDHL: $1.361576 \times 10^{-5}$ , NES: $1.419162 \times 10^{-5}$ , LHX6: $1.661077 \times 10^{-5}$ , SH2D5: $3.324459 \times 10^{-5}$ , RPL9: $-4.186648$ ).
adj.P.Val: The p-value adjusted for multiple testing (e.g., $0.08325319$ for the top genes listed).
t-statistic and B-statistic: Indicators of the strength and probability of differential expression.

Specific examples of DEGs identified in RNA-seq include:

OGDHL (oxoglutarate dehydrogenase L): logFC $3.294304$ .
NES (nestin): logFC $-4.649815$ .
LHX6 (LIM homeobox 6): logFC $4.756681$ .
SH2D5 (SH2 domain containing 5): logFC $3.815348$ .
RPL9 (ribosomal protein L9): logFC $-4.186648$ .
RAB11FIP5 (RAB11 family interacting protein 5): logFC $2.757163$ .
NPM2 (nucleophosmin/nucleoplasmin 2): logFC $3.095437$ .
GAD2 (glutamate decarboxylase 2): logFC $4.176362$ .
PPP1R1B (protein phosphatase 1 regulatory inhibitor subunit 1B): logFC $4.176593$ .
SLC25A22 (solute carrier family 25 member 22): logFC $2.784811$ .
F2R (coagulation factor II thrombin receptor): logFC $-3.427165$ .

Pathway Enrichment Analysis (PEA) Overview

Pathway enrichment analysis summarizes long gene lists into shorter, more easily interpretable lists of biological pathways. It transforms molecular data into biological insights.

Steps in Pathway Enrichment Analysis

Input Omics Data: This can include RNA-seq, GWAS, scRNA-seq, or Epigenomics data.
Choose Reference Annotation Database: Selecting a source of biological knowledge, such as Reactome, PANTHER, or the Plant Metabolic Network (PMN).
Select Methods: Choosing the appropriate statistical approach for the gene lists of interest, such as GSEA (Gene Set Enrichment Analysis), ORA (Overrepresentation Analysis), IDEA, GSA-SNP, RSS-E, GOmeth, or active-pathway methods.
Standardize Reporting: Documenting the process accurately, including the input gene list, background set selection, reference database version, statistical tests, and p-value adjustment methods.

Gene Sets and Databases

Gene sets are groups of related genes categorized by biological states or processes. Key resources include:

Gene Ontology (GO): A framework for the model of biology, covering molecular functions, cellular components, and biological processes.
Kyoto Encyclopedia of Genes and Genomes (KEGG): A database for understanding high-level functions and utilities of the biological system.
Reactome: An open-source, curated pathway database.
WikiPathways: A community-curated database of biological pathways.
Molecular Signatures Database (MSigDB): A collection of tens of thousands of annotated gene sets for use with GSEA software.

MSigDB Collections

MSigDB is divided into Human and Mouse collections. In the Human collection, gene sets are categorized as follows:

H: Hallmark gene sets representing well-defined biological states.
C1: Positional gene sets corresponding to cytogenetic bands.
C2: Curated gene sets from pathway databases, publications, and expert knowledge.
C3: Regulatory target gene sets (microRNA seeds and transcription factor binding sites).
C4: Computational gene sets defined by mining cancer expression data.
C5: Ontology gene sets (GO terms).
C6: Oncogenic signature gene sets.
C7: Immunologic signature gene sets.
C8: Cell type signature gene sets from single-cell sequencing studies.

Overrepresentation Analysis (ORA)

ORA is a statistical method that compares the number of genes in a user-provided list (usually DEGs) that belong to a specific pathway against the total number of genes in that pathway across the whole genome.

Methodology

ORA typically employs Fisher's Exact Test to determine if a pathway is significantly overrepresented. This involves a contingency table:

Differentially Expressed (DE) & In Pathway: Gene is in the DE list and the pathway of interest.
Not DE & In Pathway: Gene is in the pathway but not in the DE list.
DE & Not in Pathway: Gene is in the DE list but not in the pathway.
Not DE & Not in Pathway: Gene is neither in the DE list nor in the pathway.

Example Scenario (IL-6 Production):

Total genes in IL-6 production pathway (GO:0032635): $15$
Total genes in Gene list: $15$
DEGs in IL-6 production: $9$
DEGs NOT in IL-6 production: $6$
Non-DEGs in IL-6 production: $3$
Non-DEGs NOT in IL-6 production: $12$
Total Genes analyzed: $30$

In this example, ORA asks whether the gene list is enriched with genes involved in IL-6 production by evaluating if $9$ out of $15$ DEGs being in that pathway is statistically more likely than by chance.

Limitations of ORA

It is heavily dependent on the chosen thresholds for differential expression (e.g., $p\text{-value} < 0.05$ and $\text{Fold Change} > 2$ ). Changing these thresholds can drastically alter the results.
It ignores the magnitude of gene expression changes beyond the binary classification of "significant" or "not significant."
It treats all genes as equally independent.

Gene Set Enrichment Analysis (GSEA)

GSEA is a computational method used to determine whether a predefined set of genes shows statistically significant, concordant differences between two biological states (e.g., healthy vs. alcoholic hepatitis).

Differences from ORA

Input: Unlike ORA, which uses a list of selected DEGs, GSEA uses a Ranked List of all genes from the experiment.
Sensitivity: GSEA can identify pathways where genes show small but consistent changes across the whole set, even if individual genes do not meet the significance threshold for ORA.

Gene Ranking Metrics

Common metrics used to rank genes for GSEA include:

Signal-to-noise ratio.
$\text{sign}(\log_{2}\text{FC}) \times (-\log_{10}(p\text{-value}))$
$\log_{2}\text{FC}$
$\log_{2}\text{FC} \times (-\log_{10}(p\text{-value}))$

Enrichment Scores and Interpretation

Running Enrichment Score (ES): This reflects the degree to which a gene set is overrepresented at the top or bottom of the ranked list. The score increases when a gene in the ranked list belongs to the pathway and decreases when it does not.
Positive ES: Genes in the pathway are mostly at the top of the ranked list (upregulated in the experimental condition).
Negative ES: Genes in the pathway are mostly at the bottom of the ranked list (downregulated in the experimental condition).
Normalised Enrichment Score (NES): The ES adjusted for the size of the gene set, allowing for comparison across different pathways.
Leading Edge: The subset of genes that appear in the ranked list before (for positive ES) or after (for negative ES) the peak of the Running ES curve. These genes contribute most significantly to the enrichment score.

Multiple Testing Correction

In both ORA and GSEA, testing thousands of pathways simultaneously increases the risk of false positives. Multiple testing correction is required to adjust p-values (e.g., FDR, Bonferroni). Increasing the number of gene sets tested results in higher p-value thresholds required for significance.

Results Representation and Visualization

Interpreting standardized outputs (which include columns like setSize, enrichmentScore, NES, pvalue, p.adjust, and qvalue) can be complex. Visualization tools help identify broader biological themes.

aPEAR (advanced Pathway Enrichment Analysis Representation): This tool simplifies results by detecting clusters of similar pathways and visualizing them as an enrichment network. This reduces redundancy and highlights the most important biological themes in the data.

Questions & Discussion

Question: What's the biological meaning of all these differences in gene expression? Response: The biological meaning is derived by performing pathway enrichment analysis. This process summarizes the statistical table of differentially expressed genes into a list of biological processes (e.g., Fibrosis, Phagocytosis, TNF signaling, Apoptosis) that are significantly altered between the conditions studied, thereby revealing the underlying biological mechanisms. For instance, the OGLDL gene belongs to the oxoglutarate dehydrogenase pathway, while NES (nestin) relates to cellular structure, and their differential expression points to specific metabolic or structural changes in the samples.