Analyse de pseudotemps et de trajectoire en R avec Slingshot

Ordonnez les cellules le long d’une trajectoire de développement avec Slingshot, enracinez-la dans la vraie population de cellules souches, et trouvez les gènes qui s’activent le long de chaque lignée avec tradeSeq

Un guide pratique de l’inférence de trajectoire et du pseudotemps pour le RNA-seq single-cell en R. Regroupez et calculez l’embedding d’une population de cellules souches en cours de différenciation, ajustez des courbes de lignées ramifiées et un pseudotemps avec Slingshot (un arbre couvrant de poids minimal sur les centroïdes des regroupements lissé en courbes principales), validez la racine et les extrémités face aux étiquettes FACS connues, puis testez quels gènes changent le long de chaque lignée avec tradeSeq — sur le jeu de données canonique de cellules souches/progénitrices du sang de souris Nestorowa.

Date de publication

3 juillet 2026

Modifié

7 juillet 2026

AstucePoints clés
  • Une partie de la biologie est un continuum, pas un ensemble de cases. Quand les cellules se différencient — des cellules souches qui mûrissent vers plusieurs destins matures —, la structure intéressante est l’ordonnancement le long de ce continuum, pas une poignée d’étiquettes de regroupements discrètes. Le pseudotemps est cet ordonnancement : un nombre par cellule mesurant jusqu’où elle a progressé le long d’une trajectoire de développement.
  • Slingshot ajuste la trajectoire en deux étapes. Il construit un arbre couvrant de poids minimal (MST) sur les centroïdes des regroupements pour établir la topologie de ramification, puis ajuste des courbes principales lisses à travers les cellules pour transformer ce squelette en lignées continues et en un pseudotemps (Street et al., 2018). Un seul appel : slingshot(embedding, clusterLabels = clusters).
  • La racine fixe la direction — choisissez-la, n’acceptez pas la valeur par défaut à l’aveugle. Slingshot ordonne les cellules le long de l’arbre mais ne sait pas quelle extrémité est le départ. Pointez-le vers l’origine connue avec start.clus ; ici les étiquettes FACS nous disent quel regroupement est la cellule souche à long terme, nous y enracinons donc la trajectoire.
  • tradeSeq trouve les gènes qui changent le long d’une lignée. Une fois que Slingshot vous a donné l’ordonnancement, tradeSeq ajuste une courbe lisse (un modèle additif généralisé, GAM) de l’expression de chaque gène contre le pseudotemps et teste quels gènes varient — les gènes qui pilotent chaque destin (Van den Berge et al., 2020).
  • Le pseudotemps est une direction, pas une horloge, et Slingshot ne sait pas dans quel sens le temps s’écoule. Il récupère l’ordonnancement et la topologie ; il ne peut pas distinguer seul un début d’une fin. Obtenir la direction à partir des données (plutôt que d’une racine connue) requiert la RNA velocity, une méthode d’abord Python couverte dans une leçon compagne.

Introduction

Vous avez regroupé un jeu de données de RNA-seq single-cell (scRNA-seq) et nommé les types cellulaires. Cela fonctionne à merveille quand le tissu est un ensemble de types cellulaires stables et distincts — des monocytes ici, des lymphocytes B là. Mais une partie de la biologie la plus importante n’est pas discrète du tout. La formation du sang en est le cas d’école : une cellule souche à long terme ne saute pas à « neutrophile mature » ; elle s’écoule à travers un continuum d’états progéniteurs, s’engageant progressivement vers l’un de plusieurs destins — globules rouges, plaquettes, granulocytes, lymphocytes. Tracez des frontières de regroupements rigides là-dessus et vous avez jeté la chose même qui vous importe : l’ordonnancement.

Cette leçon répond à une autre question que ne le fait le clustering. Non pas « quel est le type de chaque cellule ? » mais « jusqu’où en est chaque cellule dans son parcours de développement, et quels gènes l’y conduisent ? » L’outil pour l’ordonnancement est Slingshot, et l’outil pour les gènes est tradeSeq. Vous regrouperez et calculerez l’embedding d’une population de cellules souches en cours de différenciation, ajusterez des courbes de lignées ramifiées et un pseudotemps avec Slingshot, vérifierez que la racine et les extrémités inférées correspondent à la biologie connue, puis testerez quels gènes changent le long de chaque lignée.

Nous travaillons sur le jeu de données de sang de souris Nestorowa canonique (Nestorowa et al., 2016) : 1 920 cellules souches et progénitrices hématopoïétiques (HSPC — les cellules souches et progénitrices qui produisent chaque cellule sanguine). Fait crucial, chaque cellule porte une étiquette FACS (fluorescence-activated cell sorting — le tri par cytométrie de flux qui a trié les cellules selon leurs marqueurs de surface) indiquant à quelle population progénitrice elle appartient. Nous ne fournissons jamais ces étiquettes à Slingshot ; nous les mettons de côté pour vérifier que la trajectoire qu’il infère est biologiquement réelle.

Ce qu’est le pseudotemps — et quand l’utiliser

Le pseudotemps est un unique nombre par cellule décrivant sa position le long d’un continuum de différenciation : les cellules au début (les moins différenciées) reçoivent une valeur basse, les cellules à la fin (pleinement engagées) reçoivent une valeur élevée. C’est un pseudo-temps parce qu’il n’est pas mesuré en heures — il est inféré de la façon dont l’expression change graduellement d’une cellule à l’autre. Une trajectoire est le chemin à travers ce continuum, et elle peut se ramifier : une population souche peut s’écouler vers plusieurs destins distincts, si bien qu’une trajectoire est souvent un arbre de lignées, une par destin.

