ACP et regroupement d’échantillons pour le RNA-seq en R (DESeq2)

Contrôlez vos échantillons avant l’expression différentielle — stabilisez la variance, réalisez une ACP et regroupez les échantillons pour repérer le groupement, les valeurs aberrantes et les effets de lot

Un guide pratique du contrôle qualité exploratoire pour le RNA-seq bulk en R avec DESeq2. Avant de tester l’expression différentielle, vérifiez que vos échantillons sont cohérents : appliquez une transformation stabilisatrice de variance (VST), tracez un graphique d’ACP pour voir si les échantillons se séparent par condition (et pour repérer les valeurs aberrantes ou les effets de lot), et une heatmap des distances entre échantillons pour voir comment ils se regroupent. Lisez les axes en pourcentage de variance en langage clair. Sur le jeu de données airway.

AstucePoints clés

• Examinez vos échantillons avant de les tester. Une vérification rapide par ACP + regroupement repère les échantillons mal étiquetés, les valeurs aberrantes et les effets de lot avant qu’ils ne ruinent votre expression différentielle.
• Stabilisez d’abord la variance. Les comptages bruts ont une variance qui croît avec la moyenne, si bien que quelques gènes bruyants domineraient. La VST de DESeq2 (vst()) place tous les gènes sur une échelle comparable pour les étapes visuelles.
• L’ACP projette les échantillons sur les directions de plus grande variation. Si vos groupes biologiques se séparent le long de PC1 ou PC2, le signal est fort ; si c’est un lot qui les sépare, vous avez un problème à modéliser.
• Une heatmap des distances entre échantillons est la vue complémentaire : elle organise les échantillons en regroupements pour voir lesquels vont ensemble et repérer ceux qui restent isolés.
• Ce sont des étapes de CQ, pas le test — elles s’appliquent aux comptages transformés ; le test d’expression différentielle s’applique aux comptages bruts (leçon DESeq2).

Introduction

Avant de demander quels gènes ont changé, posez une question plus simple : mes échantillons sont-ils cohérents ? Les expériences de séquençage sont pleines de désastres évitables — une étiquette inversée, une librairie ratée, un lot traité un jour différent. Une vérification exploratoire de deux minutes les détecte avant qu’ils ne corrompent tous les résultats en aval.

Cette vérification comporte deux vues standard, toutes deux construites sur une version à variance stabilisée des comptages : un graphique d’ACP (les échantillons se séparent-ils par condition ?) et une heatmap des distances entre échantillons (comment se regroupent-ils ?). Ce guide exécute les deux dans R avec DESeq2 sur le jeu de données airway — des cellules musculaires lisses des voies respiratoires provenant de quatre donneurs, chacune traitée et non traitée à la dexaméthasone.

Partez des données de comptage

Nous construisons l’objet DESeq2 à partir des comptages bruts — le même point de départ que la leçon DESeq2.

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

class: DESeqDataSet 
dim: 22369 8 
metadata(2): '' version
assays(1): counts
rownames(22369): ENSG00000000003 ENSG00000000419 ... ENSG00000273487
  ENSG00000273488
rowData names(10): gene_id gene_name ... seq_coord_system symbol
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample

Stabilisez la variance (VST)

Les comptages bruts sont une mauvaise entrée pour l’ACP et le regroupement : leur variance croît avec la moyenne, si bien qu’une poignée de gènes fortement exprimés domineraient uniquement parce qu’ils sont grands. La transformation stabilisatrice de variance (vst) corrige cela — elle place les gènes à faible et fort comptage sur une échelle comparable, comme un log2 plus malin.

vsd <- vst(dds, blind = TRUE)     # blind = TRUE: ignore the design, a pure QC view
head(assay(vsd))[, 1:4]

                SRR1039508 SRR1039509 SRR1039512 SRR1039513
ENSG00000000003   9.464955   9.085088   9.613950   9.349298
ENSG00000000419   8.963799   9.271415   9.154078   9.199673
ENSG00000000457   8.211542   8.117807   8.051019   8.165650
ENSG00000000460   6.610550   6.652804   6.153025   6.507349
ENSG00000000971  11.652764  12.018678  12.366560  12.641598
ENSG00000001036  10.495488  10.263000  10.563194  10.407888

Nous fixons blind = TRUE pour que la transformation n’utilise pas le plan expérimental — exactement ce qu’il vous faut pour un regard non biaisé de contrôle qualité. (Pour la modélisation en aval, vous utiliseriez blind = FALSE.) Pour de très petites études, rlog() peut être plus robuste ; pour des expériences typiques et plus grandes, VST est plus rapide et le choix standard.

