Trajectory inference using Slingshot

Seurat Toolkit

Reconstructing developmental or differentiation pathways from individual cell gene expression profiles to understand cellular transitions and relationships.
Authors

Åsa Björklund

Paulo Czarnewski

Susanne Reinsbach

Roy Francis

Published

15-Feb-2024

Note

Code chunks run R commands unless otherwise specified.

1 Loading libraries

suppressPackageStartupMessages({
  library(Seurat)
  library(plotly)
  options(rgl.printRglwidget = TRUE)
  library(Matrix)
  library(sparseMatrixStats)
  library(slingshot)
  library(tradeSeq)
  library(patchwork)
})

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

Nice function to easily draw a graph:

# Add graph to the base R graphics plot
draw_graph <- function(layout, graph, lwd = 0.2, col = "grey") {
  res <- rep(x = 1:(length(graph@p) - 1), times = (graph@p[-1] - graph@p[-length(graph@p)]))
  segments(
    x0 = layout[graph@i + 1, 1], x1 = layout[res, 1],
    y0 = layout[graph@i + 1, 2], y1 = layout[res, 2], lwd = lwd, col = col
  )
}

2 Preparing data

If you have been using the Seurat, Bioconductor or Scanpy toolkits with your own data, you need to reach to the point where you have:

  • A dimensionality reduction on which to run the trajectory (for example: PCA, ICA, MNN, harmony, Diffusion Maps, UMAP)
  • The cell clustering information (for example: from Louvain, K-means)
  • A KNN/SNN graph (this is useful to inspect and sanity-check your trajectories)

We will be using a subset of a bone marrow dataset (originally containing about 100K cells) for this exercise on trajectory inference.

The bone marrow is the source of adult immune cells, and contains virtually all differentiation stages of cell from the immune system which later circulate in the blood to all other organs.

You can download the data:

# download pre-computed data if missing or long compute
fetch_data <- TRUE

# url for source and intermediate data
path_data <- "https://export.uppmax.uu.se/naiss2023-23-3/workshops/workshop-scrnaseq"
path_file <- "data/trajectory/trajectory_seurat_filtered.rds"
if (!dir.exists(dirname(path_file))) dir.create(dirname(path_file), recursive = TRUE)
if (!file.exists(path_file)) download.file(url = file.path(path_data, "trajectory/trajectory_seurat_filtered.rds"), destfile = path_file)

We already have pre-computed and subsetted the dataset (with 6688 cells and 3585 genes) following the analysis steps in this course. We then saved the objects, so you can use common tools to open and start to work with them (either in R or Python).

In addition there was some manual filtering done to remove clusters that are disconnected and cells that are hard to cluster, which can be seen in this script

3 Reading data

obj <- readRDS("data/trajectory/trajectory_seurat_filtered.rds")

# Calculate cluster centroids (for plotting the labels later)
mm <- sparse.model.matrix(~ 0 + factor(obj$clusters_use))
colnames(mm) <- levels(factor(obj$clusters_use))
centroids2d <- as.matrix(t(t(obj@reductions$umap@cell.embeddings) %*% mm) / Matrix::colSums(mm))

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

vars <- c("batches", "dataset", "clusters_use", "Phase")
pl <- list()

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

You can check, for example, the number of cells in each cluster:

table(obj$clusters)

  1   2   5   6   7   8   9  11  12  13  14  15  16  17  18  19  20  21  22  23 
128  71  90 160 147 120 160 130 132  78  90 150 140  76 141  90  98 149  90  10 
 25  26  27  28  29  32  33  34  35  36  37  38  41  43  44  45  46  47  49  50 
 56 154  98  76 125 150 150 146 150 148 135 128 145 134 110 149 140 113 132  85 
 52  53  54  55  57  58  59  60  61 
126 129  57 129 147 127 118 120 101

4 Exploring the data

