Enrichissement GO

Prérequis

Cette section prolonge l'analyse différentielle de 6.4 : la liste de gènes qui en sort devient la liste d'intérêt. Elle suppose connue la lecture des identifiants et des annotations transférées par homologie, abordée en 2.6. Les principes de script versionné et de sortie reproductible de 0.1 s'appliquent tels quels. Les exemples numériques sont synthétiques et les commandes R sont fournies à titre illustratif.

Introduction

Une analyse différentielle transcriptomique se termine rarement par un tableau de gènes : elle pose une question biologique. Quels processus cellulaires distinguent une bactérie soumise à un stress d'une bactérie témoin ? Quelles fonctions sont communes aux gènes surexprimés ? La Gene Ontology (GO) fournit un vocabulaire contrôlé pour décrire ces fonctions, et l'enrichissement GO demande si une liste de gènes contient plus de gènes annotés à un terme que ce que le hasard produirait.

L'enjeu computationnel est double. Il faut d'abord construire une table de contingence honnête : la liste d'intérêt, l'univers de référence, la taille du terme et le chevauchement. Il faut ensuite tester des milliers de termes, corriger la multiplicité, puis démêler les termes redondants et les effets de taille. Un logiciel peut produire une liste de p-values en quelques secondes ; il ne peut pas décider à votre place si le résultat reflète une biologie ou un artefact d'annotation.

Nous partons d'une petite table calculable à la main, puis nous passons aux fonctions Python, à clusterProfiler et à Enrichr. L'objectif est de savoir refaire le calcul et de reconnaître les choix qui changent le résultat.

Théorie

Une ontologie, trois vocabulaires

La Gene Ontology est un vocabulaire structuré, organisé en trois ontologies indépendantes. Biological Process (BP) décrit les processus auxquels un gène participe, Cellular Component (CC) sa localisation, et Molecular Function (MF) son activité moléculaire. Un même gène peut recevoir plusieurs termes dans chaque ontologie.

Les termes ne forment pas une arborescence mais un graphe orienté acyclique (directed acyclic graph, DAG). Un terme enfant est plus spécifique que ses parents, et les relations incluent is_a et part_of. La règle du chemin vrai (true path rule) veut qu'un gène annoté à un terme soit implicitement annoté à ses ancêtres. Cette propriété crée une dépendance forte entre termes : si un terme enfant est enrichi, ses parents ont souvent tendance à l'être aussi.

Une annotation GO associe un produit de gène à un terme, avec une source et un code de preuve. Le code distingue preuve expérimentale, inférence de similarité, annotation électronique et curation manuelle. Une annotation n'est pas une mesure d'expression : elle décrit ce que la communauté sait ou infère sur la fonction d'un gène, à une version donnée de la base.

Liste d'intérêt et univers de référence

La liste d'intérêt, ou premier plan (foreground), contient les gènes que l'on veut interpréter : gènes différentiels, gènes d'un cluster, cibles d'un facteur de transcription. Sa taille est notée $n$. L'univers de référence, ou arrière-plan (background), contient tous les gènes qui auraient pu figurer dans cette liste. Sa taille est notée $N$. Le terme étudié contient $K$ gènes dans cet univers, et $k$ gènes de la liste appartiennent au terme.

Pour un RNA-Seq, l'univers naturel n'est pas le génome entier mais l'ensemble des gènes réellement testés : gènes exprimés, gènes passant le filtre de counts, gènes annotés dans le type d'identifiant utilisé. Si la liste a été sélectionnée parmi 2 500 gènes exprimés, prendre le génome complet ne correspond plus au modèle d'échantillonnage. Les p-values ne sont plus calibrées, et deux analyses menées avec des univers différents ne sont pas comparables.

L'univers doit contenir tous les gènes de la liste et tous les gènes des termes. Les tailles $K$ doivent être recalculées dans ce même univers : un terme de 300 gènes sur le génome peut n'en contenir que 120 parmi les gènes exprimés. Vérifiez aussi que les identifiants sont homogènes — symbole, Ensembl, Entrez — et que la conversion n'a pas perdu de gènes silencieusement.

Le test hypergéométrique

Sous l'hypothèse nulle, la liste est un échantillon aléatoire de $n$ gènes tirés sans remise dans l'univers, et l'appartenance au terme est fixée. Le nombre $X$ de gènes de la liste qui appartiennent au terme suit une loi hypergéométrique. La p-value unilatérale à droite est la probabilité d'observer au moins $k$ gènes du terme :

