Trajectory inference using Slingshot using PAGA embeddings

1. load libraries

2. Nice function to easily draw a graph:

3. Reading data

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

# 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)

📌 Step 1: Load PAGA embedding and match cells

📌 Step 2: Add PAGA UMAP as a new reduction in Seurat

📌 Step 3: Add Leiden cluster info to metadata

4. Trajectory inference using Slingshot

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

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

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

table(obj$leiden)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17 
7140 5705 5289 4709 4389 3835 3702 3308 2904 1944 1497 1268 1262  956  623  500  247   27 

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 = "leiden", label = T) + theme_void() + NoLegend())
for (i in vars) {
  pl[[i]] <- FeaturePlot(obj, features = i, order = T) + theme_void() + NoLegend()
}

Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
Please use the `layer` argument instead.

wrap_plots(pl)

NA
NA

6. compute the lineages on these dataset

# Define lineage ends
ENDS <- c("2", "6", "8")

set.seed(1)
lineages <- as.SlingshotDataSet(getLineages(
  data           = obj@reductions$paga@cell.embeddings,
  clusterLabels  = obj$leiden,
  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

Warning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observationsWarning: 'k' capped at the number of observations

# 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$paga@cell.embeddings

{
  plot(obj@reductions$paga@cell.embeddings, col = pal[obj$leiden], 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 = 100
))


# Plots
{
  plot(obj@reductions$paga@cell.embeddings, col = pal[obj$leiden], pch = 16)
  lines(curves, lwd = 2, col = "black")
  text(centroids2d, labels = levels(obj$leiden), 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$paga@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$leiden), cex = 1, font = 2)
}

7. Finding differentially expressed genes

sel_cells <- split(colnames(obj@assays$SCT@data), obj$leiden)
sel_cells <- unlist(lapply(sel_cells, function(x) {
  set.seed(1)
  return(sample(x, 20))
}))

gv <- as.data.frame(na.omit(scran::modelGeneVar(obj@assays$SCT@data[, sel_cells])))

Warning: collapsing to unique 'x' values

gv <- gv[order(gv$bio, decreasing = T), ]
sel_genes <- sort(rownames(gv)[1:500])

path_file <- "/home/bioinfo/data/trajectory/seurat_scegam.rds"

# fetch_data is defined at the top of this document
sceGAM <- fitGAM(
  counts = drop0(obj@assays$SCT@data[sel_genes, sel_cells]),
  pseudotime = pseudotime[sel_cells, ],
  cellWeights = cellWeights[sel_cells, ],
  nknots = 5, verbose = TRUE, parallel = TRUE, sce = TRUE,
  BPPARAM = BiocParallel::MulticoreParam()
)

Warning: Impossible to place a knot at all endpoints.Increase the number of knots to avoid this issue.

  |                                                                                                  
  |                                                                                            |   0%
  |                                                                                                  
  |==                                                                                          |   3%
  |                                                                                                  
  |=====                                                                                       |   5%
  |                                                                                                  
  |=======                                                                                     |   8%
  |                                                                                                  
  |==========                                                                                  |  10%
  |                                                                                                  
  |============                                                                                |  13%
  |                                                                                                  
  |==============                                                                              |  16%
  |                                                                                                  
  |=================                                                                           |  18%
  |                                                                                                  
  |===================                                                                         |  21%
  |                                                                                                  
  |======================                                                                      |  23%
  |                                                                                                  
  |========================                                                                    |  26%
  |                                                                                                  
  |==========================                                                                  |  29%
  |                                                                                                  
  |=============================                                                               |  31%
  |                                                                                                  
  |===============================                                                             |  34%
  |                                                                                                  
  |=================================                                                           |  36%
  |                                                                                                  
  |====================================                                                        |  39%
  |                                                                                                  
  |======================================                                                      |  42%
  |                                                                                                  
  |=========================================                                                   |  44%
  |                                                                                                  
  |===========================================                                                 |  47%
  |                                                                                                  
  |=============================================                                               |  49%
  |                                                                                                  
  |================================================                                            |  52%
  |                                                                                                  
  |==================================================                                          |  55%
  |                                                                                                  
  |=====================================================                                       |  57%
  |                                                                                                  
  |=======================================================                                     |  60%
  |                                                                                                  
  |=========================================================                                   |  62%
  |                                                                                                  
  |============================================================                                |  65%
  |                                                                                                  
  |==============================================================                              |  68%
  |                                                                                                  
  |=================================================================                           |  70%
  |                                                                                                  
  |==================================================================                          |  71%
  |                                                                                                  
  |====================================================================                        |  74%
  |                                                                                                  
  |======================================================================                      |  77%
  |                                                                                                  
  |=========================================================================                   |  79%
  |                                                                                                  
  |===========================================================================                 |  82%
  |                                                                                                  
  |==============================================================================              |  84%
  |                                                                                                  
  |================================================================================            |  87%
  |                                                                                                  
  |==================================================================================          |  90%
  |                                                                                                  
  |=====================================================================================       |  92%
  |                                                                                                  
  |=======================================================================================     |  95%
  |                                                                                                  
  |==========================================================================================  |  97%
  |                                                                                                  
  |============================================================================================| 100%

plotGeneCount(curves, clusters = obj$leiden, models = sceGAM)

lineages

class: SlingshotDataSet 

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

curves: 0 

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

{
  plot(obj@reductions$paga@cell.embeddings, col = pal[obj$leiden], 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")
}

NA
NA
NA
NA

Genes that change with pseudotime

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$paga@cell.embeddings, col = pal[obj$leiden], 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$paga@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]
  )
}

Error in .subscript.2ary(x, i, , drop = TRUE) : subscript out of bounds

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$paga@cell.embeddings, col = pal[obj$leiden], 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$paga@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$paga@cell.embeddings, col = pal[obj$leiden], 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$paga@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]
  )
}