Trajectory inference using Slingshot

1. load libraries

2. Nice function to easily draw a graph:

3. Reading data

All_samples_Merged <- readRDS("../../All_samples_Merged_with_STCAT_and_renamed_FINAL.rds")

Define some color palette

# Define some color palette
pal <- c(scales::hue_pal()(8), RColorBrewer::brewer.pal(9, "Set1"), RColorBrewer::brewer.pal(8, "Set2"))
set.seed(1)
pal <- rep(sample(pal, length(pal)), 200)

4. Trajectory inference using Slingshot

             used   (Mb) gc trigger    (Mb)   max used   (Mb)
Ncells    8836312  472.0   15780645   842.8   11849183  632.9
Vcells 1232081946 9400.1 1529471188 11669.0 1243595250 9487.9

Let’s visualize which clusters we have in our dataset:

vars <- c("Patient_origin", "orig.ident", "seurat_clusters", "Phase")
pl <- list()

for (i in vars) {
  pl[[i]] <- DimPlot(obj, group.by = i, label = T) + theme_void() + NoLegend()
}
wrap_plots(pl)

table(obj$seurat_clusters)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13 
6789 5275 4663 4661 4086 3634 3536 3409 3338 3273 3212 1675 1063  691 

5. Exploring the data

vars <- c("TOX", "LAG3", "CTLA4", "TIGIT", "FOXP3", "IL2RA", "CTLA4", "IKZF2", "TIGIT", "CCR7", "SELL", "IL7R", "TCF7", "CCR7", "SELL", "IL7R", "TCF7", "LEF1", "GZMB", "PRF1", "IFNG", "KLRG1", "CD69", "CXCR6", "PRF1", "GZMB", "NKG7", "GNLY", "ISG15", "IFI6", "IFIT3", "MX1", "MKI67", "TOP2A", "STAT3", "AHR", "CCR6", "BATF", "PTPRC", "GZMB", "BCL6", "ICOS", "HSPA1A", "ATF3", "IL2RA", "CD69", "HLA-DRA", "CD34")
pl <- list()

pl <- list(DimPlot(obj, group.by = "seurat_clusters", label = T) + theme_void() + NoLegend())
for (i in vars) {
  pl[[i]] <- FeaturePlot(obj, features = i, order = T) + theme_void() + NoLegend()
}
wrap_plots(pl)

NA
NA

. Exploring the data-2

# Load necessary libraries
library(Seurat)
library(patchwork)
library(glue)

# Define marker groups
marker_groups <- list(
  "CD4 Tex (Exhausted)" = c("TOX", "LAG3", "CTLA4", "TIGIT"),
  "CD4 Treg" = c("FOXP3", "IL2RA", "CTLA4", "IKZF2", "TIGIT"),
  "CD4 Tcm (Central Memory)" = c("CCR7", "SELL", "IL7R", "TCF7"),
  "CD4 Tn (Naive)" = c("CCR7", "SELL", "IL7R", "TCF7", "LEF1"),
  "CD4 Tem (Effector Memory)" = c("GZMB", "PRF1", "IFNG", "KLRG1"),
  "CD4 Trm (Tissue Resident)" = c("CD69", "CXCR6"),
  "CD4 Tc (Cytotoxic)" = c("PRF1", "GZMB", "NKG7", "GNLY"),
  "CD4 Tisg (IFN Signature)" = c("ISG15", "IFI6", "IFIT3", "MX1"),
  "CD4 Proliferation" = c("MKI67", "TOP2A"),
  "CD4 Th17" = c("STAT3", "AHR", "CCR6", "BATF"),
  "CD4 Temra (Effector Memory RA+)" = c("PTPRC", "GZMB"),
  "CD4 Tfh (Follicular Helper)" = c("BCL6", "ICOS"),
  "CD4 Tstr (Stress)" = c("HSPA1A", "ATF3"),
  "CD4 Activated" = c("IL2RA", "CD69", "HLA-DRA")
)