La table de contingence sous-jacente s'écrit :

L'espérance sous l'hypothèse nulle vaut $nK/N$. Le rapport entre l'observé et l'attendu, $k/(nK/N)$, est le fold enrichment. Une p-value petite signifie que le chevauchement observé est peu probable si la liste était un tirage aléatoire ; elle ne signifie pas que le terme est la cause de la sélection.

Prenons un univers jouet de 20 gènes, une liste de 5 gènes et un terme de 4 gènes, avec 3 gènes en commun. Alors

L'espérance vaut $5 \times 4 / 20 = 1$ et le fold enrichment vaut 3. Le résultat est nominalement significatif à 5 %, mais si les 19 autres termes ont des p-values élevées, la correction de Benjamini-Hochberg porte la q-value de ce terme à $q = 0{,}032 \times 20 = 0{,}64$. Le même chevauchement ne survit plus à la multiplicité.

Le calcul coûte $O(\min(K,n))$ termes binomiaux. Pour de grands effectifs, les combinaisons deviennent énormes et l'on travaille en logarithme avec la fonction gamma. Le nombre de termes testés peut atteindre plusieurs milliers : la correction et la redondance comptent autant que le calcul.

Odds ratio et taille d'effet

L'odds ratio (OR) résume l'association de la table 2 x 2 :

Il vaut 1 lorsque la proportion de gènes du terme est la même dans la liste et hors de la liste. Il est supérieur à 1 en cas de sur-représentation et inférieur à 1 en cas de déplétion. Ce n'est pas une probabilité : il peut dépasser 1 et n'a pas de borne. Dans l'exemple jouet, $\mathrm{OR} = 3 \times 14 / (2 \times 1) = 21$, alors que le fold enrichment vaut seulement 3 : les deux mesures ne racontent pas la même chose.

Un odds ratio est instable quand une cellule est petite ou nulle. Un terme sans gène dans la liste donne un OR de 0, un terme qui contient toute la liste peut donner un OR infini. Certains logiciels ajoutent 0,5 à chaque cellule (Haldane-Anscombe) ; signalez-le, car la valeur change. Avec de petits effectifs, l'intervalle de confiance de l'OR est large et l'ordre des termes est fragile.

La p-value et la taille d'effet répondent à deux questions différentes. La p-value demande si le chevauchement est compatible avec le hasard ; l'odds ratio et le fold demandent quelle est l'ampleur de l'association. Un terme très large peut être significatif avec un fold de 1,2, tandis qu'un terme précis peut avoir un fold de 5 sans atteindre le seuil après correction.

Multiplicité, FDR et Benjamini-Hochberg

Tester mille termes, c'est accepter environ cinquante faux positifs si l'on garde un seuil nominal de 5 % sans correction. La correction de Bonferroni contrôle la probabilité de commettre au moins une erreur (family-wise error rate, FWER) en testant chaque terme à $\alpha/m$ ; elle est très conservatrice, surtout avec des termes corrélés. La procédure de Benjamini-Hochberg (BH) contrôle le taux de fausses découvertes (false discovery rate, FDR), c'est-à-dire l'espérance de la proportion de faux positifs parmi les termes déclarés significatifs.

On trie les $m$ p-values par ordre croissant, $p_{(1)} \le \dots \le p_{(m)}$. On cherche le plus grand rang $i$ tel que $p_{(i)} \le i\alpha/m$, puis on déclare significatifs tous les termes jusqu'à ce rang. Les q-values se calculent en parcourant les rangs du plus grand au plus petit :

Avec cinq p-values $0{,}001$, $0{,}008$, $0{,}02$, $0{,}04$ et $0{,}5$, les q-values valent $0{,}005$, $0{,}02$, $0{,}0333$, $0{,}05$ et $0{,}5$. Au seuil de 5 %, les quatre premières passent. Une q-value de 0,04 signifie que, parmi les termes déclarés significatifs, environ 4 % sont attendus faux en moyenne ; ce n'est pas la probabilité que ce terme précis soit faux.

BH suppose des tests indépendants ou positivement dépendants, ce qui est souvent acceptable pour GO mais pas garanti. La correction ne règle pas la redondance : deux termes parents et enfants peuvent tous deux passer.