Utilisez l’analyse de trajectoire quand vos cellules forment un véritable continuum — différenciation, activation, une réponse à un traitement qui se déploie le long d’un gradient d’états. N’imposez pas l’analyse à des types cellulaires discrets et stables : demander un ordonnancement en pseudotemps de monocytes, de lymphocytes T et de lymphocytes B n’a aucun sens, car il n’existe aucun chemin de développement entre eux le long duquel ordonner. C’est le contraste net avec l’annotation des types cellulaires, qui répond à la question opposée — étiqueter chaque cellule comme l’un d’un ensemble de types stables. L’annotation est faite pour les cases ; l’analyse de trajectoire est faite pour les gradients. Le premier jugement dans toute analyse de trajectoire est de décider que vous avez réellement un continuum.

Les données : des cellules souches avec des étiquettes progénitrices connues

Le nestorowa_hsc.rds commité est une tranche de l’atlas HSPC de Nestorowa : une matrice de comptage brute de 2 500 gènes hautement variables (HVG) à travers 1 920 cellules, plus une table de métadonnées. La colonne broad contient la population FACS dans laquelle chaque cellule a été triée — LTHSC (cellule souche hématopoïétique à long terme, LT-HSC — la cellule souche au sommet), MPP (progéniteur multipotent), LMPP (progéniteur multipotent à orientation lymphoïde), CMP (progéniteur myéloïde commun), GMP (progéniteur granulocyte-macrophage — le destin myéloïde), et MEP (progéniteur mégacaryocyte-érythroïde — le destin globule rouge / plaquette). Ce gradient de la cellule souche au progéniteur engagé est la vérité de terrain de la différenciation face à laquelle nous validons. Les noms de gènes sont stockés comme des identifiants Ensembl avec une table symbol compagne pour que nous puissions étiqueter les figures avec des noms de gènes lisibles.

library(Seurat)

d <- readRDS("_data/nestorowa_hsc.rds")
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
sym <- d$symbol                  # named vector: Ensembl ID -> gene symbol, for labelling

dim(obj)                          # genes x cells
[1] 2500 1920
table(obj$broad, useNA = "ifany") # the FACS populations: our differentiation ground truth

  CMP   GMP  LMPP LTHSC   MEP   MPP  <NA> 
  328   123   280   323   362   368   136

Six populations triées, du sommet LT-HSC jusqu’aux progéniteurs GMP et MEP engagés — plus un groupe de cellules sans étiquette broad (les HSPC indifférenciées qui n’ont pas été assignées à une population fine). Chaque bloc ci-dessous repart exactement de ce point, de sorte que vous pouvez copier n’importe lequel d’entre eux et l’exécuter seul.

Regrouper les cellules et calculer leur embedding

Slingshot ne travaille pas sur les cellules brutes ; il travaille sur des regroupements dans un embedding de faible dimension. Le premier travail est donc le pipeline scRNA-seq standard : normaliser, choisir les gènes variables, réduire les dimensions par analyse en composantes principales (PCA), construire un UMAP (uniform manifold approximation and projection — la carte en 2-D), et regrouper. C’est exactement le workflow clustering & UMAP ; Slingshot reprend là où il s’arrête.

Nous fixons d’abord une graine aléatoire : l’UMAP et le clustering utilisent tous deux du hasard, et une trajectoire bâtie sur des numéros de regroupements changeants n’est pas reproductible. Fixez la graine et toute l’analyse est stable.

set.seed(42)
library(Seurat)
library(ggplot2)
d <- readRDS("_data/nestorowa_hsc.rds")
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)

obj <- NormalizeData(obj, verbose = FALSE)
obj <- FindVariableFeatures(obj, nfeatures = 2000, verbose = FALSE)
obj <- ScaleData(obj, verbose = FALSE)
obj <- RunPCA(obj, npcs = 30, verbose = FALSE)
obj <- RunUMAP(obj, dims = 1:30, verbose = FALSE)
obj <- FindNeighbors(obj, dims = 1:30, verbose = FALSE)
obj <- FindClusters(obj, resolution = 0.6, verbose = FALSE)

nlevels(obj$seurat_clusters)                 # how many clusters
[1] 8
DimPlot(obj, group.by = "seurat_clusters", label = TRUE) + ggtitle("Clusters on the UMAP")

A UMAP scatter plot of 1,920 blood stem and progenitor cells, coloured by the eight clusters found by the standard Seurat pipeline. The clusters form a connected, branching shape rather than well-separated islands — the continuous manifold expected of a differentiating population, where one state flows into the next instead of sitting apart.

Huit regroupements — mais remarquez la forme. Contrairement aux îlots bien séparés d’un jeu de données de tissu mixte, ces regroupements se rejoignent en une seule variété connectée et ramifiée. C’est exactement à quoi ressemble une population en cours de différenciation : chaque état se fond dans le suivant, et les ramifications sont les débuts des destins engagés. Un UMAP continu et ramifié comme celui-ci est le signal visuel que l’analyse de trajectoire est le bon outil.

Ajuster la trajectoire avec Slingshot

Slingshot transforme ces regroupements en lignées en deux étapes conceptuelles (Street et al., 2018) :

  1. Topologie — un MST sur les centroïdes des regroupements. Il calcule le centre de chaque regroupement dans l’embedding et relie les centres par un arbre couvrant de poids minimal : le plus petit ensemble d’arêtes qui relie chaque regroupement sans former de boucle. Cet arbre est le squelette de ramification — quel état s’écoule dans quel autre.
  2. Ordonnancement — des courbes principales lisses. Un arbre de segments rectilignes reliant les centroïdes est en dents de scie, donc Slingshot ajuste des courbes principales simultanées à travers les cellules réelles le long de chaque chemin de la racine à l’extrémité. Chaque courbe est une lignée lisse, et la position d’une cellule le long de sa courbe est son pseudotemps.