# OPTIONAL: Uncomment this to save all plots in a single PDF
# pdf("CD4_MarkerGroups_FeaturePlots.pdf", width = 12, height = 10)

# Main plotting loop
for (group_name in names(marker_groups)) {
  genes <- unique(marker_groups[[group_name]])
  message(glue("Plotting: {group_name}"))

  # Cleaner titles and better theme
  pl <- lapply(genes, function(gene) {
    FeaturePlot(obj, features = gene, order = TRUE) +
      ggtitle(gene) +  # shorter title
      theme_minimal(base_size = 10) +  # use minimal for spacing
      NoLegend()
  })

  # Combine plots
  combined_plot <- wrap_plots(pl, ncol = 2) +  # fewer columns = more space
                   plot_annotation(title = group_name)

  print(combined_plot)
}

# OPTIONAL: Uncomment if using pdf()
# dev.off()

. Exploring the data-3

# Required markers
progenitor_markers <- c("CD34", "KIT", "GATA2", "MKI67", "PROM1", "FLT3")

# Feature plots for each marker
FeaturePlot(obj, features = progenitor_markers, order = TRUE, ncol = 3)

NA
NA
NA
NA

6. compute the lineages on these dataset

set.seed(1)

# # Define lineage ends
ENDS <- c("1", "4", "8" , "0","13")

set.seed(1)
lineages <- as.SlingshotDataSet(getLineages(
  data           = obj@reductions$umap@cell.embeddings,
  clusterLabels  = obj$seurat_clusters,
  dist.method    = "mnn", # It can be: "simple", "scaled.full", "scaled.diag", "slingshot" or "mnn"
  end.clus       = ENDS, # You can also define the ENDS!
  start.clus     = "3"
)) # define where to START the trajectories


# IF NEEDED, ONE CAN ALSO MANULALLY EDIT THE LINEAGES, FOR EXAMPLE:
# sel <- sapply( lineages@lineages, function(x){rev(x)[1]} ) %in% ENDS
# lineages@lineages <- lineages@lineages[ sel ]
# names(lineages@lineages) <- paste0("Lineage",1:length(lineages@lineages))
# lineages


# Change the reduction to our "fixed" UMAP2d (FOR VISUALISATION ONLY)
lineages@reducedDim <- obj@reductions$umap@cell.embeddings

{
  plot(obj@reductions$umap@cell.embeddings, col = pal[obj$seurat_clusters], cex = .5, pch = 16)
  lines(lineages, lwd = 1, col = "black", cex = 2)
  text(centroids2d, labels = rownames(centroids2d), cex = 0.8, font = 2, col = "white")
}

Defining Principal Curves

# Define curves
curves <- as.SlingshotDataSet(getCurves(
  data          = lineages,
  thresh        = 1e-1,
  stretch       = 1e-1,
  allow.breaks  = F,
  approx_points = 1000
))

curves

class: SlingshotDataSet 

lineages: 7 
Lineage1: 3  5  7  9  2  6  13  
Lineage2: 3  5  7  9  2  6  8  
Lineage3: 3  5  7  9  4  
Lineage4: 3  5  7  0  
Lineage5: 3  5  7  11  
Lineage6: 3  5  7  12  
Lineage7: 3  10  1  

curves: 7 
Curve1: Length: 41.042  Samples: 23200.83
Curve2: Length: 28.57   Samples: 24480.98
Curve3: Length: 22.1    Samples: 16916.57
Curve4: Length: 26.704  Samples: 21446.99
Curve5: Length: 26.124  Samples: 13694.79
Curve6: Length: 27.465  Samples: 13510.43
Curve7: Length: 17.775  Samples: 13465.54

Plot

# Plots
{
  plot(obj@reductions$umap@cell.embeddings, col = pal[obj$seurat_clusters], pch = 16)
  lines(curves, lwd = 2, col = "black")
  text(centroids2d, labels = levels(obj$seurat_clusters), cex = 1, font = 2)
}