It is crucial that you have some understanding of the dataset being analyzed. What are the clusters you see in your data and most importantly How are the clusters related to each other?. Well, let’s explore the data a bit. With the help of this table, write down which cluster numbers in your dataset express these key markers.

Marker Cell Type
Cd34 HSC progenitor
Ms4a1 B cell lineage
Cd3e T cell lineage
Ltf Granulocyte lineage
Cst3 Monocyte lineage
Mcpt8 Mast Cell lineage
Alas2 RBC lineage
Siglech Dendritic cell lineage
C1qc Macrophage cell lineage
Pf4 Megakaryocyte cell lineage 
vars <- c("Cd34", "Ms4a1", "Cd3e", "Ltf", "Cst3", "Mcpt8", "Alas2", "Siglech", "C1qc", "Pf4")
pl <- list()

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

Another way to better explore your data is to look in higher dimensions, to really get a sense for what is right or wrong. As mentioned in the dimensionality reduction exercises, here we ran UMAP with 3 dimensions.

Important

The UMAP needs to be computed to results in exactly 3 dimensions

Since the steps below are identical to both Seurat and Scran pipelines, we will extract the matrices from both, so it is clear what is being used where and to remove long lines of code used to get those matrices. We will use them all. Plot in 3D with Plotly:

df <- data.frame(obj@reductions$umap3d@cell.embeddings, variable = factor(obj$clusters_use))
colnames(df)[1:3] <- c("UMAP_1", "UMAP_2", "UMAP_3")
p_State <- plot_ly(df, x = ~UMAP_1, y = ~UMAP_2, z = ~UMAP_3, color = ~variable, colors = pal, size = .5)  %>% add_markers()
p_State
# to save interactive plot and open in a new tab
try(htmlwidgets::saveWidget(p_State, selfcontained = T, "data/trajectory/umap_3d_clustering_plotly.html"), silent = T)
utils::browseURL("data/trajectory/umap_3d_clustering_plotly.html")

We can now compute the lineages on these dataset.

# Define lineage ends
ENDS <- c("17", "27", "25", "16", "26", "53", "49")

set.seed(1)
lineages <- as.SlingshotDataSet(getLineages(
  data           = obj@reductions$umap3d@cell.embeddings,
  clusterLabels  = obj$clusters_use,
  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     = "34"
)) # 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$clusters_use], cex = .5, pch = 16)
  lines(lineages, lwd = 1, col = "black", cex = 2)
  text(centroids2d, labels = rownames(centroids2d), cex = 0.8, font = 2, col = "white")
}

Much better!

5 Defining Principal Curves

Once the clusters are connected, Slingshot allows you to transform them to a smooth trajectory using principal curves. This is an algorithm that iteratively changes an initial curve to better match the data points. It was developed for linear data. To apply it to single-cell data, slingshot adds two enhancements:

  • It will run principal curves for each ‘lineage’, which is a set of clusters that go from a defined start cluster to some end cluster
  • Lineages with a same set of clusters will be constrained so that their principal curves remain bundled around the overlapping clusters

Since the function getCurves() takes some time to run, we can speed up the convergence of the curve fitting process by reducing the amount of cells to use in each lineage. Ideally you could all cells, but here we had set approx_points to 300 to speed up. Feel free to adjust that for your dataset.

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

curves
class: SlingshotDataSet 

 Samples Dimensions
    5828          2

lineages: 7 
Lineage1: 34  18  36  33  55  59  44  60  58  29  8  43  47  49  
Lineage2: 34  18  11  15  46  9  1  2  5  13  28  17  
Lineage3: 34  18  11  15  35  7  32  6  54  25  
Lineage4: 34  18  11  15  35  7  32  6  27  
Lineage5: 34  18  36  21  12  20  16  
Lineage6: 34  18  36  33  55  38  53  
Lineage7: 34  18  36  26  

