Trajectory analysis (Monocle3)

load libraries

1. Load Seurat Object

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

2. Prepare Data

# Seurat object
sdata <- All_samples_Merged  

DefaultAssay(sdata) <- "RNA"

# Extract counts
data <- GetAssayData(sdata, assay = "RNA", slot = "counts")

# Metadata
cell_metadata <- sdata@meta.data

# Gene annotation
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)

# Create Monocle3 CDS
cds <- new_cell_data_set(data,
                         cell_metadata = cell_metadata,
                         gene_metadata = gene_annotation)

3. Preprocessing and Dimension Reduction

cds <- preprocess_cds(cds, num_dim = 20)
cds <- reduce_dimension(cds, preprocess_method = "PCA")

3.1 Use Seurat’s UMAP for Consistency

# Extract UMAP embeddings
cds.embed <- cds@int_colData$reducedDims$UMAP
int.embed <- Embeddings(sdata, reduction = "umap")
int.embed <- int.embed[rownames(cds.embed), ]

# Overwrite UMAP in CDS
cds@int_colData$reducedDims$UMAP <- int.embed


library(monocle3)

# Create a folder for the checkpoint
dir.create("cds_checkpoint_1", showWarnings = FALSE)

# Save everything
save_monocle_objects(cds=cds, directory_path='cds_checkpoint_1', comment='This is my example cds. Stored 2025-08-22.')

3.2 Transfer useful Seurat metadata to cds

# Common metadata fields you likely have:
md_to_copy <- c("seurat_clusters","Patient_origin","cell_line","orig.ident")

for (fld in md_to_copy) {
  if (fld %in% colnames(sdata@meta.data)) {
    cds[[fld]] <- sdata@meta.data[match(colnames(cds), colnames(sdata)), fld]
  }
}

3.3 Clustering, Graph Learning, and Pseudotime Ordering

# Cluster cells
cds <- cluster_cells(cds, reduction_method = "UMAP")

# Learn trajectory graph
cds <- learn_graph(cds)

  |                                                                                                    
  |                                                                                              |   0%
  |                                                                                                    
  |==============================================================================================| 100%

  |                                                                                                    
  |                                                                                              |   0%
  |                                                                                                    
  |==============================================================================================| 100%

# Identify naive T cell markers available in CDS
naive_markers <- c("CCR7", "SELL", "LEF1", "TCF7")
naive_markers <- naive_markers[naive_markers %in% rownames(cds)]

# Extract expression matrix using logcounts assay if present, else log-transform counts
if("logcounts" %in% assayNames(cds)) {
expr_mat <- assay(cds, "logcounts")
} else {
expr_mat <- log1p(assay(cds, "counts"))
}

#Calculate mean naive marker expression (naive score) per cell
naive_score <- Matrix::colMeans(expr_mat[naive_markers, , drop = FALSE])

#Select root cells as those with naive score above 95th percentile threshold
threshold <- quantile(naive_score, 0.95)
root_cells <- names(naive_score[naive_score > threshold])


# Get closest principal graph nodes for root cells
closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex

# Map from cells → node names
root_nodes <- igraph::V(principal_graph(cds)[["UMAP"]])$name[
  as.numeric(closest_vertex[root_cells])
]

# Order cells using valid root_pr_nodes
cds <- order_cells(cds, root_cells = root_cells)

save CDS

library(monocle3)

# Create a folder for the checkpoint
dir.create("cds_checkpoint_2", showWarnings = FALSE)

# Save everything
save_monocle_objects(cds=cds, directory_path='cds_checkpoint_2', comment='This is my example cds. Stored 2025-08-22.')

4. Visualization of Trajectory

p1 <- plot_cells(cds,
                color_cells_by = "pseudotime",
                label_cell_groups = FALSE,
                label_leaves = FALSE,
                label_branch_points = FALSE,
                graph_label_size = 1.5)

p1

p2 <- plot_cells(cds,
                color_cells_by = "Patient_origin",
                label_cell_groups = FALSE,
                label_leaves = FALSE,
                label_branch_points = FALSE,
                graph_label_size = 1.5)

p2

p3 <- plot_cells(cds,
                color_cells_by = "orig.ident",
                label_cell_groups = FALSE,
                label_leaves = FALSE,
                label_branch_points = FALSE,
                graph_label_size = 1.5)

p3

Visualize Gene Expression Along Pseudotime

library(monocle3)
library(ggplot2)
library(dplyr)        # optional but commonly used for data manipulation

plot_genes_in_pseudotime(cds["CCR7", ])

plot_genes_in_pseudotime(cds["SELL", ])

plot_genes_in_pseudotime(cds["LEF1", ])

plot_genes_in_pseudotime(cds["TCF7", ])

#Output tells you: how canonical naïve markers decline/shift across pseudotime.

Explore Branch-Specific Differential Expression

## Identify genes that vary along branches
deg_by_branch <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)
deg_by_branch <- deg_by_branch %>% arrange(q_value)

save CDS

library(monocle3)

# Create a folder for the checkpoint
dir.create("cds_checkpoint_3", showWarnings = FALSE)