L’appel prend l’embedding et les étiquettes de regroupements. Il prend aussi une racine : Slingshot ordonne les cellules le long de l’arbre mais ne peut pas dire quelle extrémité est le début, donc laissé à lui-même il infère une racine — atterrissant souvent sur un regroupement indifférencié, mais sans garantie que ce soit celui que vous visez. Ici nous en savons plus : les étiquettes FACS nous disent quel regroupement est la population souche LT-HSC, donc nous y enracinons la trajectoire avec start.clus.

set.seed(42)
library(Seurat)
library(slingshot)
d <- readRDS("_data/nestorowa_hsc.rds")
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
obj <- NormalizeData(obj, verbose = FALSE) |>
  FindVariableFeatures(nfeatures = 2000, verbose = FALSE) |> ScaleData(verbose = FALSE) |>
  RunPCA(npcs = 30, verbose = FALSE) |> RunUMAP(dims = 1:30, verbose = FALSE) |>
  FindNeighbors(dims = 1:30, verbose = FALSE) |> FindClusters(resolution = 0.6, verbose = FALSE)

# Which cluster is the stem-cell root? The one most enriched for LT-HSC in the FACS labels.
lthsc_frac <- prop.table(table(obj$seurat_clusters, obj$broad), 1)[, "LTHSC"]
root <- names(which.max(lthsc_frac))
round(sort(lthsc_frac, decreasing = TRUE), 3)   # LT-HSC fraction per cluster
    5     0     6     7     3     2     1     4 
0.578 0.343 0.192 0.075 0.070 0.070 0.052 0.035 
# Fit lineages + pseudotime on the UMAP, rooted at the LT-HSC cluster
sds <- slingshot(Embeddings(obj, "umap"),
                 clusterLabels = as.character(obj$seurat_clusters),
                 start.clus = root)

slingLineages(sds)               # the branching lineages, each a path of clusters from the root
$Lineage1
[1] "5" "0" "3" "1" "4"

$Lineage2
[1] "5" "0" "3" "1" "2"

$Lineage3
[1] "5" "0" "3" "1" "7"

$Lineage4
[1] "5" "0" "6"

Le regroupement le plus enrichi en LT-HSC (regroupement 5, ~58 % de LT-HSC selon la FACS) est la racine, et Slingshot résout quatre lignées, partant toutes de lui et s’évasant vers différentes extrémités — une population souche se ramifiant en plusieurs destins, exactement l’arbre que l’hématopoïèse devrait produire. Chaque lignée est un chemin de regroupements : elles partagent un tronc précoce (5 → 0 → 3 → 1) puis se séparent vers des extrémités distinctes, plus une courte ramification (5 → 0 → 6) vers un regroupement enrichi en LMPP, à orientation lymphoïde.

Voir la trajectoire : les courbes de lignées sur l’UMAP

La figure signature est celle des courbes de lignées ajustées, tracées par-dessus les cellules. Extrayez les courbes lisses de l’objet Slingshot avec slingCurves() et posez-les sur l’UMAP : vous voyez le squelette de ramification que l’algorithme a trouvé, s’incurvant hors de la racine et se séparant vers les destins.

set.seed(42)
library(Seurat)
library(slingshot)
library(ggplot2)
d <- readRDS("_data/nestorowa_hsc.rds")
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
obj <- NormalizeData(obj, verbose = FALSE) |>
  FindVariableFeatures(nfeatures = 2000, verbose = FALSE) |> ScaleData(verbose = FALSE) |>
  RunPCA(npcs = 30, verbose = FALSE) |> RunUMAP(dims = 1:30, verbose = FALSE) |>
  FindNeighbors(dims = 1:30, verbose = FALSE) |> FindClusters(resolution = 0.6, verbose = FALSE)
root <- names(which.max(prop.table(table(obj$seurat_clusters, obj$broad), 1)[, "LTHSC"]))
sds  <- slingshot(Embeddings(obj, "umap"), clusterLabels = as.character(obj$seurat_clusters),
                  start.clus = root)

um <- Embeddings(obj, "umap")
cells <- data.frame(UMAP1 = um[, 1], UMAP2 = um[, 2], cluster = obj$seurat_clusters)
# Each fitted curve is a set of ordered points; pull them out for plotting
curves <- do.call(rbind, lapply(seq_along(slingCurves(sds)), function(i) {
  s <- slingCurves(sds)[[i]]$s[slingCurves(sds)[[i]]$ord, ]
  data.frame(UMAP1 = s[, 1], UMAP2 = s[, 2], lineage = paste0("Lineage ", i))
}))

ggplot(cells, aes(UMAP1, UMAP2)) +
  geom_point(aes(color = cluster), size = 0.6, alpha = 0.6) +
  geom_path(data = curves, aes(group = lineage), linewidth = 1.1) +
  labs(title = "Slingshot lineage curves", color = "Cluster") +
  theme_minimal(base_size = 12)

A UMAP of 1,920 blood stem and progenitor cells coloured by cluster, with four smooth black Slingshot lineage curves drawn over the points. All four curves start together at the stem-cell root cluster on the right and fan out to the left, each ending in a different committed-progenitor region. The curves trace the continuous branching of one stem population into several blood fates.

