---
title: "Tutorial 14 — Trajectory & Cell–Cell Communication"
subtitle: "Pseudotime (Slingshot/Monocle3), RNA velocity (scVelo), and ligand–receptor inference (CellChat)"
author: "Single Cell RNA-seq Workshop"
format:
  html:
    toc: true
    toc-depth: 3
    code-fold: false
    code-overflow: wrap
    highlight-style: github
    embed-resources: true
execute:
  eval: false
  echo: true
  warning: false
  message: false
editor: visual
---

::: {.callout-note title="Running 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](../Lecture_Folder/Lecture_14_Trajectory_CellCommunication.html)**. These are *optional* downstream analyses; they build on the annotated, integrated object from the core workflow. We:

::: callout-note
**Companion book chapter:** [Chapter 14 — Trajectory & Cell–Cell Communication](../Resources_Folder/Chapter_14_Trajectory_CellCommunication.html) — 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

::: {.callout-warning title="These 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).
:::

::: {.callout-tip title="How to use this page"}
The rendered HTML shows the code but does **not** execute it. To run it:

1.  Download the `.qmd` source: <a href="Tutorial_14_Trajectory_CellCommunication.qmd.txt" download="Tutorial_14_Trajectory_CellCommunication.qmd">Tutorial_14_Trajectory_CellCommunication.qmd</a>. 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](Tutorial_04_Reference_Annotation.html)** (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

```{r}
#| label: M14-packages_installed_tutorial_00
#| eval: false
# 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.

```{r}
#| label: M14-pseudotime_values_one_column
#| eval: false
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")
```

::: {.callout-tip title="Reading 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.
:::

::: {.callout-important title="Think 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?

<details><summary>Answer</summary>
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.
</details>
:::

### Monocle3 alternative

```{r}
#| label: M14-library_monocle3
#| eval: false
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")
```

::: {.callout-tip title="Reading 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:

::: callout-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`.
:::

```{python}
#| label: M14-import_scvelo_as_scv
#| eval: false
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")
```

::: {.callout-tip title="Reading 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.
:::

::: callout-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

```{r}
#| label: M14-library_cytotrace
#| eval: false
library(CytoTRACE)
cyto <- CytoTRACE(as.matrix(GetAssayData(seu, slot = "counts")))
seu$cytotrace <- cyto$CytoTRACE
FeaturePlot(seu, "cytotrace")   # higher = less differentiated (more genes detected)
```

::: {.callout-tip title="Reading 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

```{r}
#| label: M14-network_who_signals_whom
#| eval: false
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")
```

::: {.callout-tip title="Reading 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.
:::

::: {.callout-important title="Think about it"}
CellChat reports a strong predicted interaction between two clusters. List three reasons this might **not** reflect real signaling in the tissue.

<details><summary>Answer</summary>
(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.
</details>
:::

## See also

-   [Lecture 14 — Trajectory Inference & Cell–Cell Communication](../Lecture_Folder/Lecture_14_Trajectory_CellCommunication.html)
-   Saelens *et al.* 2019 — trajectory-method benchmark; Bergen *et al.* 2020 — [scVelo](https://scvelo.org/)
-   Jin *et al.* 2021 — [CellChat](http://www.cellchat.org/); Armingol *et al.* 2021 — cell–cell communication review