compute the differentiation pseudotime

pseudotime <- slingPseudotime(curves, na = FALSE)
cellWeights <- slingCurveWeights(curves)

x <- rowMeans(pseudotime)
x <- x / max(x)
o <- order(x)

{
  plot(obj@reductions$umap@cell.embeddings[o, ],
    main = paste0("pseudotime"), pch = 16, cex = 0.4, axes = F, xlab = "", ylab = "",
    col = colorRampPalette(c("grey70", "orange3", "firebrick", "purple4"))(99)[x[o] * 98 + 1]
  )
  points(centroids2d, cex = 2.5, pch = 16, col = "#FFFFFF99")
  text(centroids2d, labels = levels(obj$seurat_clusters), cex = 1, font = 2)
}

7. Finding differentially expressed genes

# 🔴 CHANGE: Use all cells
sel_cells <- colnames(obj)

# 🔴 CHANGE: Use existing HVGs
hvg_genes <- VariableFeatures(obj)

library(BiocParallel)
# Use raw counts (needed for tradeSeq fitGAM)
counts_mat <- GetAssayData(obj, assay = "RNA", slot = "counts")


pseudotime_sel <- pseudotime[sel_cells, , drop = FALSE]

all(sel_cells %in% rownames(pseudotime))    # should be TRUE

[1] TRUE

all(sel_cells %in% rownames(cellWeights))  # should be TRUE

[1] TRUE