ACP — les échantillons se séparent-ils par condition ?

L’analyse en composantes principales projette chaque échantillon sur les directions qui capturent le plus de variation. PC1 est la plus grande source de variance, PC2 la suivante. Si votre groupement biologique s’aligne sur PC1 ou PC2, l’effet est fort par rapport au bruit.

library(ggplot2)

pcaData <- plotPCA(vsd, intgroup = c("dex", "cell"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(PC1, PC2, color = dex, shape = cell)) +
  geom_point(size = 4) +
  scale_color_viridis_d(end = 0.85) +
  labs(x = paste0("PC1: ", percentVar[1], "% variance"),
       y = paste0("PC2: ", percentVar[2], "% variance"),
       color = "Treatment", shape = "Donor") +
  coord_fixed() +
  theme_minimal()

Lisez le graphique. PC1 capture 41 % de la variance et sépare nettement les échantillons traités des non traités (gauche vs droite) — l’effet de la dexaméthasone est la plus grande source de variation, un signe rassurant avant de tester. PC2 (26 %) capture la variation d’un donneur à l’autre (le plus visiblement un donneur en haut, un autre en bas), les différences interindividuelles que nous contrôlons avec cell dans le plan. Aucun échantillon ne se tient à l’écart, il n’y a donc pas de valeurs aberrantes évidentes.

AvertissementSi un lot sépare vos échantillons sur PC1

Le signe d’alerte, c’est quand PC1 s’aligne sur une variable technique — jour de séquençage, voie, kit — plutôt que sur votre biologie. Cela signifie qu’un effet de lot est le signal dominant. La solution est d’ajouter le lot à votre plan (~ batch + condition) pour que DESeq2 en tienne compte, pas de l’ignorer. L’ACP est ce qui permet de le repérer.

Heatmap des distances entre échantillons

L’ACP résume les échantillons en deux dimensions ; une heatmap des distances montre l’image complète par paires. Calculez les distances entre les échantillons sur les données VST et regroupez-les — les échantillons d’un même groupe devraient se trouver ensemble.

library(pheatmap)

sampleDists <- dist(t(assay(vsd)))
mat <- as.matrix(sampleDists)
rownames(mat) <- paste(vsd$dex, vsd$cell, sep = " · ")
colnames(mat) <- NULL

pheatmap(mat,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         main = "Sample-to-sample distances (VST)")

Le dendrogramme raconte une histoire plus subtile que l’ACP. Un donneur (N080611) domine : ses échantillons traités et non traités se regroupent et se détachent du reste en premier ; ce n’est que parmi les six restants que le traitement pilote ensuite la séparation. Pourquoi cette différence avec l’ACP ? L’ACP a utilisé les 500 gènes les plus variables — là où l’effet du traitement ressort — tandis que ces distances utilisent tous les gènes, de sorte qu’un seul donneur fortement différent peut dominer l’image globale. Que les deux vues mettent en avant une structure différente est en soi instructif : cela confirme que l’effet du donneur est réel et important, ce qui explique précisément pourquoi le plan le contrôle avec ~ cell + dex. Lisez toujours les deux — ils répondent à des questions légèrement différentes.

Problèmes fréquents

Les échantillons ne se séparent pas par condition. Soit l’effet est réellement faible (la GSEA à l’étape suivante peut tout de même trouver des décalages coordonnés), soit une variable de nuisance plus forte domine. Colorez l’ACP selon d’autres variables (lot, RIN, donneur) pour découvrir ce que PC1 suit réellement.

Un échantillon se tient loin des autres. Une valeur aberrante probable — une librairie ratée ou une inversion. Vérifiez sa taille de librairie et ses indicateurs de CQ ; envisagez de l’écarter, mais jamais en silence. Documentez la décision.

Vous avez réalisé l’ACP sur des comptages bruts ou en log simple. À éviter. Sans VST (ou rlog), les gènes à fort comptage dominent la variance et le graphique reflète la profondeur de séquençage plus que la biologie. Transformez toujours d’abord.

Questions fréquentes

NoteVST ou rlog — quelle transformation utiliser ?

Les deux stabilisent la variance pour la visualisation. vst() est rapide et le choix par défaut pour les études typiques ; rlog() est plus robuste sur les très petits jeux de données (une poignée d’échantillons) mais bien plus lent sur les grands. Pour les graphiques de CQ, l’un ou l’autre convient — utilisez VST sauf si vous avez très peu d’échantillons.

Noteblind doit-il être TRUE ou FALSE dans vst() ?

Pour le contrôle qualité, utilisez blind = TRUE : la transformation ignore le plan, offrant un regard non biaisé sur le fait que les échantillons se regroupent comme attendu. Pour des valeurs transformées que vous fournirez à la modélisation en aval, utilisez blind = FALSE afin de préserver la variance liée au plan.

NotePourquoi ne pas simplement réaliser l’ACP sur les comptages bruts ?

Parce que la variance des comptages croît avec la moyenne : quelques gènes fortement exprimés domineraient les composantes principales, et le graphique refléterait la profondeur de séquençage plus que la biologie. La VST place d’abord tous les gènes sur une échelle comparable, de sorte que l’ACP voit la vraie structure.

NoteEt si PC1 sépare par lot au lieu de la condition ?

C’est un effet de lot qui domine vos données. Ne l’ignorez pas — ajoutez la variable de lot au plan DESeq2 (~ batch + condition) pour que le test en tienne compte. L’ACP est précisément l’outil qui révèle le problème avant de tester.

Testez vos connaissances

ImportantÀ vous de jouer : quelle variable pilote PC2 ?

L’ACP ci-dessus colore par traitement. Retracez-la en colorant plutôt par donneur (cell), et déterminez sur quelle composante principale s’aligne l’effet du donneur. Est-ce PC1 ou PC2 ?

AstuceIndice

Réutilisez pcaData ; changez simplement l’attribut esthétique color = en cell. Regardez quel axe écarte les donneurs.

AstuceSolution

ggplot(pcaData, aes(PC1, PC2, color = cell, shape = dex)) +
  geom_point(size = 4) +
  coord_fixed() + theme_minimal()

Les donneurs s’étalent le long de PC2, pas de PC1 — le traitement est l’axe dominant (PC1, 41 %), le donneur le secondaire (PC2, 26 %). C’est exactement pourquoi le plan est ~ cell + dex : il retire la variance du donneur pour que le test puisse voir proprement l’effet du traitement.

NoteVérification rapide : pourquoi transformer les comptages avant l’ACP ?

A. Pour rendre le test plus significatif. B. Parce que la variance des comptages bruts croît avec la moyenne, si bien que les grands gènes domineraient les composantes. C. L’ACP ne sait pas gérer les entiers.

AstuceAfficher la réponse

B. La variance des comptages croît avec le niveau d’expression, donc sans transformation stabilisatrice de variance, quelques gènes fortement exprimés pilotent les composantes principales et le graphique suit la profondeur de séquençage, pas la biologie.

Conclusion

Vous avez réalisé la vérification exploratoire standard pour le RNA-seq bulk : stabilisez la variance des comptages (vst), puis examinez les échantillons de deux façons — un graphique d’ACP (ici PC1 a capturé 41 % de la variance et a nettement séparé les traités des non traités, avec la variation des donneurs sur PC2) et une carte des distances entre échantillons (qui, en utilisant tous les gènes, était dominée par un seul donneur — un rappel que les deux vues peuvent mettre en avant une structure différente). La règle pratique : faites toujours cela avant l’expression différentielle — c’est le moyen le moins coûteux de détecter une étiquette inversée, une valeur aberrante ou un effet de lot, et cela vous indique si votre effet est assez fort pour qu’on s’y fie. Une fois les échantillons contrôlés, vous êtes prêt à tester l’expression différentielle.

Leçons connexes

• Expression différentielle (DESeq2) — le test que cette étape de CQ précède. · Enrichissement de gènes & analyse de voies — quoi faire avec les gènes que vous trouvez. · Analyse RNA-seq bulk — la série complète. · Bioinformatique — le pilier.
• Technique connexe : Analyse en composantes principales — la méthode générale, dans le pilier Machine Learning.

Cette page vous a-t-elle été utile ?

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.

Partager cette pageXLinkedInRedditHN

Note

Cette leçon est reproductible : chaque figure a été produite par le code montré — copiez n’importe quel bloc et exécutez-le pour les reproduire. The runtime is the judge.

Références

• Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550.
• Love, M. I., Anders, S., & Huber, W. Analyzing RNA-seq data with DESeq2 — la vignette Bioconductor (transformations des données & évaluation de la qualité).
• Himes, B. E., et al. (2014). RNA-Seq transcriptome profiling identifies CRISPLD2 as a glucocorticoid responsive gene in airway smooth muscle cells (le jeu de données airway). PLoS ONE, 9(6), e99625.