Les quatre courbes émergent d’un point unique — la racine LT-HSC — et se ramifient vers l’extérieur, chacune se faufilant à travers les cellules vers un coin différent de la carte. C’est toute l’idée de l’inférence de trajectoire en une image : un chemin continu et ramifié ajusté à travers une population en cours de différenciation, remplaçant les frontières de regroupements rigides par la structure de développement qui a réellement engendré les cellules.

Colorer la carte par le pseudotemps

La distance de chaque cellule le long de sa courbe est son pseudotemps. Moyenner le pseudotemps d’une cellule à travers les lignées auxquelles elle appartient et colorer l’UMAP par cette valeur montre directement le gradient de différenciation : sombre à la racine, clair aux extrémités engagées.

set.seed(42)
library(Seurat)
library(slingshot)
library(ggplot2)
d <- readRDS("_data/nestorowa_hsc.rds")
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
obj <- NormalizeData(obj, verbose = FALSE) |>
  FindVariableFeatures(nfeatures = 2000, verbose = FALSE) |> ScaleData(verbose = FALSE) |>
  RunPCA(npcs = 30, verbose = FALSE) |> RunUMAP(dims = 1:30, verbose = FALSE) |>
  FindNeighbors(dims = 1:30, verbose = FALSE) |> FindClusters(resolution = 0.6, verbose = FALSE)
root <- names(which.max(prop.table(table(obj$seurat_clusters, obj$broad), 1)[, "LTHSC"]))
sds  <- slingshot(Embeddings(obj, "umap"), clusterLabels = as.character(obj$seurat_clusters),
                  start.clus = root)

pt <- slingPseudotime(sds)                    # cells x lineages
um <- Embeddings(obj, "umap")
cells <- data.frame(UMAP1 = um[, 1], UMAP2 = um[, 2],
                    pseudotime = rowMeans(pt, na.rm = TRUE))
curves <- do.call(rbind, lapply(seq_along(slingCurves(sds)), function(i) {
  s <- slingCurves(sds)[[i]]$s[slingCurves(sds)[[i]]$ord, ]
  data.frame(UMAP1 = s[, 1], UMAP2 = s[, 2], lineage = paste0("Lineage ", i))
}))

range(pt, na.rm = TRUE)                        # the pseudotime scale
[1]  0.00000 16.40755
ggplot(cells, aes(UMAP1, UMAP2, color = pseudotime)) +
  geom_point(size = 0.8) +
  geom_path(data = curves, aes(group = lineage), color = "black", linewidth = 0.8) +
  scale_color_viridis_c(na.value = "grey85") +
  labs(title = "Pseudotime along the trajectory") +
  theme_minimal(base_size = 12)

The same UMAP of blood stem and progenitor cells, now coloured by pseudotime on a dark-to-bright viridis scale, with the four Slingshot lineage curves overlaid in black. The stem-cell root region is dark (low pseudotime) and the colour brightens smoothly outward along every branch to the committed-progenitor tips (high pseudotime), showing a continuous differentiation gradient rather than discrete steps.

Le pseudotemps va de 0 à la racine à environ 16 aux extrémités, et le gradient est lisse : aucun saut, juste une progression régulière vers l’extérieur le long de chaque ramification. Cette continuité est le modèle — les cellules sont ordonnées selon leur degré de différenciation, et non triées en stades.

Valider la trajectoire face aux étiquettes connues

Une trajectoire qui semble plausible ne suffit pas ; vérifiez-la face à la biologie que vous avez mise de côté. Croisez les regroupements avec les étiquettes FACS broad : le regroupement racine devrait être la population de cellules souches, et les regroupements d’extrémité devraient être les progéniteurs engagés.

set.seed(42)
library(Seurat)
library(slingshot)
d <- readRDS("_data/nestorowa_hsc.rds")
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
obj <- NormalizeData(obj, verbose = FALSE) |>
  FindVariableFeatures(nfeatures = 2000, verbose = FALSE) |> ScaleData(verbose = FALSE) |>
  RunPCA(npcs = 30, verbose = FALSE) |> RunUMAP(dims = 1:30, verbose = FALSE) |>
  FindNeighbors(dims = 1:30, verbose = FALSE) |> FindClusters(resolution = 0.6, verbose = FALSE)

# Cluster vs sorted population: does the root match LT-HSC and the tips match committed progenitors?
table(cluster = obj$seurat_clusters, FACS = obj$broad)
       FACS
cluster CMP GMP LMPP LTHSC MEP MPP
      0  39   9   94   191  15 209
      1 141  31   31    20 101  58
      2  19   3    0    18 216   2
      3  19   2  133    17   6  65
      4  76  71    2     6   3  12
      5   9   4   11    63   8  14
      6   3   1    9     5   3   5
      7  22   2    0     3  10   3

L’histoire tient. Le regroupement racine (5) est dominé par les cellules LTHSC — la véritable cellule souche au sommet — tandis que les extrémités sont des progéniteurs engagés : l’extrémité érythroïde/mégacaryocytaire (regroupement 2) est à ~84 % de MEP, l’extrémité myéloïde (regroupement 4) est la plus enrichie en GMP, et le tronc lymphoïde (regroupement 3) est à ~55 % de LMPP. Slingshot a placé les cellules souches au début et les cellules engagées aux extrémités sans jamais voir les étiquettes FACS — l’ordonnancement qu’il a récupéré de la seule expression correspond à la hiérarchie de différenciation connue. Cet accord est votre preuve que la trajectoire est réelle, pas un artefact de l’embedding.

Trouver les gènes qui changent le long de la trajectoire avec tradeSeq