curves: 7 
Curve1: Length: 6.7241  Samples: 2161.05
Curve2: Length: 7.3487  Samples: 2097.02
Curve3: Length: 3.5349  Samples: 1502.25
Curve4: Length: 2.5623  Samples: 1387.66
Curve5: Length: 2.9268  Samples: 979.78
Curve6: Length: 2.8976  Samples: 1086.34
Curve7: Length: 2.1323  Samples: 644.86
# Plots
{
  plot(obj@reductions$umap@cell.embeddings, col = pal[obj$clusters_use], pch = 16)
  lines(curves, lwd = 2, col = "black")
  text(centroids2d, labels = rownames(centroids2d), cex = 1, font = 2)
}

With those results in hands, we can now 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 = rownames(centroids2d), cex = 1, font = 2)
}

Tip

The pseudotime represents the distance of every cell to the starting cluster!

6 Finding differentially expressed genes

The main way to interpret a trajectory is to find genes that change along the trajectory. There are many ways to define differential expression along a trajectory:

  • Expression changes along a particular path (i.e. change with pseudotime)
  • Expression differences between branches
  • Expression changes at branch points
  • Expression changes somewhere along the trajectory
  • …

tradeSeq is a recently proposed algorithm to find trajectory differentially expressed genes. It works by smoothing the gene expression along the trajectory by fitting a smoother using generalized additive models (GAMs), and testing whether certain coefficients are statistically different between points in the trajectory.

BiocParallel::register(BiocParallel::MulticoreParam())

The fitting of GAMs can take quite a while, so for demonstration purposes we first do a very stringent filtering of the genes.

Tip

In an ideal experiment, you would use all the genes, or at least those defined as being variable.

sel_cells <- split(colnames(obj@assays$RNA@data), obj$clusters_use)
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$RNA@data[, sel_cells])))
gv <- gv[order(gv$bio, decreasing = T), ]
sel_genes <- sort(rownames(gv)[1:500])

Fitting the model:

Caution

This is a slow compute intensive step, we will not run this now and instead use a pre-computed file in the step below.

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

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

Download the precomputed file.

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

# fetch_data is defined at the top of this document
if (fetch_data) {
  if (!file.exists(path_file)) download.file(url = file.path(path_data, "trajectory/results/seurat_scegam.rds"), destfile = path_file)
}
# read data
sceGAM <- readRDS(path_file)
plotGeneCount(curves, clusters = obj$clusters_use, models = sceGAM)

lineages
class: SlingshotDataSet 

 Samples Dimensions
    5828          2

lineages: 7 
Lineage1: 34  18  36  33  55  59  44  60  58  29  8  43  47  49  
Lineage2: 34  18  11  15  46  9  1  2  5  13  28  17  
Lineage3: 34  18  11  15  35  7  32  6  54  25  
Lineage4: 34  18  11  15  35  7  32  6  27  
Lineage5: 34  18  36  21  12  20  16  
Lineage6: 34  18  36  33  55  38  53  
Lineage7: 34  18  36  26  

curves: 0 
lc <- sapply(lineages@lineages, function(x) {
  rev(x)[1]
})
names(lc) <- gsub("Lineage", "L", names(lc))

{
  plot(obj@reductions$umap@cell.embeddings, col = pal[obj$clusters_use], pch = 16)
  lines(curves, lwd = 2, col = "black")
  points(centroids2d[lc, ], col = "black", pch = 16, cex = 4)
  text(centroids2d[lc, ], labels = names(lc), cex = 1, font = 2, col = "white")
}

6.1 Genes that change with pseudotime

We can first look at general trends of gene expression across 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, ]

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$clusters_use], cex = .5, pch = 16, axes = F, xlab = "", ylab = "")
  lines(curves, lwd = 2, col = "black")
  points(centroids2d[lc, ], col = "black", pch = 15, cex = 3, xpd = T)
  text(centroids2d[lc, ], 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$RNA@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]
  )
}

