DESeq2

Prérequis

Cette section prolonge la normalisation vue en 6.3 et le comptage par gène de 6.2. Elle suppose acquis les formats de séquençage et la lecture d'une matrice de comptages, présentés en 1.2, ainsi que les notions de moyenne, variance et logarithme. Le vocabulaire des modèles linéaires généralisés est introduit ici ; aucune connaissance préalable de R n'est nécessaire, mais les commandes doivent être lues ligne par ligne.

Introduction

Une expérience de RNA-Seq compare des populations de transcrits entre conditions. La question biologique est simple : quels gènes changent, dans quel sens, et de combien ? La réponse ne peut pas venir d'une comparaison directe des comptages. Un échantillon séquencé plus profondément produit plus de lectures pour tous les gènes ; un gène plus long capte plus de fragments ; quelques transcrits très abondants peuvent déplacer l'échelle de toute la librairie. Il faut donc un modèle qui sépare l'effet biologique de ces effets techniques, puis qui quantifie l'incertitude.

L'enjeu computationnel est double. D'abord, les comptages sont des nombres entiers, positifs, dont la variance augmente avec la moyenne. Un modèle gaussien les décrit mal, surtout aux faibles valeurs. Ensuite, une expérience teste des milliers de gènes : sans contrôle du taux de fausses découvertes, les p-values brutes produisent des centaines de faux positifs. DESeq2 rassemble une normalisation par size factors, un modèle de binomiale négative, une estimation de dispersion par gène avec shrinkage, un test de Wald et une correction de Benjamini-Hochberg.

Cette section construit chaque pièce, puis exécute l'analyse complète sur un jeu synthétique. L'objectif n'est pas de mémoriser des commandes, mais de savoir ce que chaque sortie signifie, quand elle est fiable, et pourquoi les comptages bruts restent obligatoires.

Théorie

Du comptage au modèle

Pour le gène $i$ et l'échantillon $j$, notons $c_{ij}$ le nombre de lectures assignées. Ce nombre dépend de l'abondance du transcrit, de la profondeur de séquençage et de la composition de la librairie. DESeq2 écrit :

où $s_j$ est le size factor de l'échantillon, $q_{ij}$ une concentration attendue et $x_j$ le vecteur de covariables du plan expérimental. Pour un contraste simple, $x_j$ contient un intercept et un indicateur de condition ; $\beta_{i1}$ est alors le log2 fold change entre la condition traitée et la référence. Le modèle est un GLM : la moyenne dépend des covariables par une fonction de lien, et la variance dépend de la moyenne.

La normalisation ne remplace pas les comptages. Elle fournit un offset $s_j$ dans le modèle. C'est une différence importante avec une analyse qui diviserait les comptages puis appliquerait un test ordinaire : la variabilité de comptage reste décrite par la loi du modèle.

Moyenne, variance et binomiale négative

Pour des lectures réparties au hasard, la loi de Poisson impose une variance égale à la moyenne. Les données réelles sont presque toujours plus variables : les réplicats biologiques diffèrent, la préparation de librairie introduit du bruit, et l'expression d'un gène varie entre individus. On parle de surdispersion. La binomiale négative ajoute un paramètre de dispersion $\alpha_i$ :

Quand $\alpha_i$ tend vers zéro, on retrouve Poisson. Quand il augmente, la variance croît plus vite que la moyenne. Cette relation moyenne-variance est le cœur de DESeq2 : un gène fortement exprimé peut avoir une variance élevée sans être nécessairement variable biologiquement, et un gène faiblement exprimé peut être très instable malgré une moyenne petite.

Le paramètre $\alpha_i$ n'est pas connu. L'estimer sur trois réplicats est difficile, et une estimation naïve produit souvent des valeurs extrêmes. C'est pourquoi DESeq2 ne se contente pas d'une estimation par gène.

Dispersion : par gène, tendance et shrinkage

DESeq2 estime d'abord une dispersion par gène par maximum de vraisemblance, avec une correction de Cox-Reid qui tient compte de l'incertitude des coefficients. Cette dispersion brute est bruitée, surtout avec peu de réplicats. Le logiciel ajuste ensuite une tendance de dispersion en fonction de la moyenne, puis rapproche chaque estimation de cette tendance par une approche bayésienne empirique. Le résultat est une dispersion finale, dite MAP (maximum a posteriori).

