Tutorial 14 — Trajectory & Cell–Cell Communication

Pseudotime (Slingshot/Monocle3), RNA velocity (scVelo), and ligand–receptor inference (CellChat)

Author

Single Cell RNA-seq Workshop

NoteRunning the code in this tutorial

Every code chunk here is tagged with #| eval: false, so the published page shows the code without running it. To run it yourself, open the downloaded .qmd in RStudio:

  • Run one chunk: click the green ▶ (Run Current Chunk) at the chunk’s top-right, or press Ctrl/Cmd + Enter. This runs the chunk regardless of its eval setting — the easiest way to work through the tutorial interactively.
  • Run a chunk on render: change that chunk’s #| eval: false to #| eval: true (or delete the line) so it executes when the document is rendered.
  • Render the whole tutorial: click the Render button in RStudio (or run quarto render in a terminal) to execute the chunks top-to-bottom and knit a finished HTML. To run everything on render, set eval: true once in the YAML header at the top.

About this tutorial

A bonus hands-on companion to Lecture 14 — Trajectory Inference & Cell–Cell Communication. These are optional downstream analyses; they build on the annotated, integrated object from the core workflow. We:

Note

Companion book chapter: Chapter 14 — Trajectory & Cell–Cell Communication — the long-form prose treatment of this tutorial’s material.

  • Order cells along a pseudotime with slingshot (and the monocle3 alternative)
  • Sketch an RNA velocity run with scVelo (Python), noting the spliced/unspliced input requirement
  • Predict ligand–receptor signaling between cell types with CellChat
  • Read every result with the appropriate skepticism
WarningThese analyses are easy to over-interpret

ifnb (resting + IFN-β-stimulated PBMCs) is not a developmental continuum, so a pseudotime on it is a mechanics demo, not biology. Trajectory tools always return a trajectory — only run one when a continuum is biologically plausible (development, a maturation axis, a time-course).

TipHow to use this page

The rendered HTML shows the code but does not execute it. To run it:

  1. Download the .qmd source: Tutorial_14_Trajectory_CellCommunication.qmd. If your browser saves it as ….qmd.txt, drop the trailing .txt before opening in RStudio.
  2. Make sure you’ve completed Tutorial 04 — Reference Annotation (or any tutorial that writes an annotated object) first.
  3. Install the extra packages below; work through the chunks, flipping eval: true as you go.

Setup

# Packages are installed in Tutorial 00 (Setup → bonus modules) — load them here.
library(Seurat); library(slingshot); library(SingleCellExperiment)
library(tidyverse)

seu <- readRDS("../data/ifnb_annotated.rds")
DefaultAssay(seu) <- "RNA"

Part A — Pseudotime with Slingshot

slingshot fits smooth principal curves through your clusters in a reduced space. You give it a start cluster (the root) — here we pretend one cell type is the progenitor purely to demonstrate the mechanics.

sce <- as.SingleCellExperiment(seu)
sce <- slingshot(sce,
                 clusterLabels = "seurat_annotations",
                 reducedDim    = "UMAP",
                 start.clus    = "CD14 Mono")   # <- you must justify the root

# Pseudotime values (one column per inferred lineage)
pt <- slingPseudotime(sce)
head(pt)

# Plot lineage curves over the UMAP
library(grDevices)
plot(reducedDims(sce)$UMAP, pch = 16, cex = 0.5, col = "grey70")
lines(SlingshotDataSet(sce), lwd = 2, col = "black")
TipReading the output

The scatter plot shows every cell in UMAP space (grey dots), with the inferred Slingshot lineage curves drawn over the top as black lines. Each curve is a separate trajectory branch; cells are ordered left-to-right (or whichever UMAP direction) by increasing pseudotime from the start.clus root. The head(pt) table above it shows raw pseudotime values: one column per lineage, NA for cells not assigned to that branch. Cells at the start of a lineage have values near 0; cells at the terminus have the highest values. Because ifnb is not a developmental dataset, treat these curves as a mechanics demo rather than meaningful biology — on a real differentiation dataset you would colour each cell by its pseudotime value to see the gradient.

ImportantThink about it

The start.clus argument flips the direction of the whole pseudotime. What evidence (markers, prior biology, CytoTRACE) would you need before trusting “CD14 Mono” as a root for a real dataset?

Answer You’d want an independent reason the start cells are the least-differentiated: known progenitor markers expressed there, a CytoTRACE score that’s highest at that end, and ideally agreement from RNA velocity. Without that, the ordering is arbitrary — Slingshot will happily root anywhere you tell it.

Monocle3 alternative

library(monocle3)
cds <- as.cell_data_set(seu)                 # SeuratWrappers
cds <- cluster_cells(cds)
cds <- learn_graph(cds)
cds <- order_cells(cds)                       # interactive: pick the root node
plot_cells(cds, color_cells_by = "pseudotime")
TipReading the output