sceGAM <- fitGAM(
  counts = drop0(counts_mat[hvg_genes, sel_cells]),
  pseudotime = pseudotime_sel,
  cellWeights = cellWeights[sel_cells, , drop = FALSE],
  nknots = 7,
  verbose = TRUE,
  parallel = FALSE,  # start serially to avoid memory issues
  sce = TRUE
)

  |                                                  | 0 % ~calculating  
  |+                                                 | 1 % ~17h 10m 34s  
  |+                                                 | 2 % ~13h 37m 55s  
  |++                                                | 3 % ~12h 57m 50s  
  |++                                                | 4 % ~11h 46m 57s  
  |+++                                               | 5 % ~10h 48m 20s  
  |+++                                               | 6 % ~10h 26m 10s  
  |++++                                              | 7 % ~09h 50m 33s  
  |++++                                              | 8 % ~10h 12m 47s  
  |+++++                                             | 9 % ~09h 48m 23s  
  |+++++                                             | 10% ~09h 32m 08s  
  |++++++                                            | 11% ~09h 18m 31s  
  |++++++                                            | 12% ~09h 00m 59s  
  |+++++++                                           | 13% ~08h 45m 02s  
  |+++++++                                           | 14% ~08h 31m 56s  
  |++++++++                                          | 15% ~08h 17m 02s  
  |++++++++                                          | 16% ~08h 04m 38s  
  |+++++++++                                         | 17% ~07h 53m 17s  
  |+++++++++                                         | 18% ~07h 43m 17s  
  |++++++++++                                        | 19% ~07h 33m 30s  
  |++++++++++                                        | 20% ~07h 25m 29s  
  |+++++++++++                                       | 21% ~07h 15m 59s  
  |+++++++++++                                       | 22% ~07h 05m 33s  
  |++++++++++++                                      | 23% ~06h 56m 44s  
  |++++++++++++                                      | 24% ~06h 48m 31s  
  |+++++++++++++                                     | 25% ~06h 40m 41s  
  |+++++++++++++                                     | 26% ~06h 32m 36s  
  |++++++++++++++                                    | 27% ~06h 26m 28s  
  |++++++++++++++                                    | 28% ~06h 19m 01s  
  |+++++++++++++++                                   | 29% ~06h 10m 51s  
  |+++++++++++++++                                   | 30% ~06h 03m 35s  
  |++++++++++++++++                                  | 31% ~05h 56m 42s  
  |++++++++++++++++                                  | 32% ~05h 49m 59s  
  |+++++++++++++++++                                 | 33% ~05h 44m 32s  
  |+++++++++++++++++                                 | 34% ~05h 38m 33s  
  |++++++++++++++++++                                | 35% ~05h 32m 03s  
  |++++++++++++++++++                                | 36% ~05h 26m 07s  
  |+++++++++++++++++++                               | 37% ~05h 20m 02s  
  |+++++++++++++++++++                               | 38% ~05h 13m 43s  
  |++++++++++++++++++++                              | 39% ~05h 08m 29s  
  |++++++++++++++++++++                              | 40% ~05h 01m 40s  
  |+++++++++++++++++++++                             | 41% ~04h 55m 33s  
  |+++++++++++++++++++++                             | 42% ~04h 49m 10s  
  |++++++++++++++++++++++                            | 43% ~04h 43m 31s  
  |++++++++++++++++++++++                            | 44% ~04h 37m 44s  
  |+++++++++++++++++++++++                           | 45% ~04h 31m 58s  
  |+++++++++++++++++++++++                           | 46% ~04h 25m 55s  
  |++++++++++++++++++++++++                          | 47% ~04h 20m 30s  
  |++++++++++++++++++++++++                          | 48% ~04h 14m 41s  
  |+++++++++++++++++++++++++                         | 49% ~04h 09m 51s  
  |+++++++++++++++++++++++++                         | 50% ~04h 04m 17s  
  |++++++++++++++++++++++++++                        | 51% ~03h 58m 55s  
  |++++++++++++++++++++++++++                        | 52% ~03h 53m 48s  
  |+++++++++++++++++++++++++++                       | 53% ~03h 48m 30s  
  |+++++++++++++++++++++++++++                       | 54% ~03h 43m 03s  
  |++++++++++++++++++++++++++++                      | 55% ~03h 37m 37s  
  |++++++++++++++++++++++++++++                      | 56% ~03h 32m 17s  
  |+++++++++++++++++++++++++++++                     | 57% ~03h 27m 06s  
  |+++++++++++++++++++++++++++++                     | 58% ~03h 21m 49s  
  |++++++++++++++++++++++++++++++                    | 59% ~03h 16m 41s  
  |++++++++++++++++++++++++++++++                    | 60% ~03h 11m 23s  
  |+++++++++++++++++++++++++++++++                   | 61% ~03h 06m 18s  
  |+++++++++++++++++++++++++++++++                   | 62% ~03h 01m 13s  
  |++++++++++++++++++++++++++++++++                  | 63% ~02h 56m 37s  
  |++++++++++++++++++++++++++++++++                  | 64% ~02h 51m 23s  
  |+++++++++++++++++++++++++++++++++                 | 65% ~02h 46m 42s  
  |+++++++++++++++++++++++++++++++++                 | 66% ~02h 41m 47s  
  |++++++++++++++++++++++++++++++++++                | 67% ~02h 37m 09s  
  |++++++++++++++++++++++++++++++++++                | 68% ~02h 32m 06s  
  |+++++++++++++++++++++++++++++++++++               | 69% ~02h 27m 07s  
  |+++++++++++++++++++++++++++++++++++               | 70% ~02h 22m 07s  
  |++++++++++++++++++++++++++++++++++++              | 71% ~02h 17m 05s  
  |++++++++++++++++++++++++++++++++++++              | 72% ~02h 12m 06s  
  |+++++++++++++++++++++++++++++++++++++             | 73% ~02h 07m 04s  
  |+++++++++++++++++++++++++++++++++++++             | 74% ~02h 02m 11s  
  |++++++++++++++++++++++++++++++++++++++            | 75% ~01h 57m 29s  
  |++++++++++++++++++++++++++++++++++++++            | 76% ~01h 52m 37s  
  |+++++++++++++++++++++++++++++++++++++++           | 77% ~01h 47m 47s  
  |+++++++++++++++++++++++++++++++++++++++           | 78% ~01h 42m 52s  
  |++++++++++++++++++++++++++++++++++++++++          | 79% ~01h 38m 04s  
  |++++++++++++++++++++++++++++++++++++++++          | 80% ~01h 33m 16s  
  |+++++++++++++++++++++++++++++++++++++++++         | 81% ~01h 28m 27s  
  |+++++++++++++++++++++++++++++++++++++++++         | 82% ~01h 23m 43s  
  |++++++++++++++++++++++++++++++++++++++++++        | 83% ~01h 18m 58s  
  |++++++++++++++++++++++++++++++++++++++++++        | 84% ~01h 14m 13s  
  |+++++++++++++++++++++++++++++++++++++++++++       | 85% ~01h 09m 32s  
  |+++++++++++++++++++++++++++++++++++++++++++       | 86% ~01h 04m 45s  
  |++++++++++++++++++++++++++++++++++++++++++++      | 87% ~01h 00m 07s  
  |++++++++++++++++++++++++++++++++++++++++++++      | 88% ~55m 24s      
  |+++++++++++++++++++++++++++++++++++++++++++++     | 89% ~50m 46s      
  |+++++++++++++++++++++++++++++++++++++++++++++     | 90% ~46m 12s      
  |++++++++++++++++++++++++++++++++++++++++++++++    | 91% ~41m 33s      
  |++++++++++++++++++++++++++++++++++++++++++++++    | 92% ~36m 57s      
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 93% ~32m 18s      
  |+++++++++++++++++++++++++++++++++++++++++++++++   | 94% ~27m 38s      
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 95% ~23m 01s      
  |++++++++++++++++++++++++++++++++++++++++++++++++  | 96% ~18m 24s      
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~13m 48s      
  |+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~09m 11s      
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~04m 35s      
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07h 38m 38s

