Agrégez les cellules de chaque donneur en un seul profil, exécutez DESeq2 avec les donneurs comme réplicats biologiques, et trouvez les gènes pilotés par la condition sans les faux positifs que produisent les tests par cellule
Un guide pratique pour trouver les gènes qui changent entre conditions en RNA-seq single-cell multi-échantillons — la bonne façon. Voyez pourquoi les tests par cellule comme le FindMarkers de Seurat gonflent les faux positifs sur les comparaisons inter-conditions (pseudoréplication), puis agrégez les cellules de chaque donneur en un profil pseudobulk et exécutez DESeq2 avec les donneurs comme réplicats biologiques, sur le jeu de données PBMC à l’interféron-β de Kang (8 donneurs, contrôle vs stimulé).
Date de publication
3 juillet 2026
Modifié
7 juillet 2026
AstucePoints clés
L’expression différentielle (DE) inter-conditions en données single-cell a une seule méthode correcte, et ce n’est pas FindMarkers. Tester chaque cellule comme un réplicat indépendant (un test de Wilcoxon par cellule) traite des milliers de cellules comme des milliers d’échantillons — mais vos réplicats biologiques sont les donneurs, pas les cellules. Cette pseudoréplication gonfle le taux de fausses découvertes (FDR) : les p-values deviennent astronomiquement petites et des centaines de gènes ressortent « significatifs » même quand rien n’est réellement différent.
La solution, c’est le pseudobulk. Sommez les comptages bruts de chaque donneur au sein d’un type cellulaire en un profil par donneur et par condition, puis exécutez une méthode standard de RNA-seq bulk — DESeq2, edgeR ou limma-voom — avec les donneurs comme réplicats. C’est la bonne pratique établie du domaine (Squair et al., 2021; Crowell et al., 2020).
Les donneurs sont les réplicats, donc une expérience valide en nécessite plusieurs. Nous utilisons l’étude Kang à l’interféron-β — 8 donneurs, contrôle vs stimulé — ce qui explique exactement pourquoi elle fonctionne là où les données à 2 échantillons de la leçon d’intégration ne le pouvaient pas : le pseudobulk ne peut pas fabriquer une réplication que vous n’avez pas collectée.
Bloquez sur le donneur. Dans le plan DESeq2 ~ ind + stim, ind (le donneur) est un facteur qui absorbe les différences de base d’un donneur à l’autre, de sorte que l’effet stim est estimé au sein de chaque donneur. Oublier factor() — quand l’identifiant de donneur est un nombre — ajuste silencieusement le mauvais modèle.
Le résultat est une table propre et interprétable : le changement d’expression en log2, une p-value ajustée par FDR, et le familier volcano. Sur les monocytes CD14+, les meilleurs résultats sont les gènes stimulés par l’interféron (ISG) canoniques — CCL8, ISG20, IFITM2, IFIT3, ISG15 — surexprimés dans la condition stimulée, exactement comme la biologie le prédit.
Introduction
Vous avez intégré vos échantillons, les avez regroupés et avez nommé les types cellulaires. Vient maintenant la question que l’expérience a été bâtie pour trancher : au sein d’un type cellulaire, quels gènes ont changé avec le traitement ? Pour du sang stimulé à l’interféron-β, quels gènes les monocytes ont-ils activés en réponse ?
Le raccourci tentant, c’est de saisir l’outil que vous connaissez déjà. Sous-ensemblez les monocytes, exécutez le FindMarkers() de Seurat comparant stimulé versus contrôle, et vous obtenez une liste bien nette — souvent plus d’un millier de gènes, beaucoup avec des p-values de l’ordre de 1e-80. Cela ressemble à un résultat triomphal. C’est en grande partie un artefact.
FindMarkers() exécute un test par cellule (un test de rangs de Wilcoxon par défaut), et il traite chaque cellule comme un réplicat indépendant. Mais deux monocytes du même donneur ne sont pas des mesures indépendantes de « la réponse des monocytes » — ils partagent la génétique, le lot et le traitement de ce donneur. Vos vrais réplicats biologiques sont les donneurs. Compter les cellules comme des réplicats gonfle la taille d’échantillon effective d’une poignée de donneurs à des milliers de cellules, ce qui écrase les p-values et gonfle le FDR. C’est la pseudoréplication, et elle est désormais bien documentée comme la mauvaise façon de faire de la DE inter-conditions en données single-cell (Squair et al., 2021; Crowell et al., 2020; Murphy & Skene, 2022).
Cette leçon enseigne la méthode qui corrige cela : l’agrégation en pseudobulk. Vous condensez les cellules de chaque donneur en un profil d’expression unique, puis vous confiez ces profils à un moteur de RNA-seq bulk — ici DESeq2 — avec les donneurs comme réplicats et le donneur comme facteur de blocage. Nous travaillons sur le jeu de données PBMC à l’interféron-β Kang (cellules mononucléées du sang périphérique), qui compte huit donneurs mesurés dans les deux conditions — la réplication biologique qui rend l’ensemble valide.
Les données : huit donneurs, deux conditions
Le kang_subset.rds commité est une tranche de l’étude canonique Kang 2018 (Kang et al., 2018) : des PBMC de 8 donneurs, chacun profilé au repos (ctrl) et après stimulation à l’interféron-β (IFN-β) (stim), à travers trois types cellulaires. Elles sont stockées sous forme d’une matrice de comptage brute plus une table de métadonnées à trois colonnes — ind (le donneur), stim (la condition), et cell (le type cellulaire).
library(Matrix)d <-readRDS("_data/kang_subset.rds")counts <- d$counts # genes x cells, raw counts (a sparse dgCMatrix)meta <- d$meta # ind (donor), stim (ctrl/stim), cell (cell type)# The design: donors x conditions x cell typestable(meta$cell, meta$stim)
ctrl stim
B cells 870 888
CD14+ Monocytes 1176 1200
CD4 T cells 1200 1200
CD8 T cells 0 0
Dendritic cells 0 0
FCGR3A+ Monocytes 0 0
Megakaryocytes 0 0
NK cells 0 0
length(unique(meta$ind)) # the biological replicates
[1] 8
Trois types cellulaires, chacun présent dans les deux conditions, et — le nombre qui compte — 8 donneurs. Ces donneurs sont les réplicats biologiques. La leçon d’intégration utilisait un jeu de données à deux échantillons, ce qui suffit pour aligner les échantillons mais pas pour tester un effet de condition : avec un seul donneur par bras, vous ne pouvez pas séparer le traitement des idiosyncrasies de ce seul donneur. La DE inter-conditions requiert de la réplication, et ce jeu de données en a.
Nous exécuterons l’exemple travaillé sur les monocytes CD14+, le répondeur à l’interféron le plus fort et le plus net.
Voir le problème : le test par cellule déclare presque tout significatif
Avant de corriger quoi que ce soit, regardez ce que le raccourci fait réellement. Prenez les monocytes CD14+, log-normalisez exactement comme le fait Seurat, et exécutez un test de Wilcoxon par gène comparant stimulé versus contrôle — le même test de rangs que FindMarkers() exécute par défaut — en traitant chaque cellule comme un réplicat.
library(Matrix)d <-readRDS("_data/kang_subset.rds")counts <- d$counts; meta <- d$metamono <- meta$cell =="CD14+ Monocytes"cm <- counts[, mono]; mm <- meta[mono, ]# Log-normalize (what Seurat's LogNormalize does), keep genes expressed in >=10% of cellscs <- Matrix::colSums(cm)ln <-log1p(sweep(cm, 2, cs, "/") *1e4)ln <-as.matrix(ln[Matrix::rowMeans(cm >0) >=0.10, ])# One Wilcoxon rank-sum test per gene, ctrl vs stim, treating each CELL as a replicatepv <-apply(ln, 1, function(g) wilcox.test(g ~ mm$stim)$p.value)padj <-p.adjust(pv, method ="BH")c(cells_tested =ncol(ln),genes_tested =length(pv),significant =sum(padj <0.05),pct_significant =round(100*mean(padj <0.05), 1),at_p_below_1e50 =sum(pv <1e-50))
Le verdict est absurde à première vue : sur 1 710 gènes testés, 1 422 — soit environ 83 % — ressortent « significatifs » à un FDR de 5 %, et 254 d’entre eux se situent sous une p-value brute de 1e-50. Aucune comparaison biologique réelle ne rend 83 % du transcriptome différentiellement exprimé. Le moteur ne mesure pas la biologie ; il mesure les ~2 400 cellules que vous lui avez fournies. Avec autant de « réplicats », le moindre souffle de différence franchit n’importe quel seuil. FindMarkers() donne la même image — il exécute exactement ce test.
Prouvez que c’est un artefact : un cas nul sans différence réelle
« 83 %, c’est trop » est un argument. Voici la preuve. Construisez une comparaison où il n’existe aucune différence biologique — ne prenez que les monocytes contrôle et répartissez les 8 donneurs en deux groupes aléatoires. Tout gène qui ressort « significatif » est, par construction, un faux positif. Exécutez les deux méthodes sur ce cas nul.
library(Matrix); library(DESeq2); library(ggplot2)d <-readRDS("_data/kang_subset.rds")counts <- d$counts; meta <- d$meta# A TRUE NULL: only control monocytes, donors split into two fake groups -> no real DE existssel <- meta$cell =="CD14+ Monocytes"& meta$stim =="ctrl"cm <- counts[, sel]; mm <- meta[sel, ]set.seed(1)groupA <-sample(sort(unique(mm$ind)), 4)mm$fake <-ifelse(mm$ind %in% groupA, "A", "B")# Per-cell Wilcoxon: cells as replicatescs <- Matrix::colSums(cm)ln <-log1p(sweep(cm, 2, cs, "/") *1e4)ln <-as.matrix(ln[Matrix::rowMeans(cm >0) >=0.10, ])wp <-p.adjust(apply(ln, 1, function(g) wilcox.test(g ~ mm$fake)$p.value), "BH")# Pseudobulk DESeq2: donors as replicates (one profile per donor)pb <-t(rowsum(t(as.matrix(cm)), group =paste0("d", mm$ind)))si <-unique(data.frame(sample =paste0("d", mm$ind), grp = mm$fake))rownames(si) <- si$sample; si <- si[colnames(pb), ]si$grp <-factor(si$grp, levels =c("B", "A"))pb <- pb[rowSums(pb) >=10, ]dds <-DESeq(DESeqDataSetFromMatrix(round(pb), si, design =~ grp), quiet =TRUE)dp <-results(dds, contrast =c("grp", "A", "B"))$padjfp <-c("Per-cell Wilcoxon\n(cells as replicates)"=sum(wp <0.05, na.rm =TRUE),"Pseudobulk DESeq2\n(donors as replicates)"=sum(dp <0.05, na.rm =TRUE))print(fp)
Per-cell Wilcoxon\n(cells as replicates)
86
Pseudobulk DESeq2\n(donors as replicates)
0
bar <-data.frame(method =factor(names(fp), levels =names(fp)), fp =as.integer(fp))ggplot(bar, aes(method, fp, fill = method)) +geom_col(width =0.6) +geom_text(aes(label = fp), vjust =-0.3, size =5) +scale_fill_manual(values =c("#d1495b", "#3a86d4"), guide ="none") +labs(x =NULL, y ="False positives (genes at 5% FDR)",title ="A null comparison: no real difference exists") +theme_minimal(base_size =12)
Il n’y a rien à trouver — même condition, répartition aléatoire des donneurs — et pourtant le test de Wilcoxon par cellule rapporte 86 faux positifs, tandis que le pseudobulk DESeq2 en rapporte correctement 0. C’est la pseudoréplication rendue visible : le test par cellule a tort avec assurance parce qu’il a pris les cellules pour des réplicats. Aucune méthode qui fabrique 86 découvertes à partir de pur bruit ne peut être digne de confiance sur une vraie comparaison. Construisez maintenant la méthode qui, elle, l’est.
La solution, étape 1 : agréger en pseudobulk
« Pseudobulk » veut dire exactement ce que le mot suggère — vous reconstruisez un échantillon de RNA-seq bulk pour chaque donneur en sommant les comptages bruts de toutes ses cellules (au sein d’un même type cellulaire). Huit donneurs × deux conditions donnent 16 échantillons pseudobulk, chacun un profil de comptage unique, et désormais le donneur est l’unité de réplication — comme il se doit.
L’agrégation tient en une ligne de R de base : rowsum() somme les lignes selon un facteur de regroupement, donc sommer les colonnes (cellules) par donneur-condition revient à rowsum() sur la matrice transposée, retransposée ensuite.
library(Matrix)d <-readRDS("_data/kang_subset.rds")counts <- d$counts; meta <- d$metamono <- meta$cell =="CD14+ Monocytes"cm <- counts[, mono]; mm <- meta[mono, ]# One label per donor per condition, e.g. "101_ctrl", "101_stim", ...sample_id <-paste(mm$ind, mm$stim, sep ="_")# PSEUDOBULK: sum each donor-condition's cells into one raw-count profilepb <-t(rowsum(t(as.matrix(cm)), group = sample_id))dim(pb) # genes x 16 pseudobulk samples
[1] 15870 16
pb[c("ISG15", "CCL8"), 1:6] # summed counts for two ISGs, first 6 samples
La matrice est passée de ~2 400 cellules à 16 échantillons — 8 donneurs dans chaque condition. Si vous gardez vos données dans un objet Seurat, AggregateExpression(obj, group.by = c("ind", "stim")) effectue la sommation identique ; le rowsum() en R de base montre ici exactement ce que cet appel calcule.
Construisez maintenant la table d’échantillons (le colData) qui indique à quel donneur et à quelle condition appartient chaque colonne pseudobulk :
Seize lignes : chaque donneur apparaît deux fois, une fois par condition. Cette structure appariée est tout l’enjeu — chaque donneur est son propre contrôle.
La solution, étape 2 : DESeq2 avec les donneurs comme réplicats
Confiez la matrice pseudobulk à DESeq2 (Love et al., 2014), le moteur même que vous utiliseriez pour une expérience de RNA-seq bulk. La formule de plan porte toute la logique :
~ ind + stim — teste l’effet stim (la condition) après avoir bloqué sur ind (le donneur). Comme chaque donneur est mesuré dans les deux conditions, ind absorbe le niveau d’expression de base de chaque donneur, et stim est alors estimé au sein de chaque donneur. C’est un plan apparié, et il est bien plus puissant qu’en ignorant le donneur.
si$ind <- factor(si$ind) n’est pas optionnel. Les identifiants de donneur sont ici stockés comme des nombres (101, 107, 1015, …). Si vous les laissez numériques, DESeq2 traite le donneur comme une covariable continue — « l’expression augmente avec le numéro de donneur » — ce qui n’a aucun sens, et il vous avertit même (« did you mean for this to be a factor? »). Envelopper le donneur dans factor() en fait la variable de blocage qu’il doit être.
library(Matrix); library(DESeq2)d <-readRDS("_data/kang_subset.rds")counts <- d$counts; meta <- d$metamono <- meta$cell =="CD14+ Monocytes"cm <- counts[, mono]; mm <- meta[mono, ]sample_id <-paste(mm$ind, mm$stim, sep ="_")pb <-t(rowsum(t(as.matrix(cm)), group = sample_id))si <-unique(data.frame(sample = sample_id, ind = mm$ind, stim = mm$stim))rownames(si) <- si$sample; si <- si[colnames(pb), ]# CRITICAL: donor is a FACTOR (a blocking variable), never a numbersi$ind <-factor(si$ind)si$stim <-factor(si$stim, levels =c("ctrl", "stim"))# Drop near-empty genes, then fit: block on donor, test the conditionpb <- pb[rowSums(pb) >=10, ]dds <-DESeqDataSetFromMatrix(round(pb), colData = si, design =~ ind + stim)dds <-DESeq(dds)# Contrast: stimulated vs controlres <-results(dds, contrast =c("stim", "stim", "ctrl"))res <- res[order(res$padj), ]head(as.data.frame(res), 8)
sum(res$padj <0.05, na.rm =TRUE) # genes at 5% FDR
[1] 2275
Interpréter la sortie — lisez chaque colonne
DESeq2 renvoie une ligne par gène. Lisez les colonnes, ne vous contentez pas de jeter un œil à la p-value :
baseMean — le comptage normalisé moyen sur les 16 échantillons ; un point d’ancrage grossier du niveau d’expression (un gène à baseMean = 5 est à peine exprimé, un à 9000 est fortement exprimé).
log2FoldChange (le LFC) — la taille d’effet, en unités log2. +1 signifie que le gène est 2× plus élevé dans les cellules stimulées, +3 signifie 8× plus élevé, -1 signifie divisé par deux. C’est ce que vous classez et rapportez, pas la p-value.
lfcSE — l’erreur standard de cette estimation ; la matière première de la statistique de test.
stat, pvalue — la statistique de Wald et sa p-value brute.
padj — la p-value ajustée pour les tests multiples (FDR de Benjamini-Hochberg). C’est celle sur laquelle vous fixez un seuil, jamais la pvalue brute, parce que vous avez testé des milliers de gènes.
Ici, 7 543 gènes étaient testables et 2 275 sont différentiellement exprimés à un FDR de 5 % (1 208 en hausse, 1 067 en baisse) — une réponse large mais crédible, et chacun des meilleurs résultats est un gène que l’interféron est connu pour activer. Comparez cela aux 83 %-de-tout du test par cellule : mêmes cellules, mais une réplication honnête produit un résultat que vous pouvez défendre en relecture. Beaucoup des meilleurs gènes — CCL8, ISG20, PLSCR1, IFITM2, OASL, IRF7, IFIT3, ISG15 — sont les gènes stimulés par l’interféron (ISG) classiques, tous fortement en hausse dans la condition stimulée. Quand le sommet de votre liste est la biologie connue, le pipeline fonctionne.
Une phrase de rapport que vous pouvez coller dans un article :« Au sein des monocytes CD14+, l’expression différentielle pseudobulk (DESeq2, plan ~ donor + condition, n = 8 donneurs par condition) a identifié 2 275 gènes à un FDR de 5 % ; des gènes stimulés par l’interféron dont ISG15, CCL8 et IFIT3 étaient fortement surexprimés après stimulation à l’IFN-β (log2FC > 5). »
Visualiser : le volcano
La figure signature d’un résultat de DE est le volcano plot — la taille d’effet (log2FoldChange) contre la significativité (-log10(padj)) — pour voir d’un coup l’ampleur et la confiance. Les gènes qui dépassent à la fois un seuil de changement d’expression et un seuil de FDR sont mis en évidence, et les ISG les plus forts sont étiquetés.
library(Matrix); library(DESeq2); library(ggplot2); library(ggrepel)d <-readRDS("_data/kang_subset.rds")counts <- d$counts; meta <- d$metamono <- meta$cell =="CD14+ Monocytes"cm <- counts[, mono]; mm <- meta[mono, ]sample_id <-paste(mm$ind, mm$stim, sep ="_")pb <-t(rowsum(t(as.matrix(cm)), group = sample_id))si <-unique(data.frame(sample = sample_id, ind = mm$ind, stim = mm$stim))rownames(si) <- si$sample; si <- si[colnames(pb), ]si$ind <-factor(si$ind)si$stim <-factor(si$stim, levels =c("ctrl", "stim"))pb <- pb[rowSums(pb) >=10, ]dds <-DESeq(DESeqDataSetFromMatrix(round(pb), colData = si, design =~ ind + stim), quiet =TRUE)res <-results(dds, contrast =c("stim", "stim", "ctrl"))df <-as.data.frame(res); df$gene <-rownames(df); df <- df[!is.na(df$padj), ]df$dir <-ifelse(df$padj <0.05&abs(df$log2FoldChange) >1,ifelse(df$log2FoldChange >0, "Up in stim", "Down in stim"), "n.s.")top <-head(df[order(df$padj), ], 12)ggplot(df, aes(log2FoldChange, -log10(padj), color = dir)) +geom_point(size =0.7, alpha =0.6) +scale_color_manual(values =c("Up in stim"="#d1495b", "Down in stim"="#3a86d4", "n.s."="grey80"),name =NULL) +geom_text_repel(data = top, aes(label = gene), color ="black", size =3, max.overlaps =20) +geom_vline(xintercept =c(-1, 1), linetype ="dashed", color ="grey60") +geom_hline(yintercept =-log10(0.05), linetype ="dashed", color ="grey60") +labs(x ="log2 fold change (stim / ctrl)", y ="-log10 adjusted p-value") +theme_minimal(base_size =12)
Le graphique se lit d’un coup d’œil : un haut mur de rouge à droite — la réponse à l’interféron, une surexpression coordonnée de nombreux gènes — avec les ISG qui pointent le plus haut (à la fois de très grands changements d’expression et les plus petites p-values). La biologie n’est pas subtile, ce qui est exactement pourquoi c’est un cas d’enseignement si net : une méthode correcte trouve un signal fort, cohérent et nommé, pas 83 % du transcriptome.
Mises en garde honnêtes : ce dont le pseudobulk a besoin et où sont ses limites
Le pseudobulk est la valeur par défaut correcte, mais ce n’est pas magique. Gardez trois choses à l’esprit :
Il faut une vraie réplication — habituellement au moins 3 donneurs par condition. Le pseudobulk utilise les réplicats biologiques dont vous disposez ; il ne peut pas les inventer. Avec un seul échantillon par bras (le jeu de données d’intégration) il n’y a rien à tester, et DESeq2 refusera ou ne renverra rien. Si votre plan a trop peu de réplicats, c’est un problème de plan expérimental qu’aucune statistique ne peut sauver.
Il faut assez de cellules par échantillon à agréger. Un profil pseudobulk construit à partir d’une poignée de cellules est bruité. Une règle empirique courante est d’exiger un nombre minimum de cellules (disons ~10) par donneur-condition-type cellulaire avant de faire confiance à cet échantillon pseudobulk ; écartez les types cellulaires où un donneur ne contribue quasiment aucune cellule.
DESeq2 est un moteur valide, pas le seul.edgeR et limma-voom sont tout aussi acceptés sur des comptages pseudobulk et s’accordent souvent étroitement. Et le package muscat rationalise l’ensemble du workflow sur tous les types cellulaires à la fois — muscat::aggregateData() construit les matrices pseudobulk et pbDS() exécute la DE par type cellulaire — c’est ce vers quoi vous vous tournerez sur une vraie étude avec une douzaine de types cellulaires plutôt que le seul montré ici.
Quel test quand : marqueurs vs DE inter-conditions
Le test de Wilcoxon par cellule n’est pas faux — il est ici employé pour la mauvaise question. Gardez les deux distincts :
Gènes marqueurs — utilisez le test par cellule (FindAllMarkers / FindMarkers). « Quels gènes distinguent le regroupement 3 du reste ? » est une question sur les cellules, comparant des groupes de cellules que vous avez définis ; traiter chaque cellule comme une observation est approprié, et c’est exactement ce que fait la leçon sur les gènes marqueurs.
DE inter-conditions — utilisez le pseudobulk. « Quels gènes changent entre conditions au sein d’un type cellulaire ? » est une question sur les échantillons/donneurs, donc les donneurs sont les réplicats et le test par cellule pseudoréplique. Cette leçon.
La distinction unique : comparez-vous des cellules (marqueurs → test par cellule) ou des conditions à travers les donneurs (DE → pseudobulk) ? Faites-le correctement et vous ne gonflerez plus jamais un FDR. C’est le même workflow DESeq2 que vous exécuteriez sur du RNA-seq bulk — voir le tutoriel DESeq2 bulk pour les détails de la modélisation des comptages — l’agrégation pseudobulk étant l’unique étape supplémentaire qui l’adapte aux données single-cell.
Problèmes fréquents
DESeq2 dit que presque tous les gènes sont significatifs. Vous lui avez presque certainement fourni des données par cellule au lieu du pseudobulk, ou vous avez agrégé mais oublié de bloquer sur le donneur. Agrégez d’abord (rowsum() / AggregateExpression()) pour que les colonnes soient des échantillons, pas des cellules, et utilisez le plan ~ ind + stim. Si le nombre de colonnes entrant dans DESeqDataSetFromMatrix() se compte en milliers, vous n’avez pas agrégé.
Un modèle ~ ind + stim avec un donneur numérique donne un résultat étrange ou un avertissement. DESeq2 traite l’identifiant de donneur comme une covariable continue. Convertissez-le : si$ind <- factor(si$ind). L’avertissement « the design formula contains one or more numeric variables … did you mean for this to be a factor? » est le signe révélateur.
Peu ou pas de gènes différentiellement exprimés. Habituellement trop peu de réplicats, ou trop peu de cellules par échantillon à agréger. Vérifiez combien de donneurs chaque condition a réellement (table(si$stim)) et combien de cellules sont entrées dans chaque colonne pseudobulk. Le pseudobulk avec 2 contre 2 donneurs a très peu de puissance ; c’est une limite de plan, pas un bug.
Un seul échantillon par condition. Alors vous ne pouvez pas du tout faire d’expression différentielle — il n’y a aucune réplication pour estimer la variabilité. Le pseudobulk n’aide pas ici ; il vous faut plus d’échantillons. (C’est pourquoi le jeu de données d’intégration à 2 échantillons convenait pour aligner les cellules mais était inadapté pour un test de condition.)
Questions fréquentes
NotePourquoi l’expression différentielle par cellule (FindMarkers) est-elle inadaptée à la comparaison de conditions ?
Parce qu’elle traite chaque cellule comme un réplicat indépendant, or les cellules d’un même donneur ne sont pas indépendantes — elles partagent la génétique et le lot de ce donneur. Compter les cellules comme des réplicats (pseudoréplication) gonfle la taille d’échantillon effective de quelques donneurs à des milliers de cellules, ce qui réduit les p-values et gonfle le taux de fausses découvertes. Sur une comparaison nulle sans signal réel, un test de Wilcoxon par cellule déclare tout de même des dizaines de gènes « significatifs » tandis que le pseudobulk n’en déclare correctement aucun.
NoteQu’est-ce que l’expression différentielle pseudobulk ?
Le pseudobulk somme les comptages bruts de toutes les cellules d’un échantillon (un donneur, au sein d’un type cellulaire) en un unique profil d’expression, retransformant une expérience single-cell en une expérience de type bulk avec les donneurs comme réplicats. Vous exécutez ensuite une méthode standard de RNA-seq bulk — DESeq2, edgeR ou limma-voom — sur ces profils par donneur. C’est la méthode recommandée du domaine pour la DE single-cell inter-conditions (Squair 2021 ; Crowell 2020).
NoteComment fabriquer une matrice pseudobulk en R à partir d’un objet Seurat ?
Utilisez Seurat::AggregateExpression(obj, group.by = c("donor", "condition")), qui somme les comptages de chaque groupe en une colonne. Sans Seurat, c’est une ligne de R de base sur la matrice de comptage : pb <- t(rowsum(t(as.matrix(counts)), group = paste(donor, condition, sep = "_"))). Les deux donnent une matrice gènes × échantillons que vous passez à DESeqDataSetFromMatrix().
NoteCombien de réplicats me faut-il pour une DE pseudobulk ?
Les donneurs sont vos réplicats, il vous en faut donc plusieurs par condition — un minimum courant est d’au moins 3 par groupe, et plus c’est mieux. Le pseudobulk ne peut pas créer une réplication que vous n’avez pas collectée : avec un seul échantillon par condition il n’y a aucune variabilité à estimer et aucun test valide. Vous voulez aussi assez de cellules par donneur-condition-type cellulaire (environ ≥10) pour que chaque profil pseudobulk ne soit pas dominé par le bruit.
NoteDois-je utiliser DESeq2, edgeR ou limma-voom pour le pseudobulk ?
Les trois sont valides sur des comptages pseudobulk et s’accordent généralement étroitement ; DESeq2 et edgeR modélisent les comptages directement, limma-voom les transforme d’abord. Choisissez celui que vous connaissez. Pour une vraie étude avec de nombreux types cellulaires, le package muscat enveloppe tout le workflow (aggregateData() + pbDS()) pour que vous exécutiez la DE pseudobulk sur chaque type cellulaire d’un coup plutôt qu’un à la fois.
Testez vos connaissances
ImportantÀ vous de jouer : exécuter une DE pseudobulk sur les lymphocytes B
Refaites le workflow, mais sur les lymphocytes B au lieu des monocytes. Sous-ensemblez avec meta$cell == "B cells", agrégez en un profil pseudobulk par donneur et par condition, exécutez DESeq2 avec le plan ~ ind + stim (rappelez-vous factor(ind)), et rapportez combien de gènes sont significatifs à un FDR de 5 %. Les meilleurs résultats sont-ils toujours des gènes stimulés par l’interféron ?
AstuceIndice
Le seul changement par rapport au code des monocytes est la ligne de sous-ensemblage : bcell <- meta$cell == "B cells". Tout ce qui vient ensuite — paste(ind, stim), rowsum(), factor(ind), le plan ~ ind + stim, le contraste stim vs ctrl — est identique. Comptez avec sum(res$padj < 0.05, na.rm = TRUE) et lisez les premières lignes avec head(res[order(res$padj), ]).
AstuceSolution
library(Matrix); library(DESeq2)d <-readRDS("_data/kang_subset.rds")counts <- d$counts; meta <- d$metabcell <- meta$cell =="B cells"cm <- counts[, bcell]; mm <- meta[bcell, ]sample_id <-paste(mm$ind, mm$stim, sep ="_")pb <-t(rowsum(t(as.matrix(cm)), group = sample_id))si <-unique(data.frame(sample = sample_id, ind = mm$ind, stim = mm$stim))rownames(si) <- si$sample; si <- si[colnames(pb), ]si$ind <-factor(si$ind)si$stim <-factor(si$stim, levels =c("ctrl", "stim"))pb <- pb[rowSums(pb) >=10, ]dds <-DESeq(DESeqDataSetFromMatrix(round(pb), si, design =~ ind + stim), quiet =TRUE)res <-results(dds, contrast =c("stim", "stim", "ctrl"))sum(res$padj <0.05, na.rm =TRUE) # significant genes at 5% FDRhead(res[order(res$padj), ]) # top hits — still ISGs (ISG20, ISG15, LY6E, ...)
Les lymphocytes B montent le même programme interféron : le workflow est inchangé, et le sommet de la liste est de nouveau dominé par des gènes stimulés par l’interféron. Cette reproductibilité à travers les types cellulaires est la biologie — chaque type cellulaire répond à l’IFN-β, donc une méthode correcte retrouve la même signature ISG dans chacun.
NoteVérification rapide : votre matrice pseudobulk a 2 400 colonnes et DESeq2 signale 80 % des gènes comme significatifs. Qu’est-ce qui a mal tourné ?
A. Rien — une vraie réponse à l’interféron peut toucher 80 % des gènes. B. Vous n’avez jamais agrégé : 2 400 colonnes signifie que la matrice est encore à une colonne par cellule, donc DESeq2 pseudoréplique tout comme FindMarkers. C. Vous avez utilisé la mauvaise direction de contraste.
AstuceAfficher la réponse
B. Une matrice pseudobulk devrait avoir une colonne par échantillon (donneur × condition) — ici c’est 16, pas 2 400. Des milliers de colonnes signifient que vous avez passé la matrice par cellule directement à DESeq2 sans sommer, donc il traite les cellules comme des réplicats et la significativité gonflée est l’artefact de pseudoréplication que cette leçon existe pour prévenir. Agrégez d’abord avec rowsum() / AggregateExpression().
Conclusion
L’expression différentielle inter-conditions en données single-cell a une bonne et une mauvaise façon de faire, et la différence tient à quelle unité vous traitez comme un réplicat. La mauvaise façon — un test par cellule comme FindMarkers — compte des milliers de cellules comme des milliers d’échantillons, pseudoréplique, et gonfle le FDR au point d’inventer des dizaines de résultats à partir de pur bruit. La bonne façon, c’est le pseudobulk : sommez les cellules de chaque donneur en un profil, puis exécutez DESeq2 (ou edgeR / limma-voom) avec le plan ~ donor + condition et le donneur comme facteur, de sorte que l’effet de condition soit estimé à travers de véritables réplicats biologiques. Sur les monocytes IFN-β de Kang, cela donne un résultat fort et interprétable — 2 275 gènes, menés par les gènes stimulés par l’interféron canoniques — qui survivra à la relecture. Rappelez-vous la distinction unique qui décide de la méthode : comparer des cellules est une question de marqueurs (test par cellule) ; comparer des conditions à travers les donneurs est une question de DE (pseudobulk).
Demandez à Prova« comment exécuter une expression différentielle pseudobulk sur mon propre objet Seurat multi-échantillons avec DESeq2 ? » — elle répond avec du code R que vous pouvez exécuter sur vos propres données : agréger les cellules de chaque donneur avec AggregateExpression(), construire le colData, ajuster DESeq2 avec le donneur comme facteur de blocage, et lire le résultat — au lieu de pseudorépliquer avec FindMarkers. The runtime is the judge.Ask Prova →
Cette page vous a-t-elle été utile ?
Merci pour votre retour !
Recevez les nouvelles leçons R & Python par e-mail
Pratique, reproductible, sans spam. Désinscription à tout moment.
Double opt-in. Nous ne partageons jamais votre e-mail.
Cette leçon est reproductible : chaque figure et chaque nombre ont été produits par le code montré — copiez n’importe quel bloc et exécutez-le pour les reproduire. The runtime is the judge.
Crowell, H. L., et al. (2020). muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nature Communications, 11, 6077. https://www.nature.com/articles/s41467-020-19894-4
Murphy, A. E., & Skene, N. G. (2022). A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis. Nature Communications, 13, 7851. https://www.nature.com/articles/s41467-022-35519-4
Kang, H. M., et al. (2018). Multiplexed droplet single-cell RNA-sequencing using natural genetic variation — the source of the Kang IFN-β PBMC dataset. Nature Biotechnology, 36(1), 89–94. https://www.nature.com/articles/nbt.4042