Le shrinkage n'est pas un lissage cosmétique. Il réduit la variance des estimations et évite qu'un gène à dispersion extrême domine les tests. Les gènes dont la dispersion dépasse nettement la tendance peuvent être signalés comme outliers et conserver leur estimation brute, car un excès de variabilité peut refléter une biologie réelle ou un problème de qualité. Le diagnostic plotDispEsts(dds) montre les points, la tendance et les outliers.

Le choix de la tendance est paramétrable : parametric, local ou mean. La valeur par défaut convient à la plupart des jeux, mais un nombre de gènes très faible ou une relation irrégulière peut justifier local. Ce choix modifie les dispersions, donc les p-values ; il doit être rapporté.

Size factors et offset

Le size factor corrige l'échelle relative d'un échantillon sans supposer que tous les gènes changent. Pour chaque gène positif dans tous les échantillons, on calcule la moyenne géométrique :

Pour chaque échantillon $j$, on prend la médiane des rapports :

Un gène nul dans un échantillon est exclu, car son rapport est nul ou infini. La médiane résiste à une minorité de gènes fortement régulés. Elle suppose que la majorité des gènes restent stables, ou que les changements se compensent. Les facteurs ne mesurent pas la profondeur totale : deux librairies de même taille peuvent avoir des size factors différents si leur composition diffère. DESeq2 calcule cette médiane dans l'espace logarithmique pour limiter les débordements numériques ; avec un nombre impair de gènes, les deux formulations coïncident.

Prenons trois gènes et deux échantillons. Le gène A vaut 20 et 40, le gène B 10 et 20, le gène C 30 et 30. Les moyennes géométriques valent $g_A=28{,}28$, $g_B=14{,}14$ et $g_C=30$. Les rapports de S1 valent 0,707, 0,707 et 1, de médiane 0,707. Ceux de S2 valent 1,414, 1,414 et 1, de médiane 1,414. S2 est donc environ deux fois plus « profond » au sens de la médiane. Après division, les comptages de A deviennent 28,3 et 28,3 : la hausse apparente disparaît.

GLM, contrastes et test de Wald

Le modèle ajuste simultanément tous les gènes, mais chaque gène possède ses propres coefficients et sa propre dispersion. Pour un plan à deux groupes, la matrice de design contient une colonne d'intercept et une colonne indicatrice. Avec le codage de R, le premier niveau du facteur est la référence ; condition_treated_vs_control correspond au contraste traité moins contrôle.

Le test de Wald compare l'estimation du coefficient à zéro :

où $\Phi$ est la fonction de répartition de la loi normale centrée réduite. DESeq2 utilise par défaut une approximation normale ; l'option useT=TRUE permet une approximation de Student, plus prudente avec peu de réplicats. La statistique est signée : un signe positif indique une hausse dans la condition alternative.

Avec $\hat\beta=1{,}5$ et $\mathrm{SE}=0{,}5$, on obtient $W=3$. La p-value bilatérale vaut environ $0{,}0027$. Ce n'est pas la probabilité que le gène soit différentiel : c'est la probabilité d'observer une statistique au moins aussi extrême si l'effet vrai est nul, sous le modèle.

Un contraste peut comparer deux niveaux d'un facteur, une interaction ou une combinaison linéaire. Dans tous les cas, vérifiez l'ordre des niveaux et le nom du coefficient avec resultsNames(dds). Une erreur de référence inverse le signe du log2 fold change sans modifier la p-value.

Filtrage indépendant et FDR

Tester des milliers de gènes pose un problème de comparaisons multiples. Si l'on teste 10 000 gènes au seuil de 5 %, environ 500 p-values seront inférieures à 0,05 même sans effet réel. La correction de Benjamini-Hochberg contrôle le taux de fausses découvertes (false discovery rate, FDR), c'est-à-dire la proportion attendue de faux positifs parmi les gènes déclarés significatifs. On trie les p-values $p_{(1)} \le \dots \le p_{(n)}$ et on calcule :

puis on borne le résultat à 1. La correction est monotone : une p-value ajustée ne peut pas diminuer quand la p-value brute augmente.