plot_cells() returns a UMAP where each cell is coloured by its Monocle3 pseudotime value: cells at the root (the node you interactively chose with order_cells()) are dark/purple and cells at the terminus are yellow/light. The trajectory graph (learned by learn_graph()) is drawn as a set of connected edges overlaid on the UMAP; branch points appear where lineages diverge. Unlike Slingshot, Monocle3 can handle tree-shaped trajectories with multiple branch points, making it better suited to complex developmental hierarchies. UMAP coordinates are inherited from the Seurat object, so the cell layout is the same as your standard Seurat UMAP.

Part B — RNA velocity (scVelo, Python sketch)

Velocity needs a spliced/unspliced count matrix, which the standard Cell Ranger filtered matrix does not contain. Generate it first with velocyto run or kb count --workflow lamanno, then:

Note

This Python step is also packaged as a downloadable script — 14_rna_velocity_scvelo.py (Materials → Module 14 → Script · py). On Talapas, submit it with the Python SLURM wrapper: sbatch run_python.sbatch 14_rna_velocity_scvelo.py.

import scvelo as scv, scanpy as sc
adata = scv.read("ifnb_velocyto.loom")        # spliced + unspliced layers
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata, mode="dynamical")      # dynamical model (Bergen 2020)
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding_stream(adata, basis="umap")
TipReading the output

velocity_embedding_stream overlays a streamline plot on the UMAP: each arrow or flow line points in the direction a cell is predicted to move in transcriptional state space, based on the ratio of unspliced to spliced RNA for each gene. Cells with many unspliced transcripts relative to spliced (i.e., genes that are being actively induced) have arrows pointing toward states that over-express those genes. Coherent arrow fields pointing consistently from one cluster to another support a directional trajectory; disorganized or circular arrows suggest the dataset lacks a clear velocity signal (common in steady-state or stimulation data like ifnb). A flat or near-zero field usually means the unspliced layer was too sparse or the dynamical model failed to converge.

Warning

If you only have the filtered matrix, you cannot run velocity — there is no unspliced layer to model. This is the most common reason a velocity run “fails” silently or returns flat arrows.

CytoTRACE cross-check

library(CytoTRACE)
cyto <- CytoTRACE(as.matrix(GetAssayData(seu, slot = "counts")))
seu$cytotrace <- cyto$CytoTRACE
FeaturePlot(seu, "cytotrace")   # higher = less differentiated (more genes detected)
TipReading the output

FeaturePlot colours each cell in UMAP space by its CytoTRACE score (0–1): high scores (yellow/light) indicate less-differentiated cells that express a larger number of genes, while low scores (purple/dark) indicate terminally differentiated cells with a more restricted transcriptome. Use this as an independent cross-check on your Slingshot root: if CytoTRACE scores are highest in the cluster you chose as start.clus, the root is consistent with transcriptome breadth. A mismatch (CytoTRACE peaks in a different cluster) suggests the root assignment may need revisiting. Note that CytoTRACE is sensitive to sequencing depth; cells with systematically more reads will score higher regardless of biology.

Part C — Cell–cell communication with CellChat

library(CellChat)
cc  <- createCellChat(object = seu, group.by = "seurat_annotations")
cc@DB <- CellChatDB.human                       # human L–R database
cc <- subsetData(cc)
cc <- identifyOverExpressedGenes(cc)
cc <- identifyOverExpressedInteractions(cc)
cc <- computeCommunProb(cc)                      # permutation-based probabilities
cc <- filterCommunication(cc, min.cells = 10)
cc <- aggregateNet(cc)

# Network: who signals to whom, weighted by strength
netVisual_circle(cc@net$weight, weight.scale = TRUE,
                 title.name = "Aggregated interaction strength")

# Drill into one pathway
netVisual_aggregate(cc, signaling = "MIF", layout = "circle")
TipReading the output

The first circle plot (from netVisual_circle) shows every cell type as a node around the ring; directed arcs between nodes represent predicted signaling, with arc width proportional to the aggregated interaction strength (number of interactions × probability). Thick arcs indicate strong predicted cross-talk; self-loops mean the cell type signals to itself (autocrine or juxtacrine). The second circle plot drills into a single pathway (here MIF), showing only the ligand–receptor pairs in that pathway. Because ifnb is a two-condition PBMC stimulation dataset rather than a tissue, any CellChat output here is a technical demonstration: real CellChat results should be interpreted alongside spatial transcriptomics or perturbation data to validate that the predicted sender and receiver cells are actually co-localised. netVisual_bubble() is an alternative that shows individual L–R pairs as bubbles sized by probability.

ImportantThink about it

CellChat reports a strong predicted interaction between two clusters. List three reasons this might not reflect real signaling in the tissue.

Answer
  1. The two cell types may never be spatially adjacent — co-expression ignores location. (2) mRNA ≠ secreted, active protein — the ligand may not be translated/secreted. (3) The database (human/mouse-curated) may mis-map or omit the real partners, especially for non-model systems. Validate with spatial data or perturbation before claiming the circuit.

See also