6.2 Genes that change between two pseudotime points

We can define custom pseudotime values of interest if we’re interested in genes that change between particular point in pseudotime. By default, we can look at differences between start and end:

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]

You can see now that there are several more columns, one for each lineage. This table represents the differential expression within each lineage, to 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$clusters_use], cex = .5, pch = 16, axes = F, xlab = "", ylab = "")
  lines(curves, lwd = 2, col = "black")
  points(centroids2d[lc, ], col = "black", pch = 15, cex = 3, xpd = T)
  text(centroids2d[lc, ], labels = names(lc), cex = 1, font = 2, col = "white", xpd = T)
}

for (i in vars) {
  x <- drop0(obj@assays$RNA@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]
  )
}

6.3 Genes that are different between lineages

More interesting are genes that are different between two branches. We may have seen some of these genes already pop up in previous analyses of pseudotime. There are several ways to define “different between branches”, and each have their own functions:

  • Different at the end points, using diffEndTest()
  • Different at the branching point, using earlyDETest()
  • Different somewhere in pseudotime the branching point, using patternTest()

Note that the last function requires that the pseudotimes between two lineages are aligned.

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, ]

You can see now that there are even more columns, one for the 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$clusters_use], cex = .5, pch = 16, axes = F, xlab = "", ylab = "")
  lines(curves, lwd = 2, col = "black")
  points(centroids2d[lc, ], col = "black", pch = 15, cex = 3, xpd = T)
  text(centroids2d[lc, ], labels = names(lc), cex = 1, font = 2, col = "white", xpd = T)
}

for (i in vars) {
  x <- drop0(obj@assays$RNA@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]
  )
}

Check out this vignette for a more in-depth overview of tradeSeq and many other differential expression tests.

7 Generating batch-corrected data for differential gene expression

Before computing differential gene expression, sometimes it is a good idea to make sure our dataset is somewhat homogeneous (without very strong batch effects). In this dataset, we actually used data from 4 different technologies (Drop-seq, SmartSeq2 and 10X) and therefore massive differences in read counts can be observed:

If you want to know more about how to control for this issue, please have a look at batch_corrected_counts.Rmd

8 References

Cannoodt, Robrecht, Wouter Saelens, and Yvan Saeys. 2016. “Computational Methods for Trajectory Inference from Single-Cell Transcriptomics.” European Journal of Immunology 46 (11): 2496–2506. doi.

Saelens, Wouter, Robrecht Cannoodt, Helena Todorov, and Yvan Saeys. 2019. “A Comparison of Single-Cell Trajectory Inference Methods.” Nature Biotechnology 37 (5): 547–54. doi.

9 Session info

Click here 
sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] patchwork_1.1.2             tradeSeq_1.16.0            
 [3] slingshot_2.10.0            TrajectoryUtils_1.10.0     
 [5] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
 [7] Biobase_2.62.0              GenomicRanges_1.54.1       
 [9] GenomeInfoDb_1.38.5         IRanges_2.36.0             