Prenons cinq p-values : 0,001 ; 0,008 ; 0,039 ; 0,041 ; 0,2. Les valeurs ajustées valent 0,005 ; 0,020 ; 0,051 ; 0,051 ; 0,2. La troisième p-value brute est inférieure à 0,05, mais sa valeur ajustée ne l'est pas : elle n'est pas déclarée significative au FDR de 5 %. Le rang et le nombre de tests comptent autant que la p-value brute.

DESeq2 ajoute un filtrage indépendant (independentFiltering=TRUE par défaut). Les gènes à très faible moyenne normalisée ont peu de chances d'être détectés ; le logiciel choisit un seuil de moyenne qui maximise le nombre de découvertes pour un FDR cible. Les gènes éliminés reçoivent padj = NA, et non une p-value de 1. Il faut distinguer un gène non testé d'un gène testé et non significatif.

Résultats, seuils et visualisations

results() renvoie notamment baseMean, log2FoldChange, lfcSE, stat, pvalue et padj. baseMean est la moyenne des comptages normalisés. lfcSE est l'erreur standard du log2 fold change. stat est la statistique de Wald. Une table triée par padj croissant est un point de départ, pas une conclusion.

Les seuils usuels combinent une significance et un effet : padj < 0,05 et |log2FoldChange| >= 1. Ils doivent être annoncés avant l'analyse, ou justifiés après. Un gène à fort effet mais faible expression peut avoir une p-value élevée ; un gène très significatif peut avoir un effet biologiquement négligeable. Le seuil de deux fois est une convention, pas une loi.

lfcShrink() fournit des log2 fold changes rétrécis vers zéro, utiles pour classer les gènes et dessiner un volcano plot. Le shrinkage réduit les effets des gènes faiblement exprimés ou très variables. Il ne remplace pas le test : results() reste la table des p-values, et lfcShrink() doit être utilisé avec le même coefficient. Mélanger une p-value non rétrécie et un effet rétréci est fréquent et rarement documenté.

Le volcano plot place le log2 fold change en abscisse et $-\log_{10}(\mathrm{padj})$ en ordonnée. Les seuils forment une croix ; les gènes intéressants sont dans les coins supérieurs. Un axe tronqué, un seuil non indiqué ou des étiquettes trop nombreuses rendent la figure trompeuse. La heatmap représente les gènes différentiels sur une échelle de couleur, souvent après transformation stabilisante et centrage par gène. Elle sert à voir des groupes d'échantillons et des blocs de gènes, pas à prouver une significativité.

La transformation stabilisante de la variance (VST) et le rlog produisent des valeurs dont la variance dépend moins de la moyenne. VST est rapide sur de grands jeux ; rlog est plus coûteux mais utile sur de petits effectifs. Pour un contrôle qualité ou une classification non supervisée, blind=TRUE évite d'utiliser le plan expérimental. Pour une heatmap de gènes différentiels, blind=FALSE est cohérent avec l'analyse. Les comptages bruts ne doivent jamais être colorés directement : l'échelle serait dominée par quelques gènes.

Réplicats, puissance et limites

La puissance d'une analyse différentielle dépend de l'effet réel, de la dispersion, de la profondeur et du nombre de réplicats. Avec deux réplicats par groupe, la dispersion est estimée avec très peu de degrés de liberté ; le shrinkage aide, mais ne crée pas d'information. Trois réplicats constituent un minimum prudent pour des cellules en culture, davantage pour des organismes entiers ou des populations hétérogènes.

La composition est une limite structurelle. Si un petit nombre de gènes monopolise les lectures et change massivement, la médiane des rapports peut être biaisée, et une répression apparente peut toucher des gènes inchangés. Les spike-in ou une référence externe apportent une correction, mais restent des hypothèses. Un facteur de taille n'est pas une mesure absolue.

Le batch est l'autre piège. Si la condition est confondue avec la série de séquençage, aucun modèle ne peut séparer les deux effets. Ajoutez le batch au design (~ batch + condition) quand il est identifiable, et vérifiez la colinéarité. Les valeurs extrêmes, les librairies de mauvaise qualité et les gènes chevauchants doivent être examinés avant l'interprétation.

Enfin, le modèle suppose des comptages indépendants et une dispersion correctement spécifiée. Il n'est pas conçu pour des TPM, des données de puces, des proportions ou des comptages de cellules uniques sans adaptation. Une analyse DESeq2 valide reste une analyse conditionnelle au plan, aux annotations et aux seuils choisis.