# Save everything
save_monocle_objects(cds=cds, directory_path='cds_checkpoint_3', comment='This is my example cds. Stored 2025-08-22.')

Explore Branch-Specific Differential Expression

## Select top DEGs per branch (q_value < 0.05)
top_degs <- deg_by_branch %>%
  filter(q_value < 0.05) %>%
  arrange(q_value)

# Optional: take top 50 for plotting
top50_genes <- head(top_degs$gene_short_name, 50)

## Plot pseudotime expression of top DEGs
plot_cells(cds[top50_genes, ], 
           genes=top50_genes, 
           show_trajectory_graph = TRUE, 
           label_cell_groups = FALSE)

#Output: Highlights branch-specific transcriptional programs and cell fate biases.

Explore Branch-Specific Differential Expression (genes one by one)

## Identify genes that vary along branches
deg_by_branch <- graph_test(cds, neighbor_graph = "principal_graph", cores = 4)

## Select top DEGs per branch (q_value < 0.05)
top_degs <- deg_by_branch %>%
  filter(q_value < 0.05) %>%
  arrange(q_value)

# Optional: take top 50 for plotting
top50_genes <- head(top_degs$gene_short_name, 50)

library(ggplot2)

for (gene in top50_genes) {
  p <- plot_cells(cds[gene, ], 
                  genes = gene,
                  show_trajectory_graph = TRUE,
                  label_cell_groups = FALSE) +
       ggtitle(gene)
  print(p)
}

Summarize Lineages and Cell Fate Branches

cds <- cluster_cells(cds)
p_branches <- plot_cells(cds, color_cells_by = "cluster")
print(p_branches)

Visualize Gene Modules Along Pseudotime

## Identify co-expressed gene modules
gene_modules <- find_gene_modules(cds, resolution = 1e-2)

## Explore first module along trajectory
plot_cells(cds, 
           genes = gene_modules[[1]], 
           show_trajectory_graph = FALSE, 
           label_cell_groups = FALSE) +
  ggtitle("Gene Module 1 Expression Along Trajectory")

## Optional: loop over modules to visualize multiple
lapply(seq_along(gene_modules)[1:5], function(i) {
  plot_cells(cds, genes = gene_modules[[i]], show_trajectory_graph = FALSE)
})

Multi-Gene Pseudotime Plot with Facets

## Multi-Gene Pseudotime Plots by Functional Marker Groups

library(ggplot2)

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

# Loop through each marker group and plot
for (group_name in names(marker_groups)) {
  
  genes_to_plot <- marker_groups[[group_name]]
  # Keep only genes present in CDS
  genes_to_plot <- genes_to_plot[genes_to_plot %in% rownames(cds)]
  
  if (length(genes_to_plot) > 0) {
    p <- plot_cells(cds[genes_to_plot, ], 
                    genes = genes_to_plot, 
                    show_trajectory_graph = FALSE) +
         facet_wrap(~ gene_short_name, scales = "free_y") +
         ggtitle(paste0(group_name, " Marker Dynamics Along Pseudotime"))
    
    print(p)
  }
}

# ✅ Output:
#   One faceted pseudotime plot per functional T-cell program.
#   Each figure shows expression dynamics of the group’s markers across pseudotime.

Multi-Gene Pseudotime Plot with Facets

library(ggplot2)

# Example gene list
genes_to_plot 

# Keep only genes present in CDS
genes_to_plot <- genes_to_plot[genes_to_plot %in% rownames(cds)]

plot_cells(cds[genes_to_plot, ], 
           genes = genes_to_plot, 
           show_trajectory_graph = FALSE) +
  facet_wrap(~ gene_short_name, scales = "free_y") +
  ggtitle("Naive T Cell Marker Dynamics Along Pseudotime")

#Output: Clear visualization of each gene’s dynamics along pseudotime.

Optional: Overlay CITE-seq Protein Expression

# Example: if your Seurat object has CITE-seq (ADT) data
adt_markers <- c("adt_CD45RA", "adt_CD62L")  # adjust to your panel
for(marker in adt_markers){
  plot_cells(cds, 
             color_cells_by = marker, 
             show_trajectory_graph = TRUE) +
    ggtitle(paste("ADT Expression:", marker))
}

Pathway Enrichment of Top DEGs or Gene Modules

library(clusterProfiler)
library(org.Hs.eg.db)

# Convert gene symbols to Entrez IDs
entrez_ids <- bitr(top50_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")

# KEGG enrichment
kegg_res <- enrichKEGG(entrez_ids$ENTREZID, organism="hsa")

# Visualize top pathways
dotplot(kegg_res, showCategory=10) + ggtitle("KEGG Pathways for Top DEGs")

# Output: Biological interpretation of branches (e.g., proliferation, T-cell activation, immune evasion).

Branch Probability / Fate Bias

library(igraph)

# Identify branch points
branch_nodes <- branch_nodes(cds)

# Map each cell to branch node
cell_branches <- sapply(branch_nodes, function(node){
  igraph::shortest_paths(principal_graph(cds)[["UMAP"]], 
                         from = node, to = colnames(cds))$vpath
})

# Visualize branch probabilities along trajectory
plot_cells(cds, color_cells_by = "pseudotime", label_branch_points = TRUE)