Le pseudotemps est un moyen, pas une fin. Le gain biologique est quels gènes pilotent la différenciation — les gènes qui s’activent ou s’éteignent à mesure que les cellules parcourent chaque lignée. tradeSeq répond à cela : pour chaque gène il ajuste un modèle additif généralisé (GAM) lisse de l’expression contre le pseudotemps, par lignée, puis teste si le gène varie significativement le long de la trajectoire (Van den Berge et al., 2020).

fitGAM() fait l’ajustement et associationTest() fait le test. Ajuster un GAM par gène n’est pas gratuit, donc — comme les auteurs de tradeSeq le recommandent pour une première passe — nous ajustons les 150 gènes les plus variables avec 5 nœuds (la régularité de la courbe ajustée) ; une vraie analyse monte en échelle. Le pseudotemps et les poids de lignée par cellule proviennent tous deux directement de l’objet Slingshot.

set.seed(42)
library(Seurat)
library(slingshot)
library(tradeSeq)
d <- readRDS("_data/nestorowa_hsc.rds")
sym <- d$symbol
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
obj <- NormalizeData(obj, verbose = FALSE) |>
  FindVariableFeatures(nfeatures = 2000, verbose = FALSE) |> ScaleData(verbose = FALSE) |>
  RunPCA(npcs = 30, verbose = FALSE) |> RunUMAP(dims = 1:30, verbose = FALSE) |>
  FindNeighbors(dims = 1:30, verbose = FALSE) |> FindClusters(resolution = 0.6, verbose = FALSE)
root <- names(which.max(prop.table(table(obj$seurat_clusters, obj$broad), 1)[, "LTHSC"]))
sds  <- slingshot(Embeddings(obj, "umap"), clusterLabels = as.character(obj$seurat_clusters),
                  start.clus = root)

pt    <- slingPseudotime(sds, na = FALSE)      # cells x lineages, no NAs
cw    <- slingCurveWeights(sds)                # each cell's weight on each lineage
genes <- head(VariableFeatures(obj), 150)      # a fast first pass on the top 150 HVGs
set.seed(42)                                   # fitGAM initialization uses randomness — seed for a stable result
gam   <- fitGAM(counts = as.matrix(GetAssayData(obj, layer = "counts")[genes, ]),
                pseudotime = pt, cellWeights = cw, nknots = 5, verbose = FALSE)

assoc <- associationTest(gam)                  # test each gene for change along the lineages
assoc$symbol <- sym[rownames(assoc)]           # readable names
assoc <- assoc[order(assoc$pvalue), ]

sum(p.adjust(assoc$pvalue, "fdr") < 0.05, na.rm = TRUE)   # dynamic genes at 5% FDR
[1] 33
head(assoc[!is.na(assoc$symbol), c("symbol", "waldStat")], 10)
                    symbol waldStat
ENSMUSG00000051504 Siglech 127.2230
ENSMUSG00000105906   Iglc1 514.0457
ENSMUSG00000029373     Pf4 328.7270
ENSMUSG00000060586  H2-Eb1 206.4824
ENSMUSG00000036782  Klhl13 188.3286
ENSMUSG00000046826 Fam187b 123.0724
ENSMUSG00000061353  Cxcl12 196.5018
ENSMUSG00000038179  Slamf7 151.0669
ENSMUSG00000046818  Ddit4l 191.9911
ENSMUSG00000027524    Edn3 133.5374

Sur les 150 gènes testés, 33 varient significativement le long de la trajectoire à un taux de fausses découvertes de 5 % — des gènes dont l’expression est couplée à la différenciation plutôt que plate. Le haut de la liste est dominé par des gènes d’engagement de lignée et immunitaires : le gène d’immunoglobuline Iglc1 (marquant l’engagement B-lymphoïde), le facteur plaquettaire Pf4 (un gène d’engagement mégacaryocytaire), et le gène de présentation d’antigène H2-Eb1 (CMH classe II). Ce sont exactement les programmes que vous attendez de voir s’activer à mesure que des cellules multipotentes s’engagent vers des destins sanguins spécifiques.

Tracer les gènes dynamiques le long du pseudotemps

La façon la plus claire de voir un gène dynamique est de tracer son expression contre le pseudotemps, colorée par lignée. Choisissez un marqueur par destin — un gène érythroïde, un gène mégacaryocytaire, un gène myéloïde, un gène B-lymphoïde — et chacun devrait s’allumer dans sa propre lignée et rester silencieux dans les autres, la signature de l’engagement spécifique à un destin.

set.seed(42)
library(Seurat)
library(slingshot)
library(ggplot2)
d <- readRDS("_data/nestorowa_hsc.rds")
sym <- d$symbol
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
obj <- NormalizeData(obj, verbose = FALSE) |>
  FindVariableFeatures(nfeatures = 2000, verbose = FALSE) |> ScaleData(verbose = FALSE) |>
  RunPCA(npcs = 30, verbose = FALSE) |> RunUMAP(dims = 1:30, verbose = FALSE) |>
  FindNeighbors(dims = 1:30, verbose = FALSE) |> FindClusters(resolution = 0.6, verbose = FALSE)
root <- names(which.max(prop.table(table(obj$seurat_clusters, obj$broad), 1)[, "LTHSC"]))
sds  <- slingshot(Embeddings(obj, "umap"), clusterLabels = as.character(obj$seurat_clusters),
                  start.clus = root)