Implémentation

Les commandes ci-dessous utilisent R 4.4 et DESeq2 1.38 ou ultérieur. Le fichier counts.tsv contient une colonne gene_name et quatre colonnes d'échantillons. Il faut donc séparer les annotations de la matrice numérique, puis construire un tableau de métadonnées dont les noms de lignes correspondent aux colonnes. Le code complet est commenté ; les sorties R ne sont pas exécutées dans cet environnement et sont donc illustratives.

# 1. Lire les comptages bruts et séparer annotations et matrice
raw <- read.delim("datasets/rnaseq/counts.tsv", check.names = FALSE)
gene_names <- raw$gene_name
count_matrix <- as.matrix(raw[, c("S1", "S2", "S3", "S4")])
rownames(count_matrix) <- raw$gene_id
storage.mode(count_matrix) <- "integer"
 
# 2. Construire les métadonnées : l'ordre des lignes doit suivre les colonnes
coldata <- data.frame(
  sample = c("S1", "S2", "S3", "S4"),
  condition = factor(
    c("control", "control", "treated", "treated"),
    levels = c("control", "treated")
  )
)
rownames(coldata) <- coldata$sample
stopifnot(identical(colnames(count_matrix), rownames(coldata)))
 
# 3. Construire l'objet DESeq2 et filtrer les gènes trop faibles
library(DESeq2)
dds <- DESeqDataSetFromMatrix(
  countData = count_matrix,
  colData = coldata,
  design = ~ condition
)
smallest_group <- min(table(coldata$condition))
keep <- rowSums(counts(dds) >= 10) >= smallest_group
dds <- dds[keep, ]
 
# 4. Estimer size factors, dispersions, coefficients et tests de Wald
dds <- DESeq(dds)
resultsNames(dds)
# Exemple de sortie illustrative :
# [1] "Intercept"                  "condition_treated_vs_control"
 
# 5. Table de résultats pour le contraste traité - contrôle
res <- results(
  dds,
  contrast = c("condition", "treated", "control"),
  alpha = 0.05
)
summary(res)
head(res[order(res$padj), ])
 
# 6. Log2 fold changes rétrécis pour le classement et les figures
coef_name <- resultsNames(dds)[2]
res_shrunk <- lfcShrink(dds, coef = coef_name, type = "apeglm")
head(res_shrunk[order(res_shrunk$padj), ])
 
# 7. Diagnostic de la relation moyenne-dispersion
plotDispEsts(dds, main = "Dispersion estimée et tendance")
 
# 8. MA plot : moyenne normalisée contre log2 fold change rétréci
plotMA(res_shrunk, ylim = c(-3, 3), main = "MA plot")
 
# 9. Volcano plot avec seuils explicites
library(ggplot2)
volcano_data <- as.data.frame(res_shrunk)
volcano_data$gene <- rownames(volcano_data)
volcano_data$significant <- with(
  volcano_data,
  !is.na(padj) & padj < 0.05 & abs(log2FoldChange) >= 1
)
ggplot(volcano_data, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = significant), size = 2, na.rm = TRUE) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  scale_color_manual(values = c("grey70", "#b2182b")) +
  labs(x = "log2 fold change (shrunken)", y = "-log10(padj)") +
  theme_minimal()
 
# 10. Heatmap des gènes les plus significatifs après VST
library(pheatmap)
vsd <- vst(dds, blind = FALSE)
mat <- assay(vsd)
top <- head(order(res$padj, na.last = NA), 30)
if (length(top) >= 2) {
  mat_top <- mat[top, , drop = FALSE]
  mat_scaled <- t(scale(t(mat_top)))
  annotation_col <- data.frame(condition = coldata$condition)
  rownames(annotation_col) <- colnames(mat_top)
  pheatmap(
    mat_scaled,
    annotation_col = annotation_col,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    show_rownames = TRUE,
    show_colnames = TRUE
  )
}

Trois points critiques. La matrice doit rester en nombres entiers : storage.mode(count_matrix) <- "integer" évite un échec de DESeq() avec les versions récentes. Le contraste c("condition", "treated", "control") fixe le sens du log2 fold change ; sans lui, R utilise le dernier niveau et la référence peut surprendre. Enfin, padj = NA après filtrage indépendant n'est pas une p-value de 1 : un gène non testé ne doit pas être compté comme non significatif dans un bilan naïf.

