30 Trajectory Analysis (Monocle3)
Monocle3 orders cells along a learned principal graph (pseudotime) rooted
at a chosen starting cluster, visualising differentiation-like structure as
a single branching tree over the WNN embedding. cluster_6 is used as the
root, as in the original analysis; clusters 6 and 7 are excluded before
learning the trajectory, also matching the original analysis. The
trajectory is learned once on the full dataset, and again separately on the
PASC-only and Healthy-only subsets, to compare pseudotime progression and
branch usage between the two conditions on a matched cluster structure.
30.1 Load object and build a 3D companion embedding
A 3D UMAP is computed alongside the existing 2D integration.wnn embedding
from the same weighted nearest-neighbour graph, for the 3D trajectory view
further down.
library(SeuratWrappers)
library(monocle3)
library(ggplot2)
library(dplyr)
library(tidyr)
set.seed(1234)
combined@assays[["regulon"]]@key <- "regulon_"
combined <- RunUMAP(combined, nn.name = "weighted.nn", assay = "RNA",
reduction.name = "integration.wnn.3d", n.components = 3)
# this was with the first run of the object
combined <- combined[, !combined$cell_labels %in% c("cluster_5", "cluster_7", "cluster_8")]30.2 Build trajectory helper
as.cell_data_set() embeds Seurat reductions into the monocle
cell_data_set by copying whatever reducedDims slots the intermediate
SingleCellExperiment coercion produced. The original script relied on
those slots’ numeric position (reducedDims(sce)[[5]],
reducedDims(sce)[[6]], …), which is fragile and breaks silently if
Seurat’s internal reduction ordering changes; the UMAP embeddings needed
for cluster_cells()/learn_graph() are assigned directly by name
instead.
build_trajectory <- function(obj, root_cluster = "cluster_6") {
cds <- as.cell_data_set(obj, assay = "RNA")
umap2d <- Embeddings(obj, "integration.wnn")[colnames(cds), , drop = FALSE]
umap3d <- Embeddings(obj, "integration.wnn.3d")[colnames(cds), , drop = FALSE]
reducedDim(cds, "UMAP") <- umap2d
reducedDim(cds, "UMAP3D") <- umap3d
root_cells <- colnames(obj)[obj$cell_labels == root_cluster]
cds <- cluster_cells(cds, reduction_method = "UMAP")
cds <- learn_graph(cds, use_partition = TRUE)
cds <- order_cells(cds, reduction_method = "UMAP", root_cells = root_cells)
cds
}
build_trajectory_3d <- function(obj, root_cluster = "cluster_6") {
cds <- as.cell_data_set(obj, assay = "RNA")
umap3d <- Embeddings(obj, "integration.wnn.3d")[colnames(cds), , drop = FALSE]
reducedDim(cds, "UMAP") <- umap3d # must be literally "UMAP" for plot_cells_3d
root_cells <- colnames(obj)[obj$cell_labels == root_cluster]
cds <- cluster_cells(cds, reduction_method = "UMAP")
cds <- learn_graph(cds, use_partition = TRUE)
cds <- order_cells(cds, reduction_method = "UMAP", root_cells = root_cells)
cds
}
summarise_trajectory <- function(cds) {
pseudo_df <- data.frame(
celllabels = cds$cell_labels,
condition = cds$individual_condition,
pseudo = pseudotime(cds)
)
print(
pseudo_df %>%
group_by(celllabels, condition) %>%
summarise(mean_pseudo = mean(pseudo, na.rm = TRUE), .groups = "drop")
)
print(
ggplot(pseudo_df, aes(x = celllabels, y = pseudo, fill = condition)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
)
closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
count_node_df <- data.frame(table(data.frame(
celllabels = cds$cell_labels,
condition = cds$individual_condition,
node = as.matrix(closest_vertex[colnames(cds), ])
)))
print(
ggplot(count_node_df, aes(y = Freq, x = node, fill = celllabels)) +
geom_bar(position = "fill", stat = "identity") +
xlab("Trajectory Nodes") + ylab("Percentage of Cells")+
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
)
print(
plot_cells(cds, color_cells_by = "pseudotime", show_trajectory_graph = TRUE,
cell_size = 1, trajectory_graph_segment_size = 0.8, graph_label_size = 5)
)
print(
plot_cells(cds, color_cells_by = "pseudotime", show_trajectory_graph = TRUE,
cell_size = 1, trajectory_graph_segment_size = 0.8, graph_label_size = 5,
label_leaves = FALSE, label_branch_points = FALSE) +
scale_colour_gradient(low = "grey96", high = "grey10")
)
}30.3 Full dataset
cds_all <- build_trajectory(combined)
#> | | | 0% | |==================================================| 100%
summarise_trajectory(cds_all)
#> # A tibble: 12 × 3
#> celllabels condition mean_pseudo
#> <chr> <chr> <dbl>
#> 1 cluster_0 Healthy 4.32
#> 2 cluster_0 PASC 3.42
#> 3 cluster_1 Healthy 7.83
#> 4 cluster_1 PASC 8.14
#> 5 cluster_2 Healthy 1.85
#> 6 cluster_2 PASC 1.92
#> 7 cluster_3 Healthy 2.72
#> 8 cluster_3 PASC 2.86
#> 9 cluster_4 Healthy 3.72
#> 10 cluster_4 PASC 3.86
#> 11 cluster_6 Healthy 0.185
#> 12 cluster_6 PASC 0.135
30.3.1 Top markers by cell type
marker_test_res <- top_markers(cds_all, group_cells_by = "cell_labels",
reference_cells = 1000, cores = 8)
top_specific_marker_ids <- marker_test_res %>%
dplyr::filter(fraction_expressing >= 0.10) %>%
group_by(cell_group) %>%
top_n(1, pseudo_R2) %>%
pull(gene_id) %>%
unique()
rowData(cds_all)$gene_short_name <- row.names(rowData(cds_all))
plot_genes_by_group(cds_all, top_specific_marker_ids, group_cells_by = "cell_labels",
ordering_type = "maximal_on_diag", max.size = 3)
30.3.2 3D trajectory
cds_all_3d <- build_trajectory_3d(combined)
#> | | | 0% | |==================================================| 100%
plot_cells_3d(cds_all_3d, reduction_method = "UMAP", color_cells_by = "pseudotime")30.4 PASC only
cds_pasc <- build_trajectory(combined[, combined$individual_condition == "PASC"])
#> | | | 0% | |==================================================| 100%
summarise_trajectory(cds_pasc)
#> # A tibble: 6 × 3
#> celllabels condition mean_pseudo
#> <chr> <chr> <dbl>
#> 1 cluster_0 PASC 2.94
#> 2 cluster_1 PASC 8.03
#> 3 cluster_2 PASC 0.0325
#> 4 cluster_3 PASC 1.11
#> 5 cluster_4 PASC 2.95
#> 6 cluster_6 PASC 0.0198
30.5 Healthy only
cds_healthy <- build_trajectory(combined[, combined$individual_condition == "Healthy"])
#> | | | 0% | |==================================================| 100%
summarise_trajectory(cds_healthy)
#> # A tibble: 6 × 3
#> celllabels condition mean_pseudo
#> <chr> <chr> <dbl>
#> 1 cluster_0 Healthy 4.41
#> 2 cluster_1 Healthy 4.78
#> 3 cluster_2 Healthy 0.982
#> 4 cluster_3 Healthy 2.41
#> 5 cluster_4 Healthy 2.40
#> 6 cluster_6 Healthy 1.05
30.6 Save results
combined$pseudotime_all <- NA
combined$pseudotime_all[colnames(cds_all)] <- pseudotime(cds_all)
combined$pseudotime_pasc <- NA
combined$pseudotime_pasc[colnames(cds_pasc)] <- pseudotime(cds_pasc)
combined$pseudotime_healthy <- NA
combined$pseudotime_healthy[colnames(cds_healthy)] <- pseudotime(cds_healthy)
SaveSeuratRds(combined, "data/seurat_object_cross_modal_integration_peakslinks_motifs_grn_regulons_monoclepseudotime.rds")
saveRDS(
list(cds_all = cds_all, cds_pasc = cds_pasc, cds_healthy = cds_healthy,
marker_test_res = marker_test_res, top_specific_marker_ids = top_specific_marker_ids),
"data/monocle3_trajectory_results.rds"
)