# Assign each cell to its strongest lineage and take that lineage's pseudotime
cw  <- slingCurveWeights(sds); pt <- slingPseudotime(sds, na = FALSE)
lin <- apply(cw, 1, which.max)
cell_pt <- pt[cbind(seq_len(nrow(pt)), lin)]

# One recognizable marker per fate (Ensembl IDs from the symbol map)
markers <- c(Slc4a1 = "ENSMUSG00000006574",  # erythroid
             Pf4    = "ENSMUSG00000029373",  # megakaryocyte
             Lgals3 = "ENSMUSG00000050335",  # myeloid
             Cd79a  = "ENSMUSG00000003379")  # B-lymphoid
logn <- GetAssayData(obj, layer = "data")
plotdf <- do.call(rbind, lapply(names(markers), function(g)
  data.frame(pseudotime = cell_pt, lineage = paste0("Lineage ", lin),
             expression = as.numeric(logn[markers[g], ]), gene = g)))

ggplot(plotdf, aes(pseudotime, expression, color = lineage)) +
  geom_point(size = 0.4, alpha = 0.25) +
  geom_smooth(se = FALSE, method = "loess", formula = y ~ x, linewidth = 0.9) +
  facet_wrap(~ gene, scales = "free_y") +
  scale_color_viridis_d(end = 0.9) +
  labs(x = "Pseudotime", y = "Log-normalized expression", color = NULL) +
  theme_minimal(base_size = 12)

A four-panel plot of gene expression against pseudotime, one panel per marker gene, with points and smoothed lines coloured by lineage. Slc4a1 (erythroid) rises late in the MEP-tipped lineage; Lgals3 (myeloid) rises in the GMP-tipped lineage; Pf4 and Cd79a each peak in a different branch. Different fate markers rise in different lineages, tracking the progenitor identity at each branch tip.

Chaque marqueur s’active dans une lignée différente, et les cas les plus nets suivent l’identité du progéniteur à l’extrémité de la ramification : le gène érythroïde Slc4a1 grimpe tard dans la lignée qui se termine dans le regroupement dominé par les MEP, et le gène myéloïde Lgals3 monte dans la lignée qui se termine dans le regroupement enrichi en GMP. Le gène mégacaryocytaire Pf4 et le gène B-lymphoïde Cd79a culminent chacun dans une ramification encore différente — un rappel que toute extrémité inférée n’est pas un destin engagé de manuel : deux de ces lignées se terminent dans des regroupements progéniteurs mixtes (un regroupement à pluralité CMP et un regroupement à orientation LMPP), donc les marqueurs rapportent un biais de lignée, pas un type cellulaire terminal net. C’est à la fois la puissance et la limite honnête de l’analyse de trajectoire : elle révèle non seulement que les cellules se différencient, mais quels gènes montent le long de chaque ramification et quand — aussi nettement que votre clustering résout les ramifications (voir les mises en garde ci-dessous).

Mises en garde honnêtes : ce que le pseudotemps est et n’est pas

L’inférence de trajectoire est puissante et facile à sur-interpréter. Gardez quatre choses au clair :

  • Le pseudotemps n’est pas le temps réel. Il ordonne les cellules par progression transcriptionnelle, pas par heures d’horloge. Deux cellules aux pseudotemps 5 et 10 diffèrent en maturité, mais rien ne dit que la seconde a pris deux fois plus de temps. Ne citez jamais le pseudotemps comme une durée.
  • Il faut un véritable continuum. Slingshot ajustera volontiers un arbre à n’importe quel ensemble de regroupements, y compris des types cellulaires discrets et sans lien — et vous rendra un ordonnancement dénué de sens. La méthode ne peut pas vous dire si une trajectoire existe ; c’est votre décision, prise avant de l’exécuter. Si votre UMAP est fait d’îlots bien séparés, pas d’une variété connectée, il n’y a aucune trajectoire à inférer.
  • La racine fixe la direction, et la choisir de travers renverse l’histoire. Slingshot récupère la topologie et l’ordonnancement le long d’elle, mais pas quelle extrémité est le début. Laissé à inférer la racine, il peut choisir le mauvais regroupement, inversant une lignée de sorte que des gènes semblent s’éteindre pendant la différenciation alors qu’ils s’allument. Ancrez la racine dans des connaissances réelles — ici l’étiquette FACS de cellule souche — chaque fois que vous le pouvez (start.clus).
  • Topologie et ordonnancement ne sont pas la même chose que la directionnalité. Même correctement enraciné, Slingshot infère la direction à partir de votre choix de racine, pas des données. Lire la direction dans les données elles-mêmes requiert la RNA velocity (à partir du rapport des transcrits non épissés sur les transcrits épissés) — une méthode d’abord Python (scVelo) sans équivalent R de première classe, couverte dans une leçon compagne sur la RNA velocity et la cartographie des destins. Une leçon compagne distincte couvre Monocle3, une alternative fondée sur les graphes qui apprend un graphe principal pour des topologies plus complexes ou déconnectées. Slingshot est la valeur par défaut propre et native R pour le cas courant de la ramification.

Quelle méthode quand

Vous avez… Utilisez Pourquoi
Une population en cours de différenciation avec une structure de ramification claire, en R Slingshot (+ tradeSeq) MST + courbes principales sur votre embedding existant ; simple, robuste, natif Seurat/SCE
Une topologie complexe ou déconnectée (boucles, nombreuses ramifications) Monocle3 (leçon compagne) Apprend un graphe principal, pas seulement un arbre ; son propre graph_test pour les gènes dynamiques
Un besoin d’inférer la direction à partir des données, pas d’une racine connue RNA velocity (scVelo, Python ; leçon compagne) Utilise les rapports épissés/non épissés pour pointer la flèche de la différenciation
Des types cellulaires discrets et stables (aucun continuum) Aucune — annotez plutôt Il n’y a aucun chemin de développement le long duquel ordonner ; utilisez l’annotation des types cellulaires