Le coût computationnel est dominé par l'ajustement des GLM et l'estimation des dispersions, pas par la normalisation. Pour $n$ gènes et $m$ échantillons, les opérations de taille sont en $O(nm)$ ; l'optimisation par gène est itérative. Sur quelques milliers de gènes, l'analyse prend quelques secondes à quelques minutes.

Étude de cas

Nous utilisons datasets/rnaseq/counts.tsv, un jeu synthétique créé pour ce cours. Il contient quatre gènes et quatre échantillons, avec deux réplicats par condition : S1 et S2 sont des contrôles, S3 et S4 des traités. Les comptages sont volontairement petits pour permettre le calcul à la main ; ils ne représentent aucune expérience réelle et ne doivent pas servir à une conclusion biologique.

La question est de savoir quels gènes changent entre contrôle et traité. Les étapes attendues sont : lire les comptages, construire les métadonnées, vérifier l'ordre des échantillons, filtrer les gènes trop faibles, estimer les size factors, ajuster le modèle, extraire le contraste traité moins contrôle, rétrécir les effets et produire un volcano plot et une heatmap. Le filtrage ne retire probablement aucun gène ici, car les quatre gènes ont des comptages suffisants dans au moins deux échantillons.

Les size factors restent proches de 1, car les librairies ont des totaux voisins et les gènes médians ne changent pas d'échelle. Le signal attendu est porté par deux gènes : geneA diminue dans les échantillons traités, geneB augmente, tandis que geneC et geneD restent stables. Cette lecture qualitative suffit pour vérifier le pipeline ; aucune p-value ni aucun log2 fold change n'est avancé, car deux réplicats ne permettent pas une estimation fiable de la dispersion et les sorties R ne sont pas exécutées ici.

L'interprétation doit rester prudente pour trois raisons. D'abord, geneD représente une part importante des lectures : sa stabilité est heureuse, mais un changement de sa part relative déplacerait les size factors. Ensuite, la médiane porte sur quatre gènes seulement ; elle est beaucoup moins robuste que sur un transcriptome entier. Enfin, le plan à deux réplicats est insuffisant pour distinguer un effet biologique d'une variation technique. Une analyse réelle demanderait plus de réplicats, un contrôle de qualité, une vérification du batch et une annotation complète.

Erreurs fréquentes

Quiz

Exercices

python -m pytest exercises/module-6/6.4-deseq2 -q

python -m pytest exercises/module-6/6.4-deseq2 -q

python -m pytest exercises/module-6/6.4-deseq2/projet -q

Résumé

• DESeq2 modélise les comptages bruts par une binomiale négative, avec une variance $\mu + \alpha\mu^2$ qui décrit la surdispersion.
• Le GLM relie la moyenne aux covariables par un lien logarithmique et utilise un size factor comme offset.
• La dispersion est estimée par gène, ajustée à une tendance selon la moyenne, puis rétrécie ; les outliers peuvent être conservés.
• Les size factors sont les médianes des rapports entre comptages et moyennes géométriques ; ils corrigent la composition, pas la profondeur totale.
• Le test de Wald compare le coefficient à zéro ; Benjamini-Hochberg contrôle le FDR et le filtrage indépendant produit des padj = NA.
• Les seuils combinent padj et effet ; lfcShrink() sert au classement et aux figures, pas à la décision statistique.
• VST et rlog stabilisent la variance pour les heatmaps ; les comptages bruts restent indispensables au modèle.
• Deux réplicats ne suffisent pas à une estimation fiable de la dispersion ; composition, batch et puissance limitent toute conclusion.

Références

• Anders, S. et Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biology, 11(10), R106.
• 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, Series B, 57(1), 289–300.
• Love, M. I., Huber, W. et Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550.
• Robinson, M. D. et Oshlack, A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology, 11(3), R25.
• Zhu, A., Ibrahim, J. G. et Love, M. I. (2019). Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics, 35(12), 2084–2092.
• DESeq2. Analyzing RNA-seq data with DESeq2, vignette et documentation officielle. https://bioconductor.org/packages/DESeq2