# Visualize gene counts per lineage
plotGeneCount(curves, clusters = obj$seurat_clusters, models = sceGAM)

# Plot lineages
lc <- sapply(lineages@lineages, function(x) {rev(x)[1]})
names(lc) <- gsub("Lineage", "L", names(lc))
lc.idx = match(lc, levels(obj$seurat_clusters))

plot(obj@reductions$umap@cell.embeddings, col = pal[obj$seurat_clusters], pch = 16)
lines(curves, lwd = 2, col = "black")

points(centroids2d[lc.idx, ], col = "black", pch = 16, cex = 4)
text(centroids2d[lc.idx, ], labels = names(lc), cex = 1, font = 2, col = "white")

Genes that change with pseudotime

gc()

             used    (Mb) gc trigger    (Mb)    max used    (Mb)
Ncells    9541130   509.6   19618421  1047.8    18617056   994.3
Vcells 1659947594 12664.4 9093818920 69380.4 11367270202 86725.4

set.seed(8)
res <- na.omit(associationTest(sceGAM, contrastType = "consecutive"))
res <- res[res$pvalue < 1e-3, ]
res <- res[res$waldStat > mean(res$waldStat), ]
res <- res[order(res$waldStat, decreasing = T), ]
res[1:10, ]

NA

We can plot their expression

par(mfrow = c(4, 4), mar = c(.1, .1, 2, 1))
{
  plot(obj@reductions$umap@cell.embeddings, col = pal[obj$seurat_clusters], cex = .5, pch = 16, axes = F, xlab = "", ylab = "")
  lines(curves, lwd = 2, col = "black")
  points(centroids2d[lc.idx, ], col = "black", pch = 15, cex = 3, xpd = T)
  text(centroids2d[lc.idx, ], labels = names(lc), cex = 1, font = 2, col = "white", xpd = T)
}

vars <- rownames(res[1:15, ])
vars <- na.omit(vars[vars != "NA"])

