Transformez une liste de gènes différentiellement exprimés en les voies biologiques qui la sous-tendent — analyse de surreprésentation et GSEA, interprétées et visualisées
Un guide pratique de l’analyse d’enrichissement fonctionnel en R avec clusterProfiler. Partez d’un résultat d’expression différentielle DESeq2 et retrouvez les voies biologiques qui le sous-tendent de deux façons : l’analyse de surreprésentation (ORA) sur votre liste restreinte de gènes significatifs avec enrichGO, et l’analyse d’enrichissement de jeux de gènes (GSEA) sur la liste de gènes classée complète avec gseGO. Comprenez quand utiliser l’une ou l’autre, lisez le GeneRatio, le NES et les p-value ajustées en langage clair, et visualisez le résultat avec des diagrammes à points et le graphique du score cumulé GSEA. Sur le jeu de données airway.
padj < 0.05) et ne donne aucune direction ; la GSEA n’a besoin d’aucun seuil et donne un score signé (NES) — à la hausse ou à la baisse. L’accord entre elles rend un résultat convaincant.GeneRatio + p.adjust pour l’ORA, NES (signe = direction) + p.adjust pour la GSEA, et visualisez avec un diagramme à points et le graphique du score cumulé GSEA.L’expression différentielle (la leçon DESeq2) vous laisse avec une longue liste de gènes significatifs. La parcourir gène par gène est sans espoir et passe à côté de l’essentiel : la biologie se joue dans les voies, pas dans les gènes isolés. L’analyse d’enrichissement transforme la liste de gènes en une liste bien plus courte de programmes biologiques — réponse immunitaire, cycle cellulaire, organisation de la matrice extracellulaire — qui se sont déplacés ensemble.
Il existe deux stratégies standard, qui répondent à des questions légèrement différentes, d’où l’intérêt d’exécuter les deux. Ce guide utilise clusterProfiler, la boîte à outils R standard pour les deux, sur les gènes différentiellement exprimés de l’expérience airway (cellules musculaires lisses des voies respiratoires traitées à la dexaméthasone vs non traitées).
L’enrichissement prend un résultat d’expression différentielle en entrée. Nous le reconstruisons ici pour que la leçon soit autonome — c’est le même flux de travail que la leçon DESeq2, condensé.
library(airway)
library(DESeq2)
data("airway")
airway$dex <- relevel(airway$dex, ref = "untrt")
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- dds[rowSums(counts(dds)) >= 10, ]
dds <- DESeq(dds)
res <- results(dds, contrast = c("dex", "trt", "untrt"))
nrow(res) # genes tested[1] 22369sum(res$padj < 0.05, na.rm = TRUE) # significant at 5% FDR[1] 4000Environ quatre mille gènes sont significatifs — bien trop pour les interpréter directement. C’est exactement à cela que sert l’analyse d’enrichissement.
Les deux méthodes prennent des entrées différentes et répondent à des questions différentes. Exécutez les deux : l’accord entre elles (et entre les bases de données) est ce qui rend un résultat d’enrichissement fiable.
| ORA (surreprésentation) | GSEA (enrichissement de jeux de gènes) | |
|---|---|---|
| Entrée | votre liste restreinte de gènes significatifs (une liste oui/non) | tous les gènes, classés selon une statistique |
| Question | Une voie est-elle surreprésentée dans ma liste de hits, plus que par hasard ? | Les gènes d’une voie se regroupent-ils en haut ou en bas du classement ? |
| Seuil | en a besoin d’un (filtre padj / fold-change) |
aucun — utilise chaque gène |
| Direction | non (juste « enrichi ») | oui — un score signé (à la hausse ou à la baisse) |
| Optimal pour | une liste de hits forte et nette | des déplacements faibles mais coordonnés de nombreux gènes ensemble |
En bref : l’ORA demande « cette voie est-elle dans ma liste de hits ? » ; la GSEA demande « cette voie se déplace-t-elle de façon cohérente sur toute l’expérience, même à travers de nombreux petits changements ? »
L’ORA utilise un test hypergéométrique (équivalent au test exact de Fisher). Avec un univers de \(N\) gènes, \(K\) dans la voie, \(n\) gènes dans votre liste restreinte et \(k\) de ceux-ci dans la voie, la p-value est la probabilité de voir au moins \(k\) chevauchements par hasard : \[P(X \ge k) = \sum_{i=k}^{\min(n,K)} \frac{\binom{K}{i}\binom{N-K}{n-i}}{\binom{N}{n}}\]
La GSEA parcourt la liste de gènes classée en tenant un score d’enrichissement cumulé : il monte à chaque membre de la voie (pondéré par la métrique de classement du gène) et descend à chaque non-membre. Le score d’enrichissement (ES) est l’écart maximal de ce score cumulé par rapport à zéro ; le NES le normalise selon la taille du jeu de gènes pour que les jeux soient comparables, et son signe donne la direction. La significativité provient de la permutation des étiquettes de gènes un grand nombre de fois. Les deux méthodes appliquent ensuite la correction de Benjamini–Hochberg sur tous les jeux de gènes testés.
L’ORA teste votre liste restreinte de gènes significatifs pour chaque terme GO. Deux entrées : la liste de gènes (gènes significatifs) et l’univers (tous les gènes que vous avez testés) — l’univers compte, car « enrichi » signifie « plus qu’attendu compte tenu de ce qui aurait pu être détecté ».
library(clusterProfiler)
library(org.Hs.eg.db)
sig_genes <- rownames(subset(res, padj < 0.05))
universe <- rownames(res)
ego <- enrichGO(gene = sig_genes,
universe = universe,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
ont = "BP", # Biological Process
pAdjustMethod = "BH",
qvalueCutoff = 0.05,
readable = TRUE) # show gene symbols, not Ensembl IDs
head(as.data.frame(ego)[, c("Description", "GeneRatio", "p.adjust")], 8) Description GeneRatio p.adjust
GO:0001525 angiogenesis 197/3382 8.126049e-13
GO:0030198 extracellular matrix organization 130/3382 1.869361e-10
GO:0043062 extracellular structure organization 130/3382 1.869361e-10
GO:0045229 external encapsulating structure organization 130/3382 1.927672e-10
GO:0003013 circulatory system process 191/3382 3.447507e-10
GO:0003012 muscle system process 148/3382 7.136804e-10
GO:0030335 positive regulation of cell migration 188/3382 1.325191e-09
GO:2000147 positive regulation of cell motility 192/3382 3.042036e-09Lisez le résultat : GeneRatio est la fraction de vos gènes significatifs annotés qui tombent dans le terme ; p.adjust est sa significativité corrigée par FDR. Les termes enrichis sont des programmes de remodelage tissulaire et vasculaires — angiogenèse, motilité et migration cellulaires, processus des systèmes circulatoire et musculaire, et organisation de la matrice extracellulaire — une biologie cohérente pour du muscle lisse des voies respiratoires traité aux glucocorticoïdes, pas une liste aléatoire. Cette cohérence est le contrôle de bon sens.
Le diagramme à points est la figure d’enrichissement standard : une ligne par terme, le ratio de gènes en abscisse, la taille des points = le nombre de gènes, la couleur = la p-value ajustée.
library(enrichplot)
dotplot(ego, showCategory = 12) + ggplot2::ggtitle("ORA — GO biological process")
Le diagramme à points montre quels termes sont enrichis ; un réseau gènes-concepts montre les gènes qui les sous-tendent et quels termes ils partagent — utile pour repérer une poignée de gènes qui pilotent plusieurs voies à la fois.
cnetplot(ego, showCategory = 5)
La GSEA n’a besoin d’aucun seuil. Au lieu d’une liste restreinte, vous lui donnez chaque gène, classé selon une statistique — ici la statistique de Wald stat de DESeq2, qui porte à la fois la direction et l’intensité du changement. La GSEA demande alors si les gènes de chaque voie dérivent vers le haut (à la hausse) ou le bas (à la baisse) de ce classement.
# A named, decreasing-sorted vector: the ranked gene list.
ranks <- res$stat
names(ranks) <- rownames(res)
ranks <- sort(ranks[!is.na(ranks)], decreasing = TRUE)
set.seed(42) # GSEA permutes — seed for a reproducible result
gse <- gseGO(geneList = ranks,
OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL",
ont = "BP",
pAdjustMethod = "BH",
verbose = FALSE)
gse <- setReadable(gse, org.Hs.eg.db, keyType = "ENSEMBL")
gdf <- as.data.frame(gse)
gdf <- gdf[order(-abs(gdf$NES)), ]
head(gdf[, c("Description", "NES", "p.adjust")], 8) Description NES
GO:0046324 regulation of D-glucose import 2.217309
GO:0046323 D-glucose import 2.208967
GO:0031069 hair follicle morphogenesis -2.179956
GO:0048730 epidermis morphogenesis -2.166981
GO:0010827 regulation of D-glucose transmembrane transport 2.103370
GO:1904659 D-glucose transmembrane transport 2.099115
GO:0008645 hexose transmembrane transport 2.083159
GO:0003382 epithelial cell morphogenesis 2.066909
p.adjust
GO:0046324 0.0005099558
GO:0046323 0.0001763789
GO:0031069 0.0010133368
GO:0048730 0.0012759375
GO:0010827 0.0009633465
GO:1904659 0.0006022138
GO:0008645 0.0005099558
GO:0003382 0.0033167697Lisez le NES (score d’enrichissement normalisé) : sa magnitude est l’intensité, son signe la direction — positif = à la hausse dans les cellules traitées, négatif = à la baisse. Les mouvements les plus forts sont des programmes métaboliques et de transport qui évoluent à la hausse (import/transport du glucose) face à des programmes de morphogenèse épidermique et des follicules pileux qui évoluent à la baisse — une image signée et cohérente que l’ORA fondée sur un seuil ne peut pas vous donner.
La figure emblématique de la GSEA montre le score d’enrichissement cumulé au fur et à mesure que vous descendez la liste classée, avec les gènes de la voie marqués par des traits. Un pic près du sommet = une voie enrichie parmi les gènes surexprimés.
top_id <- gdf$ID[which.max(gdf$NES)] # the most up-regulated set
gseaplot2(gse, geneSetID = top_id, title = gdf$Description[gdf$ID == top_id])
GO est un vocabulaire, rarement le seul que vous voulez. Deux collections soigneusement constituées offrent des vues complémentaires de la même liste classée, et toutes deux s’exécutent à partir de packages d’annotation locaux — sans API web — ce qui les garde reproductibles :
Exécuter la GSEA sur ces collections repose sur la même machinerie de liste classée que gseGO — seule la source du jeu de gènes change.
Les jeux Hallmark sont indexés par symbole de gène, donc nous réétiquetons la liste classée des identifiants Ensembl vers les symboles, puis nous la passons à la fonction générique GSEA() avec Hallmark comme table de jeux de gènes.
library(msigdbr)
library(clusterProfiler)
sym <- mapIds(org.Hs.eg.db, names(ranks), "SYMBOL", "ENSEMBL")
rsym <- ranks[!is.na(sym)]
names(rsym) <- sym[!is.na(sym)]
rsym <- sort(rsym[!duplicated(names(rsym))], decreasing = TRUE)
hallmark <- msigdbr(species = "Homo sapiens", collection = "H") # the 50 Hallmark sets
set.seed(42)
gsea_h <- GSEA(geneList = rsym,
TERM2GENE = hallmark[, c("gs_name", "gene_symbol")],
pvalueCutoff = 0.25,
eps = 0,
verbose = FALSE)
hdf <- as.data.frame(gsea_h)
head(hdf[order(hdf$p.adjust), c("Description", "NES", "p.adjust")], 8) Description NES
HALLMARK_ADIPOGENESIS HALLMARK_ADIPOGENESIS 1.957457
HALLMARK_TNFA_SIGNALING_VIA_NFKB HALLMARK_TNFA_SIGNALING_VIA_NFKB 1.810453
HALLMARK_ANDROGEN_RESPONSE HALLMARK_ANDROGEN_RESPONSE 1.725946
HALLMARK_OXIDATIVE_PHOSPHORYLATION HALLMARK_OXIDATIVE_PHOSPHORYLATION 1.589944
HALLMARK_P53_PATHWAY HALLMARK_P53_PATHWAY -1.575926
HALLMARK_XENOBIOTIC_METABOLISM HALLMARK_XENOBIOTIC_METABOLISM 1.554954
HALLMARK_E2F_TARGETS HALLMARK_E2F_TARGETS -1.477450
HALLMARK_IL2_STAT5_SIGNALING HALLMARK_IL2_STAT5_SIGNALING 1.537419
p.adjust
HALLMARK_ADIPOGENESIS 1.209687e-06
HALLMARK_TNFA_SIGNALING_VIA_NFKB 1.101194e-04
HALLMARK_ANDROGEN_RESPONSE 6.748773e-03
HALLMARK_OXIDATIVE_PHOSPHORYLATION 6.748773e-03
HALLMARK_P53_PATHWAY 6.748773e-03
HALLMARK_XENOBIOTIC_METABOLISM 7.212627e-03
HALLMARK_E2F_TARGETS 1.690343e-02
HALLMARK_IL2_STAT5_SIGNALING 1.709731e-02Une trentaine de jeux Hallmark passent FDR < 0.25. Les mouvements les plus forts et les plus significatifs sont des programmes métaboliques et de signalisation qui montent — HALLMARK_ADIPOGENESIS, HALLMARK_TNFA_SIGNALING_VIA_NFKB, HALLMARK_ANDROGEN_RESPONSE, HALLMARK_OXIDATIVE_PHOSPHORYLATION — aux côtés de HALLMARK_P53_PATHWAY qui descend à la baisse (NES négatif). Cela colle à la biologie : la dexaméthasone est un glucocorticoïde, et un aperçu grossier de 50 jeux devrait faire ressortir exactement ces signatures métaboliques et anti-inflammatoires.
Un diagramme à points montre le même résultat d’un coup d’œil, réparti par direction :
library(enrichplot)
library(ggplot2)
dotplot(gsea_h, showCategory = 10, split = ".sign") + facet_grid(. ~ .sign)
Le diagramme à points répartit les jeux par direction — activés (NES positif, à la hausse dans les cellules traitées) d’un côté, réprimés de l’autre — avec la taille des points donnant le nombre de gènes du leading edge et la couleur la p-value ajustée. Lire les deux panneaux côte à côte est la façon la plus rapide de voir quels programmes sont montés et lesquels sont descendus.
Reactome est indexé par identifiants Entrez, donc nous réétiquetons la liste classée une fois de plus, puis nous exécutons la GSEA sur les voies Reactome avec gsePathway().
library(ReactomePA)
ent <- mapIds(org.Hs.eg.db, names(ranks), "ENTREZID", "ENSEMBL")
rent <- ranks[!is.na(ent)]
names(rent) <- ent[!is.na(ent)]
rent <- sort(rent[!duplicated(names(rent))], decreasing = TRUE)
set.seed(42)
reactome <- gsePathway(geneList = rent,
organism = "human",
minGSSize = 15,
maxGSSize = 500,
pvalueCutoff = 0.25,
eps = 0,
verbose = FALSE)
reactome <- setReadable(reactome, org.Hs.eg.db, keyType = "ENTREZID") # Entrez -> symbols
rdf <- as.data.frame(reactome)
head(rdf[order(rdf$p.adjust), c("Description", "NES", "p.adjust")], 8) Description NES p.adjust
R-HSA-186763 Downstream signal transduction 2.120125 0.003182459
R-HSA-6798695 Neutrophil degranulation 1.645348 0.003182459
R-HSA-379716 Cytosolic tRNA aminoacylation -2.173355 0.003325519
R-HSA-397014 Muscle contraction 1.805247 0.003325519
R-HSA-445355 Smooth Muscle Contraction 2.016220 0.007075624
R-HSA-400685 Sema4D in semaphorin signaling 2.002347 0.007633437
R-HSA-1483257 Phospholipid metabolism 1.702325 0.007633437
R-HSA-5627117 RHO GTPases Activate ROCKs 1.995327 0.008961273Les voies de Reactome sont plus fines, donc il renvoie bien plus de hits — ici près de deux cents à FDR < 0.25. Les voies de tête sont cliniquement lisibles et biologiquement pertinentes pour les cellules musculaires lisses des voies respiratoires : contraction musculaire et des muscles lisses, transduction du signal et dégranulation des neutrophiles qui montent, face à l’aminoacylation cytosolique des ARNt qui descend. C’est le niveau de détail que vous citeriez dans les résultats d’un article.
Le même diagramme à points réparti, maintenant sur les voies Reactome — activées versus réprimées, côte à côte :
library(enrichplot)
library(ggplot2)
dotplot(reactome, showCategory = 10, split = ".sign") + facet_grid(. ~ .sign)
Une collection grossière peut placer un programme sans rapport près du sommet par hasard, et aucune voie unique n’est jamais toute l’histoire. Faites confiance au résultat qui apparaît de façon constante — à travers GO, Hallmark et Reactome, et dans la direction que prédit la biologie. Trois vocabulaires qui s’accordent, voilà ce qui rend un résultat d’enrichissement convaincant.
L’ORA ne trouve rien, la GSEA trouve beaucoup (ou inversement). Elles répondent à des questions différentes, donc c’est normal. L’ORA a besoin d’une liste restreinte forte et nette ; la GSEA détecte des déplacements faibles mais coordonnés qui ne passent jamais un seuil par gène. Le désaccord est une information, pas un bug — rapportez les deux.
« No gene can be mapped » / résultat vide. Le keyType doit correspondre à vos identifiants de gènes. airway utilise des identifiants Ensembl, donc keyType = "ENSEMBL". Si vos gènes sont des symboles, utilisez "SYMBOL", s’ils sont Entrez utilisez "ENTREZID" — ou convertissez d’abord avec bitr().
Vous avez oublié l’univers dans l’ORA. Sans universe = (tous les gènes testés), enrichGO compare votre liste au génome entier, gonflant la significativité. Passez toujours les gènes qui auraient pu être détectés dans votre expérience.
L’ORA teste une liste restreinte de gènes significatifs (un seuil oui/non) pour la surreprésentation d’une voie à l’aide d’un test hypergéométrique — sans direction. La GSEA utilise tous les gènes classés selon une statistique et demande si une voie se regroupe en haut ou en bas, donnant un NES signé (à la hausse ou à la baisse) et sans besoin d’un seuil. Utilisez les deux ; l’accord est convaincant.
Cette leçon exécute déjà Reactome (gsePathway) et MSigDB Hallmark (GSEA + msigdbr) ci-dessus — tous deux à partir d’annotations locales, entièrement hors ligne. KEGG suit le même flux de travail (enrichKEGG() / gseKEGG()) avec un hic : les fonctions KEGG ont besoin d’identifiants Entrez (convertissez avec bitr()) et interrogent l’API web KEGG au moment de l’exécution, donc elles nécessitent une connexion internet — ce qui explique pourquoi cette leçon reproductible hors ligne n’exécute pas de cellule KEGG en direct. Pour l’ORA sur une liste restreinte, ReactomePA::enrichPathway() est l’équivalent de surreprésentation de la GSEA gsePathway() montrée ci-dessus.
Faites correspondre keyType à vos données. DESeq2 sur airway donne des identifiants Ensembl (keyType = "ENSEMBL"). De nombreux outils attendent des identifiants Entrez (KEGG, MSigDB), et les gens lisent les symboles. Convertissez entre eux avec clusterProfiler::bitr() ou AnnotationDbi::mapIds(), et définissez readable = TRUE pour que les résultats montrent les symboles.
C’est le score d’enrichissement GSEA normalisé selon la taille du jeu de gènes pour que les jeux soient comparables. Sa magnitude est l’intensité de l’enrichissement ; son signe est la direction — positif signifie que la voie est enrichie parmi les gènes surexprimés, négatif parmi les sous-exprimés. Filtrez sur p.adjust < 0.05, puis classez par NES.
Oui. La GSEA estime la significativité par permutation, donc les résultats varient légèrement d’une exécution à l’autre. Appelez set.seed() avant gseGO()/GSEA() pour un résultat reproductible — sinon les p-value (et parfois les termes limites) changent entre les exécutions.
L’ORA sur toute la liste significative mélange gènes induits et réprimés. Séparez-les : exécutez enrichGO sur les gènes avec log2FoldChange > 0 puis de nouveau sur ceux avec log2FoldChange < 0, et comparez les termes de tête. Quels programmes biologiques sont induits par la dexaméthasone versus réprimés ?
Construisez deux vecteurs de gènes à partir de res — rownames(subset(res, padj < 0.05 & log2FoldChange > 0)) et la version < 0 — et appelez enrichGO sur chacun, en gardant le même universe.
up <- rownames(subset(res, padj < 0.05 & log2FoldChange > 0))
down <- rownames(subset(res, padj < 0.05 & log2FoldChange < 0))
ego_up <- enrichGO(up, universe = rownames(res), OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL", ont = "BP", readable = TRUE)
ego_down <- enrichGO(down, universe = rownames(res), OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL", ont = "BP", readable = TRUE)
head(ego_up$Description, 5)
head(ego_down$Description, 5)Séparer par direction est souvent plus informatif qu’une seule liste combinée — cela vous dit ce qui est monté à la hausse versus à la baisse, ce qui est exactement la direction que la GSEA vous donne gratuitement.
A. La GSEA — elle ne regarde que les gènes significatifs. B. L’ORA — elle teste une liste restreinte définie par un seuil padj/fold-change. C. Aucune — les deux utilisent tous les gènes.
B. L’ORA teste une liste restreinte que vous définissez avec un seuil (p. ex. padj < 0.05). La GSEA utilise la liste de gènes classée complète et n’a besoin d’aucun seuil — c’est son principal avantage pour détecter de petits déplacements coordonnés.
Vous avez transformé un résultat d’expression différentielle de quatre mille gènes en une courte liste de voies biologiques, de deux façons complémentaires : l’ORA (enrichGO) a testé la liste restreinte de gènes significatifs et fait ressortir les programmes de matrice extracellulaire, d’angiogenèse et musculaires ; la GSEA (gseGO) a utilisé la liste classée complète pour donner une image signée de ce qui est monté versus descendu. Vous avez ensuite exécuté la même GSEA à travers trois vocabulaires — GO, MSigDB Hallmark et Reactome — tous hors ligne, et vu les signatures du muscle lisse des voies respiratoires et des glucocorticoïdes se répéter d’un vocabulaire à l’autre. Les règles pratiques : passez toujours l’univers à l’ORA, classez selon une statistique signée (la statistique de Wald stat) pour la GSEA, fixez une graine aléatoire, faites correspondre le keyType à vos identifiants (ou convertissez avec bitr()), et lisez GeneRatio/NES + p.adjust — puis faites confiance aux résultats qui s’accordent entre méthodes et bases de données.
enrichKEGG/gseKEGG) mais nécessite une API web — voir la FAQ.Cette leçon est reproductible : chaque tableau et chaque figure a été produit par le code montré — copiez n’importe quel bloc et exécutez-le pour les reproduire. The runtime is the judge.
@online{2026,
author = {},
title = {Analyse d’enrichissement de jeux de gènes et de voies en R
avec clusterProfiler},
date = {2026-06-29},
url = {https://www.datanovia.com/learn/bioinformatics/bulk-rna-seq/gene-set-enrichment-analysis},
langid = {fr}
}