[11] S4Vectors_0.40.2            BiocGenerics_0.48.1        
[13] princurve_2.1.6             sparseMatrixStats_1.14.0   
[15] MatrixGenerics_1.14.0       matrixStats_1.0.0          
[17] Matrix_1.5-4                plotly_4.10.2              
[19] ggplot2_3.4.2               SeuratObject_4.1.3         
[21] Seurat_4.3.0               

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.20          splines_4.3.0            
  [3] later_1.3.1               bitops_1.0-7             
  [5] tibble_3.2.1              polyclip_1.10-4          
  [7] lifecycle_1.0.3           edgeR_4.0.7              
  [9] globals_0.16.2            lattice_0.21-8           
 [11] MASS_7.3-58.4             crosstalk_1.2.0          
 [13] magrittr_2.0.3            limma_3.58.1             
 [15] rmarkdown_2.22            yaml_2.3.7               
 [17] metapod_1.10.1            httpuv_1.6.11            
 [19] sctransform_0.3.5         sp_1.6-1                 
 [21] spatstat.sparse_3.0-1     reticulate_1.30          
 [23] cowplot_1.1.1             pbapply_1.7-0            
 [25] RColorBrewer_1.1-3        abind_1.4-5              
 [27] zlibbioc_1.48.0           Rtsne_0.16               
 [29] purrr_1.0.1               RCurl_1.98-1.12          
 [31] GenomeInfoDbData_1.2.11   ggrepel_0.9.3            
 [33] irlba_2.3.5.1             listenv_0.9.0            
 [35] spatstat.utils_3.0-3      goftest_1.2-3            
 [37] dqrng_0.3.0               spatstat.random_3.1-5    
 [39] fitdistrplus_1.1-11       parallelly_1.36.0        
 [41] DelayedMatrixStats_1.24.0 leiden_0.4.3             
 [43] codetools_0.2-19          DelayedArray_0.28.0      
 [45] scuttle_1.12.0            tidyselect_1.2.0         
 [47] farver_2.1.1              ScaledMatrix_1.10.0      
 [49] viridis_0.6.3             spatstat.explore_3.2-1   
 [51] jsonlite_1.8.5            BiocNeighbors_1.20.2     
 [53] ellipsis_0.3.2            progressr_0.13.0         
 [55] ggridges_0.5.4            survival_3.5-5           
 [57] tools_4.3.0               ica_1.0-3                
 [59] Rcpp_1.0.10               glue_1.6.2               
 [61] gridExtra_2.3             SparseArray_1.2.3        
 [63] xfun_0.39                 mgcv_1.8-42              
 [65] dplyr_1.1.2               withr_2.5.0              
 [67] fastmap_1.1.1             bluster_1.12.0           
 [69] fansi_1.0.4               digest_0.6.31            
 [71] rsvd_1.0.5                R6_2.5.1                 
 [73] mime_0.12                 colorspace_2.1-0         
 [75] scattermore_1.2           tensor_1.5               
 [77] spatstat.data_3.0-1       utf8_1.2.3               
 [79] tidyr_1.3.0               generics_0.1.3           
 [81] data.table_1.14.8         httr_1.4.6               
 [83] htmlwidgets_1.6.2         S4Arrays_1.2.0           
 [85] uwot_0.1.14               pkgconfig_2.0.3          
 [87] gtable_0.3.3              lmtest_0.9-40            
 [89] XVector_0.42.0            htmltools_0.5.5          
 [91] scales_1.2.1              png_0.1-8                
 [93] scran_1.30.0              knitr_1.43               
 [95] rstudioapi_0.14           reshape2_1.4.4           
 [97] nlme_3.1-162              zoo_1.8-12               
 [99] stringr_1.5.0             KernSmooth_2.23-20       
[101] parallel_4.3.0            miniUI_0.1.1.1           
[103] pillar_1.9.0              grid_4.3.0               
[105] vctrs_0.6.2               RANN_2.6.1               
[107] promises_1.2.0.1          BiocSingular_1.18.0      
[109] beachmat_2.18.0           xtable_1.8-4             
[111] cluster_2.1.4             evaluate_0.21            
[113] cli_3.6.1                 locfit_1.5-9.8           
[115] compiler_4.3.0            rlang_1.1.1              
[117] crayon_1.5.2              future.apply_1.11.0      
[119] labeling_0.4.2            plyr_1.8.8               
[121] stringi_1.7.12            viridisLite_0.4.2        
[123] deldir_1.0-9              BiocParallel_1.36.0      
[125] munsell_0.5.0             lazyeval_0.2.2           
[127] spatstat.geom_3.2-1       future_1.32.0            
[129] statmod_1.5.0             shiny_1.7.4              
[131] ROCR_1.0-11               igraph_1.4.3