Le fil conducteur : Slingshot ordonne les cellules le long d’une trajectoire que vous avez jugée réelle ; tradeSeq trouve les gènes qui changent le long d’elle ; et la racine — issue de la biologie — fixe la direction.

Problèmes fréquents

Les lignées pointent dans le mauvais sens, ou la racine est sur un regroupement mature. Slingshot a inféré la racine et a deviné de travers. Il ordonne les cellules le long de l’arbre mais ne connaît pas la biologie — alors dites-le-lui : passez start.clus = le regroupement que vous savez être l’origine (ici le regroupement enrichi en LT-HSC, identifié à partir des étiquettes FACS). Avec la bonne racine, chaque lignée court de la cellule souche à l’engagée, et la dynamique des gènes se lit dans la bonne direction.

Trop ou trop peu de lignées. Le nombre de ramifications est piloté par le clustering, car le MST est construit sur les centroïdes des regroupements. Un sur-clustering (une résolution FindClusters élevée) fragmente un destin en plusieurs extrémités et invente des ramifications fallacieuses ; un sous-clustering fusionne des destins réels. Si le nombre de lignées ne correspond pas à la biologie connue, ajustez la résolution du clustering avant d’accuser Slingshot — la trajectoire ne vaut que ce que valent les regroupements sous-jacents.

Le pseudotemps est NA pour certaines cellules. Une cellule qui appartient à une seule lignée n’a pas de pseudotemps sur les autres, donc slingPseudotime() y renvoie NA par conception. Pour une valeur par cellule (par ex. pour colorer une figure), soit moyennez à travers les lignées avec rowMeans(pt, na.rm = TRUE), soit assignez chaque cellule à sa lignée la plus forte avec slingCurveWeights() et prenez le pseudotemps de cette lignée. Passez na = FALSE à slingPseudotime() quand un outil en aval (comme tradeSeq) a besoin d’une matrice complète.

Questions fréquentes

Le pseudotemps est un nombre attribué à chaque cellule qui mesure jusqu’où elle a progressé le long d’une trajectoire de développement, inféré de la façon dont son expression génique change graduellement par rapport aux autres cellules. Les cellules au début de la différenciation reçoivent des valeurs basses et les cellules pleinement engagées des valeurs élevées. Il n’est pas mesuré en temps réel (d’où pseudo) ; il capture l’ordonnancement et la maturité, pas la durée. Utilisez-le quand vos cellules forment un véritable continuum, comme des cellules souches se différenciant en types matures.

Slingshot travaille en deux étapes. Il construit d’abord un arbre couvrant de poids minimal sur les centroïdes des regroupements dans votre embedding de faible dimension pour trouver la topologie de ramification — quel état cellulaire se connecte à quel autre. Puis il ajuste des courbes principales lisses à travers les cellules le long de chaque chemin de la racine à l’extrémité, et la position d’une cellule le long de sa courbe est son pseudotemps. Vous l’appelez avec slingshot(embedding, clusterLabels = clusters, start.clus = root), en lui donnant un embedding (comme UMAP ou PCA), vos étiquettes de regroupements, et le regroupement de départ.

Slingshot ajuste un arbre couvrant de poids minimal sur les regroupements et des courbes lisses — simple, robuste, natif R/Seurat, et idéal pour une différenciation à ramification claire. Monocle3 apprend un graphe principal et gère des topologies plus complexes ou déconnectées (boucles, nombreuses ramifications) avec son propre graph_test pour les gènes dynamiques. Pour une trajectoire simple en forme d’arbre, Slingshot est la valeur par défaut propre ; tournez-vous vers Monocle3 quand la topologie est véritablement complexe. Une leçon compagne couvre Monocle3 en profondeur.

Utilisez tradeSeq après Slingshot. Il ajuste un modèle additif généralisé (GAM) lisse de l’expression de chaque gène contre le pseudotemps par lignée avec fitGAM(), puis associationTest() teste quels gènes varient significativement le long de la trajectoire. Fournissez-lui le pseudotemps (slingPseudotime(sds, na = FALSE)) et les poids de lignée par cellule (slingCurveWeights(sds)) de l’objet Slingshot. Les gènes significatifs sont ceux dont l’expression est couplée à la différenciation.

Non — Slingshot récupère la topologie (l’arbre de ramification) et l’ordonnancement le long d’elle, mais il ne peut pas dire quelle extrémité est le début. Vous fournissez la direction en fixant le regroupement racine (start.clus) à partir de connaissances biologiques, comme un marqueur de cellule souche connu ou une étiquette FACS. Inférer la direction à partir des données elles-mêmes requiert la RNA velocity (scVelo), qui utilise le rapport des transcrits non épissés sur les transcrits épissés et est une méthode d’abord Python couverte dans une leçon compagne.

Testez vos connaissances

En partant de nestorowa_hsc.rds, construisez les regroupements et exécutez Slingshot deux fois : une fois enraciné au regroupement enrichi en LT-HSC (comme dans la leçon), et une fois enraciné au regroupement le plus enrichi en MEP (une extrémité engagée). Comparez slingLineages() et le pseudotemps de quelques cellules LT-HSC entre les deux. Qu’arrive-t-il à l’ordonnancement, et pourquoi le choix de la racine n’est-il pas un détail cosmétique ?