Redondance et similarité sémantique

La structure en DAG fabrique de la redondance. Si « processus glycolytique » est enrichi, « processus métabolique des glucides » a de fortes chances de l'être aussi, sans apporter une seconde découverte. De même, deux termes peuvent partager la majorité de leurs gènes annotés. Compter les termes significatifs comme autant de preuves indépendantes surestime la solidité du résultat.

Plusieurs stratégies atténuent ce problème. clusterProfiler propose simplify, qui regroupe les termes selon une similarité sémantique et conserve un représentant par groupe, avec un seuil souvent fixé à 0,7. REVIGO résume une liste de termes en catégories lisibles. Les cartes d'enrichissement (enrichment maps) regroupent les termes qui partagent des gènes. Une simple similarité de Jaccard entre ensembles de gènes, ou entre gènes chevauchants, donne un diagnostic rapide.

La similarité sémantique peut se mesurer avec l'information content (IC) du plus proche ancêtre commun, comme dans la mesure de Resnik, ou avec les mesures de Lin ou de Wang. Le choix change les regroupements ; il doit être documenté.

Taille des termes, puissance et biais

La taille $K$ d'un terme gouverne la puissance et la stabilité. Un terme de 4 gènes peut produire une p-value spectaculaire avec un seul gène de plus ou de moins ; les p-values sont discrètes et l'estimation est fragile. Un terme de 800 gènes a beaucoup de puissance : un fold de 1,2 peut être significatif, mais le terme est trop général pour être informatif. Une plage de 10 à 500 gènes est fréquente, mais doit être justifiée et fixée avant de regarder les résultats.

Les tests conditionnels, comme ceux de topGO, réduisent l'effet de la hiérarchie en testant un terme conditionnellement à ses parents. Pour un RNA-Seq, GOseq modélise le biais de longueur des gènes, qui influence la détection différentielle. Ces méthodes corrigent des biais connus, mais ne remplacent pas un univers correct.

L'annotation elle-même est biaisée. Les organismes modèles, les familles de gènes très étudiées et les voies canoniques sont mieux annotées que les gènes rares ou récemment découverts. Une annotation électronique peut propager une erreur, et l'absence d'annotation ne signifie pas l'absence de fonction. Deux versions de GO ou de GOA ne contiennent pas les mêmes termes : comparer deux analyses exige de fixer la version et de la citer. L'enrichissement reflète donc la composition de l'annotation autant que la biologie de la liste.

Implémentation

Le code Python suivant reprend les trois briques du calcul et les applique à l'exemple jouet.

Aucune sortie. Exécutez le code pour voir le résultat.

La sortie vérifiée le 2026-11-01 avec Python 3.13 est :

p = 0.031992
OR = 21.0
[0.005, 0.02, 0.03333333333333333, 0.05, 0.5]

Trois points sont critiques. Le calcul de la p-value doit sommer toute la queue, pas seulement le terme $i=k$. L'odds ratio doit être calculé sur les mêmes effectifs que la p-value, sinon les deux chiffres décrivent des tables différentes. La correction BH doit conserver l'ordre d'entrée, sinon les q-values ne correspondent plus aux termes.

En R, clusterProfiler automatise la table et la correction. L'exemple ci-dessous est une sortie illustrative : les vecteurs contiennent des identifiants Entrez, et le type de clé doit être adapté à votre annotation. Le paramètre universe est l'élément le plus important ; sans lui, le paquet utilise l'univers de l'OrgDb, qui n'est pas nécessairement celui de votre expérience.

library(clusterProfiler)
library(org.EcK12.eg.db)
 
ego <- enrichGO(
  gene = interest,
  universe = background,
  OrgDb = org.EcK12.eg.db,
  keyType = "ENTREZID",
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.2
)
 
resultats <- as.data.frame(ego)
resultats[, c("Description", "GeneRatio", "BgRatio", "pvalue", "p.adjust", "Count")]
dotplot(ego, showCategory = 15, orderBy = "GeneRatio")

dotplot place les termes sur l'axe vertical, la proportion de gènes de la liste sur l'axe horizontal et code la p-value ajustée par la couleur. La taille des points représente le nombre de gènes, pas un effet biologique. Pour réduire la redondance, simplify(ego, cutoff = 0.7, by = "p.adjust", select_fun = min) conserve un représentant par groupe sémantique.