for (i in vars) {
  x <- drop0(obj@assays$SCT@data)[i, ]
  x <- (x - min(x)) / (max(x) - min(x))
  o <- order(x)
  plot(obj@reductions$umap@cell.embeddings[o, ],
    main = paste0(i), pch = 16, cex = 0.5, axes = F, xlab = "", ylab = "",
    col = colorRampPalette(c("lightgray", "grey60", "navy"))(99)[x[o] * 98 + 1]
  )
}

Genes that change between two pseudotime points

res <- na.omit(startVsEndTest(sceGAM, pseudotimeValues = c(0, 1)))
res <- res[res$pvalue < 1e-3, ]
res <- res[res$waldStat > mean(res$waldStat), ]
res <- res[order(res$waldStat, decreasing = T), ]
res[1:10, 1:6]

identify which genes go up or down. Let’s check lineage 1:

# Get the top UP and Down regulated in lineage 1
res_lin1 <- sort(setNames(res$logFClineage1, rownames(res)))
vars <- names(c(rev(res_lin1)[1:7], res_lin1[1:8]))
vars <- na.omit(vars[vars != "NA"])

par(mfrow = c(4, 4), mar = c(.1, .1, 2, 1))

{
  plot(obj@reductions$umap@cell.embeddings, col = pal[obj$seurat_clusters], cex = .5, pch = 16, axes = F, xlab = "", ylab = "")
  lines(curves, lwd = 2, col = "black")
  points(centroids2d[lc.idx, ], col = "black", pch = 15, cex = 3, xpd = T)
  text(centroids2d[lc.idx, ], labels = names(lc), cex = 1, font = 2, col = "white", xpd = T)
}

for (i in vars) {
  x <- drop0(obj@assays$SCT@data)[i, ]
  x <- (x - min(x)) / (max(x) - min(x))
  o <- order(x)
  plot(obj@reductions$umap@cell.embeddings[o, ],
    main = paste0(i), pch = 16, cex = 0.5, axes = F, xlab = "", ylab = "",
    col = colorRampPalette(c("lightgray", "grey60", "navy"))(99)[x[o] * 98 + 1]
  )
}

Genes that are different between lineages

res <- na.omit(diffEndTest(sceGAM))
res <- res[res$pvalue < 1e-3, ]
res <- res[res$waldStat > mean(res$waldStat), ]
res <- res[order(res$waldStat, decreasing = T), ]
res[1:10, ]

NA

pairwise comparison between each lineage. Let’s check lineage 1 vs lineage 2:

# Get the top UP and Down regulated in lineage 1 vs 2
res_lin1_2 <- sort(setNames(res$logFC1_2, rownames(res)))
vars <- names(c(rev(res_lin1_2)[1:7], res_lin1_2[1:8]))
vars <- na.omit(vars[vars != "NA"])

par(mfrow = c(4, 4), mar = c(.1, .1, 2, 1))
{
  plot(obj@reductions$umap@cell.embeddings, col = pal[obj$seurat_clusters], cex = .5, pch = 16, axes = F, xlab = "", ylab = "")
  lines(curves, lwd = 2, col = "black")
  points(centroids2d[lc.idx, ], col = "black", pch = 15, cex = 3, xpd = T)
  text(centroids2d[lc.idx, ], labels = names(lc), cex = 1, font = 2, col = "white", xpd = T)
}

for (i in vars) {
  x <- drop0(obj@assays$SCT@data)[i, ]
  x <- (x - min(x)) / (max(x) - min(x))
  o <- order(x)
  plot(obj@reductions$umap@cell.embeddings[o, ],
    main = paste0(i), pch = 16, cex = 0.5, axes = F, xlab = "", ylab = "",
    col = colorRampPalette(c("lightgray", "grey60", "navy"))(99)[x[o] * 98 + 1]
  )
}

save sceGAM after fitting:

saveRDS(sceGAM, "sceGAM_allHVGs.rds")

save sceGAM after fitting:

saveRDS(obj, "All_samples_Merged_after_Slingshot.rds")