Le seul changement est start.clus. Trouvez le regroupement enrichi en MEP de la même façon que vous avez trouvé celui enrichi en LT-HSC : names(which.max(prop.table(table(obj$seurat_clusters, obj$broad), 1)[, "MEP"])). Puis comparez slingPseudotime(sds_stem) à slingPseudotime(sds_mep) pour les cellules souches — une cellule qui était près du pseudotemps 0 dans un enracinement sera près du maximum dans l’autre.

set.seed(42)
library(Seurat); library(slingshot)
d <- readRDS("_data/nestorowa_hsc.rds")
obj <- CreateSeuratObject(counts = d$counts, meta.data = d$meta)
obj <- NormalizeData(obj, verbose = FALSE) |>
  FindVariableFeatures(nfeatures = 2000, verbose = FALSE) |> ScaleData(verbose = FALSE) |>
  RunPCA(npcs = 30, verbose = FALSE) |> RunUMAP(dims = 1:30, verbose = FALSE) |>
  FindNeighbors(dims = 1:30, verbose = FALSE) |> FindClusters(resolution = 0.6, verbose = FALSE)
emb <- Embeddings(obj, "umap"); cl <- as.character(obj$seurat_clusters)
fac <- prop.table(table(obj$seurat_clusters, obj$broad), 1)

stem_root <- names(which.max(fac[, "LTHSC"]))
mep_root  <- names(which.max(fac[, "MEP"]))
sds_stem  <- slingshot(emb, clusterLabels = cl, start.clus = stem_root)
sds_mep   <- slingshot(emb, clusterLabels = cl, start.clus = mep_root)

slingLineages(sds_stem)     # lineages run stem -> committed
slingLineages(sds_mep)      # lineages now run committed -> ... (reversed / re-rooted)

Le ré-enracinement inverse la lecture de la trajectoire : les cellules que Slingshot avait placées au pseudotemps 0 avec la racine souche se retrouvent à l’extrémité la plus lointaine avec la racine MEP. Les courbes sont similaires, mais la direction — et donc le fait qu’un gène se lise comme s’activant ou s’éteignant pendant la différenciation — dépend entièrement de la racine. C’est pourquoi vous ancrez la racine dans la biologie réelle, et non sur la supposition par défaut de Slingshot.

A. Oui — Slingshot s’est exécuté sans erreur, donc l’ordonnancement est valide. B. Non — ce sont des types cellulaires discrets et stables sans continuum de développement entre eux, donc il n’y a aucune trajectoire le long de laquelle ordonner ; Slingshot ajuste un arbre malgré tout, mais le pseudotemps est dénué de sens. C. Oui, mais seulement si vous fixez la racine correctement.

B. Slingshot ajustera toujours un arbre couvrant de poids minimal et renverra un pseudotemps — il ne peut pas vous dire si une trajectoire existe. Les monocytes, les lymphocytes T et les lymphocytes B sont des lignées matures séparées, pas des points sur un seul chemin de différenciation, donc les ordonner le long d’un pseudotemps est dénué de sens quelle que soit la façon dont vous l’enracinez. L’analyse de trajectoire est faite pour de véritables continuums ; décider que vous en avez un est un jugement que vous portez avant d’exécuter la méthode, en cherchant une variété connectée et ramifiée plutôt que des îlots séparés.

Conclusion

Une partie de la biologie est un gradient, et le clustering l’aplatit. L’analyse de trajectoire récupère le gradient : Slingshot construit un arbre couvrant de poids minimal sur vos regroupements et le lisse en courbes de lignées ramifiées, et la position de chaque cellule le long de sa courbe est son pseudotemps. Sur le jeu de données de sang Nestorowa il a résolu quatre lignées s’évasant d’une racine unique de cellules souches — et le recoupement avec les étiquettes FACS mises de côté a confirmé que la racine est bien la population LT-HSC et que les extrémités sont bien des progéniteurs engagés. Ensuite tradeSeq a trouvé les 33 gènes qui changent le long de la trajectoire, et tracer les marqueurs de destin contre le pseudotemps a montré chaque programme sanguin — érythroïde, mégacaryocytaire, myéloïde, lymphoïde — s’activant dans sa propre lignée. La discipline est ce qui compte : confirmez que vous avez un continuum avant de commencer, enracinez la trajectoire dans la biologie réelle pour que la direction soit correcte, et rappelez-vous que le pseudotemps est un ordonnancement, pas une horloge — et que lire la direction à partir des données elles-mêmes est le travail de la RNA velocity, pas de Slingshot.

Leçons connexes

🟢 Avec un agent IA

Demandez à Prova « comment exécuter le pseudotemps Slingshot sur mon propre objet Seurat et fixer le regroupement racine, puis trouver les gènes qui changent le long de la trajectoire avec tradeSeq ? » — elle répond avec du code R que vous pouvez exécuter sur vos propres données : regrouper et calculer l’embedding, appeler slingshot() avec start.clus, lire les lignées, et ajuster tradeSeq. The runtime is the judge. Demander à Prova →

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 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.

Références

Réutilisation

Citation

BibTeX
@online{2026,
  author = {},
  title = {Analyse de pseudotemps et de trajectoire en R avec Slingshot},
  date = {2026-07-03},
  url = {https://www.datanovia.com/learn/bioinformatics/single-cell/scrnaseq-trajectory-slingshot.html},
  langid = {fr}
}
Veuillez citer ce travail comme suit :
“Analyse de pseudotemps et de trajectoire en R avec Slingshot.” 2026. July 3. https://www.datanovia.com/learn/bioinformatics/single-cell/scrnaseq-trajectory-slingshot.html.