Enrichr fonctionne autrement. On dépose une liste dans l'interface web, on choisit une bibliothèque GO, et le service renvoie les termes, les gènes chevauchants, une p-value, une p-value ajustée et un score combiné. L'interface web classique utilise l'univers de la bibliothèque et ne permet pas toujours de fournir un univers personnalisé ; vérifiez ce point dans la version utilisée. Le score combiné est un rang heuristique, pas une probabilité. Pour la reproductibilité, notez la bibliothèque exacte, sa date et l'identifiant de la liste.

Étude de cas

Exemple synthétique. L'univers compte 2 000 gènes (g0001 à g2000), la liste 80 gènes et six termes dont les tailles et chevauchements sont fixés. Aucun gène ne correspond à un organisme réel ; les résultats servent à vérifier la chaîne de calcul.

univers = [f"g{i:04d}" for i in range(1, 2001)]
interet = univers[:20] + univers[100:106] + univers[170:171] + univers[340:393]
 
termes = {
    "division cellulaire": univers[:100],
    "réparation de l'ADN": univers[100:160],
    "glycolyse": univers[160:190],
    "ribosome": univers[190:340],
    "septum": univers[:60] + univers[500:540],
    "métabolisme": univers[0:1600:2],
}

La liste contient 80 gènes. « Division cellulaire » contient 100 gènes, dont 20 dans la liste ; « septum » partage 60 gènes avec lui et retrouve les mêmes 20 gènes. « Métabolisme » est un terme large de 800 gènes, dont 41 dans la liste. Les gènes de « ribosome » sont absents de la liste.

Les deux premiers termes ne sont pas deux découvertes indépendantes : ils partagent 60 gènes et les mêmes 20 gènes de la liste ; le FDR ne le signale pas. « Métabolisme » passe le seuil avec un fold de 1,28 : la puissance d'un terme de 800 gènes rend significatif un effet modeste, à ne pas confondre avec un fold de 5.

« Réparation de l'ADN » passe de justesse ($q = 0{,}045$). Ce résultat fragile dépend de l'univers, de l'annotation et de la liste. « Glycolyse » n'est pas enrichie. « Ribosome » a un odds ratio de 0 et une p-value de 1 : le test unilatéral à droite ne teste pas la déplétion. Pour conclure à une sous-représentation, il faudrait calculer la queue inférieure ou interpréter l'odds ratio avec un intervalle de confiance.

Aucune conclusion biologique ne doit être tirée de ce jeu : un terme large peut être significatif avec un petit effet, deux termes redondants gonflent le nombre de lignes significatives, et une p-value de 1 n'est pas une absence d'association.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-6/6.5-enrichissement-go -q

python -m pytest exercises/module-6/6.5-enrichissement-go -q

python -m pytest exercises/module-6/6.5-enrichissement-go/projet -q

Résumé

• GO réunit BP, CC et MF en DAG ; un gène annoté à un terme est implicitement annoté à ses ancêtres.
• L'univers doit représenter les gènes qui pouvaient être sélectionnés ; pour un RNA-Seq, les gènes exprimés et testés sont un choix cohérent.
• Le test hypergéométrique compare le chevauchement à un tirage sans remise ; il est unilatéral et ne teste pas la déplétion.
• L'odds ratio et le fold enrichment mesurent l'ampleur de l'association ; une p-value petite ne garantit ni un grand effet ni une importance biologique.
• Benjamini-Hochberg contrôle le FDR ; une q-value de 0,04 s'interprète sur l'ensemble des termes déclarés, pas sur un terme isolé.
• Les termes GO sont redondants : simplifiez ou regroupez avant de compter les découvertes.
• La taille des termes, le biais d'annotation et la version des bases changent le résultat ; fixez-les et documentez-les.

Références

• Ashburner, M. et al. (2000). Gene Ontology: tool for the unification of biology. Nature Genetics, 25(1), 25–29.
• Benjamini, Y. et Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B, 57(1), 289–300.
• Yu, G., Wang, L.-G., Han, Y. et He, Q.-Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS, 16(5), 284–287.
• Chen, E. Y. et al. (2013). Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics, 14, 128.
• Documentation clusterProfiler. https://bioconductor.org/packages/clusterProfiler/
• Documentation Enrichr. https://maayanlab.cloud/Enrichr/help