Neutrophil trajectory (GMP + Pre-Neut + Neut)

粒系发育：GMP → 早幼粒 → Pre-Neut → 成熟粒。CytoTRACE2 高分更早。Slingshot 拟时间是转录组距离，不是小时。

骨髓中性粒发育：GMP → Promyelocyte → Pre-Neut → 成熟 Neut。GMP 是粒单祖细胞，也是这条路上 ACI 上调的前体，所以并入轨迹，不单用 Pre-Neut / Neut。细胞集合变了，必须重新 Normalize / HVG / PCA / Harmony / UMAP。不重新聚类，亚型用 06 的 annot_sub。

本脚本只做：CytoTRACE2 潜能、diffusion map、PHATE、slingshot。对照写法见 NMOSD trajectory.R。group 以 CTRL 为对照。

pacman::p_load(
  Seurat, tidyverse, qs2, harmony, tictoc, cowplot, ggh4x, ggpubr,
  CytoTRACE2, destiny, phateR, slingshot, SingleCellExperiment
)
knitr::opts_knit$set(root.dir = "/home/han/Desktop/guoqiang/stroke_bm/stroke_bm_Rproj")
setwd("/home/han/Desktop/guoqiang/stroke_bm/stroke_bm_Rproj")
source("codes/plot_colors.R")
# CytoTRACE2 / destiny 会挂上 plyr，Rfast 会盖掉 Order。
summarise <- dplyr::summarise
mutate <- dplyr::mutate
filter <- dplyr::filter
rename <- dplyr::rename
arrange <- dplyr::arrange
group_by <- dplyr::group_by
ungroup <- dplyr::ungroup
group_modify <- dplyr::group_modify

res_path <- "results/trajectory"
obj_path <- "data/rdata/trajectory"
if (dir.exists(res_path)) {
  old_res <- list.files(res_path, full.names = TRUE)
  if (length(old_res) > 0) unlink(old_res, recursive = TRUE)
} else {
  dir.create(res_path, recursive = TRUE)
}
if (!dir.exists(obj_path)) dir.create(obj_path, recursive = TRUE)

alpha <- 0.35
pt_size <- 0.1
raster <- FALSE
raster_dpi <- c(3000, 3000)
n_hvg <- 3000
min_cells_sample <- 5

traj_levels <- c(
  "GMP",
  "Promyelocyte", "MMP8 Pre-Neut", "MMP9 Pre-Neut",
  "S100A12 Neut", "SLC8A1 Neut", "ISG15 Neut"
)
traj_col <- annot_sub_col[traj_levels]

plot_embed_short <- function(
    obj,
    reduction,
    group.by,
    cols,
    xlab,
    ylab,
    label = TRUE,
    label.size = 4,
    legend = FALSE,
    title = NULL
) {
  axis <- umap_trunc_axis()
  dim_args <- list(
    object = obj,
    reduction = reduction,
    group.by = group.by,
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi,
    label = label,
    label.size = label.size
  )
  if (!is.null(cols)) dim_args$cols <- cols
  p <- do.call(Seurat::DimPlot, dim_args) +
    theme_umap_short() +
    ggplot2::labs(x = xlab, y = ylab) +
    ggplot2::guides(x = axis, y = axis)
  if (isFALSE(legend)) {
    p <- p + ggplot2::guides(color = "none")
  }
  if (!is.null(title)) {
    p <- p +
      ggplot2::ggtitle(title) +
      ggplot2::theme(plot.title = ggplot2::element_text(size = 24))
  }
  p
}

wilcox_ctrl_aci <- function(score, group) {
  group <- factor(as.character(group), levels = c("CTRL", "ACI"))
  ctrl <- score[group == "CTRL"]
  aci <- score[group == "ACI"]
  tibble(
    n_ctrl = length(ctrl),
    n_aci = length(aci),
    median_ctrl = median(ctrl, na.rm = TRUE),
    median_aci = median(aci, na.rm = TRUE),
    delta_aci_minus_ctrl = median(aci, na.rm = TRUE) - median(ctrl, na.rm = TRUE),
    p = if (length(ctrl) >= 2 && length(aci) >= 2) {
      wilcox.test(aci, ctrl, exact = FALSE)$p.value
    } else {
      NA_real_
    }
  )
}

1. Subset GMP + neutrophil lineage

只留 GMP 和中性粒这条发育路上的细胞。

从 seurat_object_final.qs2 取 GMP 和全部 Pre-Neut / Neut 亚型。

test_genes <- qs_read("data/rdata/test_genes.qs2", nthreads = 10)
seurat_object <- qs_read("data/rdata/seurat_object_final.qs2", nthreads = 10)
stopifnot(all(c("sample", "group", "annot_main", "annot_sub") %in%
  colnames(seurat_object[[]])))
seurat_object$group <- factor(seurat_object$group, levels = c("CTRL", "ACI"))

keep <- as.character(seurat_object$annot_sub) %in% traj_levels
print(table(seurat_object$annot_sub[keep], seurat_object$group[keep]))

               
                CTRL  ACI
  HSC              0    0
  MEP              0    0
  GMP            204 1320
  MDP              0    0
  Pro-B            0    0
  Large Pre-B      0    0
  Small Pre-B      0    0
  Trans B          0    0
  Naive B          0    0
  Memory B         0    0
  Plasma           0    0
  CD4 Tnaive       0    0
  CD4 Tcm          0    0
  CD4 Treg         0    0
  CD8 Tnaive       0    0
  CD8 Tem          0    0
  CD8 CTL          0    0
  CD8 NKT          0    0
  MAIT             0    0
  gdT              0    0
  CD56bright NK    0    0
  CD56dim NK       0    0
  CD14 Mono        0    0
  CD16 Mono        0    0
  cDC2             0    0
  pDC              0    0
  Promyelocyte   154 1774
  MMP8 Pre-Neut  294 5696
  MMP9 Pre-Neut  768 4821
  S100A12 Neut   486 7570
  SLC8A1 Neut    280 8465
  ISG15 Neut     119  748
  Early Ery        0    0
  Late Ery         0    0

seurat_object <- seurat_object[, keep]
seurat_object$annot_sub <- factor(
  as.character(seurat_object$annot_sub),
  levels = traj_levels
)
dim(seurat_object)

[1] 38606 32699

if (length(Layers(seurat_object, search = "counts")) > 1) {
  seurat_object <- JoinLayers(seurat_object)
}
seurat_object[["RNA"]]$data <- NULL
seurat_object[["RNA"]]$scale.data <- NULL
DefaultLayer(seurat_object@assays$RNA)

[1] "counts"

Layers(seurat_object)

[1] "counts"

2. Normalize, HVG, Scale, PCA, Harmony

集合变了，重做降维和整合，不重新聚类，亚型沿用 06。

子集上重新高变基因、PCA、Harmony 和 UMAP。亚型沿用 06。

seurat_object <- NormalizeData(seurat_object, verbose = TRUE)
seurat_object <- FindVariableFeatures(
  seurat_object, nfeatures = n_hvg, verbose = TRUE
)
head(VariableFeatures(seurat_object), 20)

 [1] "IGKC"   "PRTN3"  "DEFA1"  "ELANE"  "DEFA3"  "H4C3"   "MPO"    "CTSG"  
 [9] "AZU1"   "DEFA4"  "HBB"    "IGHA2"  "HSPA1B" "HBG2"   "CCL4"   "CCL4L2"
[17] "PRG2"   "CAMP"   "RPL13"  "RPS4X" 

pca_features <- setdiff(
  VariableFeatures(seurat_object),
  test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr]
)
sum(VariableFeatures(seurat_object) %in%
  test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr])

[1] 154

seurat_object <- ScaleData(
  seurat_object, features = pca_features, assay = "RNA"
)
seurat_object <- RunPCA(
  seurat_object, npcs = 50, features = pca_features,
  verbose = TRUE, reduction.name = "pca"
)
ElbowPlot(seurat_object, reduction = "pca", ndims = 50)

tic("Harmony integration")
seurat_object <- RunHarmony(
  seurat_object, c("sample"), ncores = 15, plot_convergence = TRUE
)

toc()

Harmony integration: 17.88 sec elapsed

plot_harmony_umap <- function(obj, dims, file) {
  umap_temp <- RunUMAP(
    obj, dims = 1:dims, reduction = "harmony",
    reduction.name = "umap_harmony"
  )
  p <- plot_grid(
    ncol = 2,
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "annot_sub",
      cols = traj_col, label = TRUE, label.size = 3,
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1),
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "annot_main",
      cols = annot_main_col, label = TRUE,
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1),
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "sample",
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1),
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "group",
      cols = group_col,
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1)
  )
  ggsave(file, plot = p, width = 16, height = 14, dpi = 600)
  p
}

plot_harmony_umap(seurat_object, 20, file.path(res_path, "01_dim_umap20.pdf"))

plot_harmony_umap(seurat_object, 30, file.path(res_path, "01_dim_umap30.pdf"))

plot_harmony_umap(seurat_object, 40, file.path(res_path, "01_dim_umap40.pdf"))

# 看过 20/30/40 后改这里。05 亚型官方 UMAP 用 20。
n_dims <- 20
seurat_object <- RunUMAP(
  seurat_object, dims = 1:n_dims, reduction = "harmony",
  reduction.name = "umap"
)

3. Annotation UMAP

看这条路上各亚型是否按 GMP → 早幼粒 → Pre-Neut → 成熟粒排开。

plot_umap_short(
  seurat_object, group.by = "annot_sub", cols = traj_col,
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi,
  label = TRUE, label.size = 4, legend = TRUE
)

ggsave(file.path(res_path, "02_umap_annot.pdf"), width = 9, height = 8, dpi = 600)
ggsave(file.path(res_path, "02_umap_annot.png"), width = 9, height = 8, dpi = 600)

plot_grid(
  ncol = 2,
  plot_umap_short(
    subset(seurat_object, group == "CTRL"),
    group.by = "annot_sub", cols = traj_col,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi,
    label = FALSE, title = "CTRL"
  ),
  plot_umap_short(
    subset(seurat_object, group == "ACI"),
    group.by = "annot_sub", cols = traj_col,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi,
    label = FALSE, title = "ACI"
  )
)

ggsave(
  file.path(res_path, "02_umap_annot_group.pdf"),
  width = 15, height = 8, dpi = 600
)

plot_grid(
  ncol = 2,
  plot_umap_short(
    seurat_object, group.by = "annot_sub", cols = traj_col,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi,
    label = TRUE, label.size = 3
  ),
  plot_umap_short(
    seurat_object, group.by = "group", cols = group_col,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi,
    label = FALSE, legend = TRUE
  ),
  plot_umap_short(
    seurat_object, group.by = "sample", cols = NULL,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi,
    label = FALSE, legend = TRUE
  ),
  plot_umap_short(
    seurat_object, group.by = "annot_main", cols = annot_main_col,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi,
    label = TRUE, label.size = 4
  )
)

ggsave(
  file.path(res_path, "02_umap_overview.pdf"),
  width = 16, height = 14, dpi = 600
)

4. CytoTRACE2

分数高更「早」、更不成熟。成熟粒应偏低。两组颜色叠在一起，看的是发育早晚，不是疾病标签。

高分更像未成熟细胞。请同时看按亚型上色和按组别上色：若 ACI 只是「成熟粒更多」则平均潜能会下降，那是组成变化，不一定是同一种细胞变幼稚。

5. Diffusion map

另一种摊平方法，用来看是不是拉成一条轴，而不是好几团。

细胞水平 Wilcoxon 会被组成分牵着走（ACI 成熟中性粒更多）。样本水平是每个样本的中位数再比组；亚型内比较更干净。样本至少 min_cells_sample 个细胞才进样本水平。

高分更接近未成熟。组间比较用每例中位数。

cyto_qs2 <- file.path(obj_path, "cytotrace_neut.qs2")
cyto_cache <- NULL
if (file.exists(cyto_qs2)) {
  cyto_cache <- qs_read(cyto_qs2, nthreads = 10)
}
same_n <- !is.null(cyto_cache) && nrow(cyto_cache) == ncol(seurat_object)
same_cells <- same_n &&
  identical(sort(rownames(cyto_cache)), sort(colnames(seurat_object)))
if (same_cells) {
  message(
    "Reuse CytoTRACE2 from qs2 (", nrow(cyto_cache), " cells). Skip cytotrace2()."
  )
  cyto_cache <- cyto_cache[colnames(seurat_object), , drop = FALSE]
  seurat_object <- AddMetaData(seurat_object, cyto_cache)
} else {
  if (!is.null(cyto_cache)) {
    message(
      "CytoTRACE2 qs2 n=", nrow(cyto_cache),
      " vs current n=", ncol(seurat_object), "; re-run."
    )
  }
  counts_mat <- LayerData(seurat_object, assay = "RNA", layer = "counts")
  seurat_object[["RNA"]]$counts <- as(counts_mat, "dgCMatrix")
  seurat_object <- cytotrace2(
    seurat_object, species = "human", is_seurat = TRUE, ncores = 10
  )
  cyto_cols <- grep("^CytoTRACE2", colnames(seurat_object[[]]), value = TRUE)
  cyto_out <- FetchData(seurat_object, vars = cyto_cols)
  qs_save(cyto_out, cyto_qs2, nthreads = 10)
}
table(seurat_object$CytoTRACE2_Potency)

Differentiated      Unipotent    Oligopotent    Multipotent    Pluripotent 
         26434           3048           3217              0              0 
    Totipotent 
             0 

summary(seurat_object$CytoTRACE2_Score)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.002572 0.067856 0.087618 0.130437 0.128374 0.493659 

FeaturePlot(
  seurat_object, reduction = "umap", features = "CytoTRACE2_Score",
  pt.size = pt_size, raster = raster, raster.dpi = raster_dpi
) + theme(aspect.ratio = 1)

ggsave(file.path(res_path, "03_cytotrace_score_umap.pdf"), width = 8, height = 7, dpi = 600)

FeaturePlot(
  seurat_object, reduction = "umap", features = "CytoTRACE2_Relative",
  pt.size = pt_size, raster = raster, raster.dpi = raster_dpi
) + theme(aspect.ratio = 1)

ggsave(
  file.path(res_path, "03_cytotrace_relative_umap.pdf"),
  width = 8, height = 7, dpi = 600
)

cyto_plots <- tryCatch(
  plotData(
    cytotrace2_result = seurat_object,
    annotation = as.character(seurat_object$annot_sub),
    pc_dims = n_dims,
    is_seurat = TRUE
  ),
  error = function(e) {
    message("CytoTRACE2 plotData skipped: ", conditionMessage(e))
    NULL
  }
)
if (!is.null(cyto_plots)) {
  print(cyto_plots$CytoTRACE2_UMAP)
  ggsave(file.path(res_path, "03_cytotrace_umap_pkg.pdf"), width = 8, height = 6)
  print(cyto_plots$CytoTRACE2_Boxplot_byPheno)
  ggsave(file.path(res_path, "03_cytotrace_box_pkg.pdf"), width = 8, height = 6)
}

cyto_meta <- FetchData(
  seurat_object,
  vars = c(
    "sample", "group", "annot_sub", "annot_main",
    "CytoTRACE2_Score", "CytoTRACE2_Relative", "CytoTRACE2_Potency"
  )
) |>
  rownames_to_column("barcode")
write_csv(cyto_meta, file.path(obj_path, "cytotrace_cell.csv"))

亚型间潜能（不分 group）：

ggplot(cyto_meta, aes(CytoTRACE2_Score, annot_sub, fill = annot_sub)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7, width = 0.65) +
  geom_jitter(shape = 16, alpha = 0.05, size = 0.15, height = 0.2) +
  scale_fill_manual(values = traj_col) +
  scale_y_discrete(limits = rev(traj_levels)) +
  xlab("CytoTRACE2 score") +
  ggtitle("Developmental potential") +
  theme_bw() +
  theme(
    legend.position = "none",
    axis.title.y = element_blank(),
    plot.title = element_text(face = "bold", hjust = 0.5)
  )

ggsave(file.path(res_path, "03_cytotrace_by_subtype.pdf"), width = 4.5, height = 5, dpi = 600)

总体 CTRL vs ACI：细胞水平 + 样本水平。

cyto_overall_cell <- wilcox_ctrl_aci(cyto_meta$CytoTRACE2_Score, cyto_meta$group) |>
  mutate(level = "cell", stratum = "all")
cyto_sample_all <- cyto_meta |>
  group_by(sample, group) |>
  summarise(
    n_cell = dplyr::n(),
    CytoTRACE2_Score = median(CytoTRACE2_Score),
    .groups = "drop"
  ) |>
  filter(.data$n_cell >= min_cells_sample)
cyto_overall_sample <- wilcox_ctrl_aci(
  cyto_sample_all$CytoTRACE2_Score, cyto_sample_all$group
) |>
  mutate(level = "sample", stratum = "all")

p_all_cell <- ggplot(cyto_meta, aes(group, CytoTRACE2_Score, fill = group)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7, width = 0.65) +
  geom_jitter(shape = 16, alpha = 0.04, size = 0.15, width = 0.2) +
  scale_fill_manual(values = group_col) +
  stat_compare_means(method = "wilcox.test", label = "p.format") +
  ylab("CytoTRACE2 score") +
  ggtitle("All lineage cells") +
  theme_bw() +
  theme(legend.position = "none", axis.title.x = element_blank())

p_all_sample <- ggplot(
  cyto_sample_all, aes(group, CytoTRACE2_Score, fill = group)
) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7, width = 0.65) +
  geom_jitter(width = 0.12, size = 2, alpha = 0.9) +
  scale_fill_manual(values = group_col) +
  stat_compare_means(method = "wilcox.test", label = "p.format") +
  ylab("Median CytoTRACE2 / sample") +
  ggtitle("Sample medians") +
  theme_bw() +
  theme(legend.position = "none", axis.title.x = element_blank())

plot_grid(p_all_cell, p_all_sample, ncol = 2)

ggsave(file.path(res_path, "03_cytotrace_group_all.pdf"), width = 7, height = 4, dpi = 600)

每个亚型：细胞水平 + 样本水平。

cyto_sub_cell <- cyto_meta |>
  group_by(annot_sub) |>
  group_modify(function(d, key) wilcox_ctrl_aci(d$CytoTRACE2_Score, d$group)) |>
  ungroup() |>
  mutate(level = "cell", stratum = as.character(annot_sub))

cyto_sample_sub <- cyto_meta |>
  group_by(sample, group, annot_sub) |>
  summarise(
    n_cell = dplyr::n(),
    CytoTRACE2_Score = median(CytoTRACE2_Score),
    .groups = "drop"
  ) |>
  filter(.data$n_cell >= min_cells_sample)

cyto_sub_sample <- cyto_sample_sub |>
  group_by(annot_sub) |>
  group_modify(function(d, key) wilcox_ctrl_aci(d$CytoTRACE2_Score, d$group)) |>
  ungroup() |>
  mutate(level = "sample", stratum = as.character(annot_sub))

cyto_stats <- bind_rows(
  cyto_overall_cell, cyto_overall_sample, cyto_sub_cell, cyto_sub_sample
) |>
  mutate(p_adj = p.adjust(p, method = "BH"))
write_csv(cyto_stats, file.path(res_path, "03_cytotrace_wilcox.csv"))
cyto_stats

# A tibble: 16 × 10
   n_ctrl n_aci median_ctrl median_aci delta_aci_minus_ctrl        p level 
    <int> <int>       <dbl>      <dbl>                <dbl>    <dbl> <chr> 
 1   2305 30394      0.0905     0.0874             -0.00312 8.38e- 1 cell  
 2      6     5      0.0948     0.0843             -0.0105  2.35e- 1 sample
 3    204  1320      0.355      0.429               0.0747  1.07e- 3 cell  
 4    154  1774      0.237      0.410               0.173   1.14e-51 cell  
 5    294  5696      0.111      0.196               0.0847  2.27e-55 cell  
 6    768  4821      0.0784     0.0847              0.00626 9.71e- 5 cell  
 7    486  7570      0.0811     0.0710             -0.0101  1.92e-12 cell  
 8    280  8465      0.0828     0.0790             -0.00387 1.17e- 3 cell  
 9    119   748      0.0818     0.0784             -0.00343 2.65e- 3 cell  
10      6     5      0.370      0.433               0.0625  7.84e- 1 sample
11      6     5      0.239      0.410               0.172   1.37e- 2 sample
12      6     5      0.110      0.160               0.0500  1.37e- 2 sample
13      6     5      0.0772     0.0817              0.00451 4.11e- 1 sample
14      4     5      0.0791     0.0737             -0.00541 2.00e- 2 sample
15      4     5      0.0832     0.0795             -0.00371 1.78e- 1 sample
16      4     5      0.0809     0.0794             -0.00147 3.91e- 1 sample
# ℹ 3 more variables: stratum <chr>, annot_sub <fct>, p_adj <dbl>

p_lab_cell <- cyto_sub_cell |>
  mutate(lab = ifelse(is.na(p), "", sprintf("p=%.3g", p)))
ggplot(cyto_meta, aes(group, CytoTRACE2_Score, fill = group)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7, width = 0.65) +
  geom_jitter(shape = 16, alpha = 0.06, size = 0.15, width = 0.18) +
  geom_text(
    data = p_lab_cell, aes(x = 1.5, y = Inf, label = lab),
    vjust = 1.4, size = 3, inherit.aes = FALSE
  ) +
  scale_fill_manual(values = group_col) +
  facet_wrap(~annot_sub, nrow = 2, scales = "free_y") +
  ylab("CytoTRACE2 score") +
  ggtitle("Cell-level, by subtype") +
  theme_bw() +
  theme(legend.position = "none", axis.title.x = element_blank())

ggsave(
  file.path(res_path, "03_cytotrace_group_subtype_cell.pdf"),
  width = 12, height = 6, dpi = 600
)

p_lab_sample <- cyto_sub_sample |>
  mutate(lab = ifelse(is.na(p), "", sprintf("p=%.3g", p)))
ggplot(cyto_sample_sub, aes(group, CytoTRACE2_Score, fill = group)) +
  geom_boxplot(outlier.shape = NA, alpha = 0.7, width = 0.65) +
  geom_jitter(width = 0.12, size = 1.8, alpha = 0.9) +
  geom_text(
    data = p_lab_sample, aes(x = 1.5, y = Inf, label = lab),
    vjust = 1.4, size = 3, inherit.aes = FALSE
  ) +
  scale_fill_manual(values = group_col) +
  facet_wrap(~annot_sub, nrow = 2, scales = "free_y") +
  ylab("Median CytoTRACE2 / sample") +
  ggtitle("Sample-level, by subtype") +
  theme_bw() +
  theme(legend.position = "none", axis.title.x = element_blank())

ggsave(
  file.path(res_path, "03_cytotrace_group_subtype_sample.pdf"),
  width = 12, height = 6, dpi = 600
)

5. Diffusion map

另一种摊平方法，用来看是不是拉成一条轴，而不是好几团。

Harmony 前 n_dims 维，destiny DiffusionMap。

set.seed(123)
dm <- DiffusionMap(Embeddings(seurat_object, "harmony")[, 1:n_dims])
stopifnot(identical(rownames(dm@eigenvectors), colnames(seurat_object)))
seurat_object[["diffusion_map"]] <- CreateDimReducObject(
  embeddings = dm@eigenvectors,
  key = "diffusion_map_",
  assay = "RNA"
)
qs_save(dm, file.path(obj_path, "dm_neut.qs2"), nthreads = 10)
rm(dm)
gc()

            used   (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells   9558050  510.5   14660271  783.0  14660271  783.0
Vcells 199111007 1519.1  577636070 4407.1 737742890 5628.6

plot_embed_short(
  seurat_object, reduction = "diffusion_map", group.by = "annot_sub",
  cols = traj_col, xlab = "DM-1", ylab = "DM-2", label = TRUE, label.size = 4
)

ggsave(file.path(res_path, "04_dm_annot.pdf"), width = 8, height = 8, dpi = 600)

plot_embed_short(
  seurat_object, reduction = "diffusion_map", group.by = "group",
  cols = group_col, xlab = "DM-1", ylab = "DM-2", label = FALSE, legend = TRUE
)

ggsave(file.path(res_path, "04_dm_group.pdf"), width = 8, height = 8, dpi = 600)

FeaturePlot(
  seurat_object, reduction = "diffusion_map", features = "CytoTRACE2_Score",
  pt.size = pt_size, raster = raster, raster.dpi = raster_dpi
) + theme(aspect.ratio = 1)

ggsave(file.path(res_path, "04_dm_cytotrace.pdf"), width = 8, height = 7, dpi = 600)

6. PHATE

再换一种可视化，交叉核对轴是否稳定。

同样用 Harmony 前 n_dims 维。

set.seed(123)
phate_fit <- phate(
  Embeddings(seurat_object, "harmony")[, 1:n_dims],
  knn = 5, gamma = 1, n.jobs = 15, ndim = 2
)

Calculating PHATE...
  Running PHATE on 32699 observations and 20 variables.
  Calculating graph and diffusion operator...
    Calculating KNN search...
    Calculated KNN search in 7.09 seconds.
    Calculating affinities...
    Calculated affinities in 0.08 seconds.
  Calculated graph and diffusion operator in 7.18 seconds.
  Calculating landmark operator...
    Calculating SVD...
    Calculated SVD in 1.41 seconds.
    Calculating KMeans...
    Calculated KMeans in 1.61 seconds.
  Calculated landmark operator in 3.56 seconds.
  Calculating optimal t...
    Automatically selected t = 21
  Calculated optimal t in 0.70 seconds.
  Calculating diffusion potential...
  Calculated diffusion potential in 0.24 seconds.
  Calculating metric MDS...
  Calculated metric MDS in 2.59 seconds.
Calculated PHATE in 14.28 seconds.

colnames(phate_fit$embedding) <- paste0("phate_harmony_", 1:2)
rownames(phate_fit$embedding) <- colnames(seurat_object)
seurat_object[["phate_harmony"]] <- CreateDimReducObject(
  embeddings = phate_fit$embedding,
  key = "phate_harmony_",
  assay = "RNA",
  misc = phate_fit$params
)
qs_save(phate_fit, file.path(obj_path, "phate_neut.qs2"), nthreads = 10)
rm(phate_fit)
gc()

            used   (Mb) gc trigger   (Mb)  max used   (Mb)
Ncells   9592453  512.3   14660271  783.0  14660271  783.0
Vcells 200082174 1526.6  577636070 4407.1 737742890 5628.6

plot_embed_short(
  seurat_object, reduction = "phate_harmony", group.by = "annot_sub",
  cols = traj_col, xlab = "PHATE-1", ylab = "PHATE-2",
  label = TRUE, label.size = 4
)

ggsave(file.path(res_path, "05_phate_annot.pdf"), width = 8, height = 8, dpi = 600)

plot_embed_short(
  seurat_object, reduction = "phate_harmony", group.by = "group",
  cols = group_col, xlab = "PHATE-1", ylab = "PHATE-2",
  label = FALSE, legend = TRUE
)

ggsave(file.path(res_path, "05_phate_group.pdf"), width = 8, height = 8, dpi = 600)

FeaturePlot(
  seurat_object, reduction = "phate_harmony", features = "CytoTRACE2_Score",
  pt.size = pt_size, raster = raster, raster.dpi = raster_dpi
) + theme(aspect.ratio = 1)

ggsave(file.path(res_path, "05_phate_cytotrace.pdf"), width = 8, height = 7, dpi = 600)

7. Slingshot

在 UMAP 上从 GMP 出发画粒系拟时间。拟时间是转录组距离，不是小时。与通路 / TF / 细胞因子的对齐在 12。

起点为 GMP。拟时间供 12 与活性分数对齐。

sling <- list()
sling$dimred <- Embeddings(seurat_object, "umap")[, 1:2]
sling$clustering <- seurat_object$annot_sub
sling$slingshot <- slingshot(
  sling$dimred,
  sling$clustering,
  start.clus = "GMP",
  stretch = 0
)
sling$curves <- slingCurves(sling$slingshot, as.df = TRUE)
sling$pseudotime <- as.data.frame(slingPseudotime(sling$slingshot))
stopifnot(identical(rownames(sling$dimred), colnames(seurat_object)))
print(slingLineages(sling$slingshot))

$Lineage1
[1] "GMP"           "Promyelocyte"  "MMP8 Pre-Neut" "MMP9 Pre-Neut"
[5] "S100A12 Neut"  "SLC8A1 Neut"   "ISG15 Neut"   

pt_cols <- grep("^Lineage", colnames(sling$pseudotime), value = TRUE)
for (ln in pt_cols) {
  seurat_object[[paste0("sling_", ln)]] <- sling$pseudotime[[ln]]
}

curve_xy <- sling$curves
xy <- intersect(colnames(curve_xy), c("umap_1", "UMAP_1", "umap_2", "UMAP_2"))
stopifnot(length(xy) >= 2)
curve_xy <- curve_xy |>
  dplyr::rename(x = dplyr::all_of(xy[1]), y = dplyr::all_of(xy[2]))
lab_df <- curve_xy |>
  dplyr::group_by(Lineage) |>
  dplyr::slice_max(order_by = .data$Order, n = 1, with_ties = FALSE) |>
  dplyr::ungroup() |>
  dplyr::mutate(lab = paste0("L", Lineage))

plot_df <- as.data.frame(sling$dimred) |>
  mutate(annot_sub = sling$clustering)
x_umap <- colnames(plot_df)[1]
y_umap <- colnames(plot_df)[2]

ggplot(plot_df, aes(.data[[x_umap]], .data[[y_umap]])) +
  geom_point(aes(color = annot_sub), shape = 19, size = 0.1, alpha = 0.3) +
  geom_path(
    data = curve_xy |> dplyr::arrange(.data$Lineage, .data$Order),
    aes(x = x, y = y, group = Lineage),
    arrow = arrow(type = "closed", length = unit(0.4, "cm")),
    linewidth = 1.2
  ) +
  geom_text(
    data = lab_df, aes(x = x, y = y, label = lab),
    size = 5, fontface = "bold", vjust = 1.2
  ) +
  scale_color_manual(values = traj_col) +
  theme_classic() +
  theme_umap_short() +
  theme(legend.position = "none") +
  xlab("UMAP-1") +
  ylab("UMAP-2") +
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  guides(x = umap_trunc_axis(), y = umap_trunc_axis())

ggsave(file.path(res_path, "06_slingshot_umap.pdf"), width = 8, height = 8, dpi = 600)

if (length(pt_cols) > 0) {
  FeaturePlot(
    seurat_object, reduction = "umap", features = paste0("sling_", pt_cols[1]),
    pt.size = pt_size, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
  ggsave(
    file.path(res_path, "06_slingshot_pseudotime.pdf"),
    width = 8, height = 7, dpi = 600
  )
}

qs_save(sling, file.path(obj_path, "slingshot_neut.qs2"), nthreads = 10)
qs_save(
  seurat_object,
  file.path(obj_path, "seurat_object_neut_traj.qs2"),
  nthreads = 10
)

sessionInfo()

R version 4.4.3 (2025-02-28)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.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: Asia/Shanghai
tzcode source: system (glibc)

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

other attached packages:
 [1] slingshot_2.14.0            TrajectoryUtils_1.14.0     
 [3] SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
 [5] Biobase_2.66.0              GenomicRanges_1.58.0       
 [7] GenomeInfoDb_1.42.0         IRanges_2.40.0             
 [9] S4Vectors_0.44.0            BiocGenerics_0.52.0        
[11] MatrixGenerics_1.18.0       matrixStats_1.4.1          
[13] princurve_2.1.6             phateR_1.0.7               
[15] destiny_3.20.0              CytoTRACE2_1.1.0           
[17] RSpectra_0.16-2             Rfast_2.1.0                
[19] RcppParallel_5.1.9          RcppZiggurat_0.1.6         
[21] plyr_1.8.9                  Matrix_1.7-3               
[23] magrittr_2.0.3              HiClimR_2.2.1              
[25] doParallel_1.0.17           iterators_1.0.14           
[27] foreach_1.5.2               data.table_1.16.2          
[29] ggpubr_0.6.0                ggh4x_0.2.8                
[31] cowplot_1.1.3               tictoc_1.2.1               
[33] harmony_1.2.3               Rcpp_1.0.13-1              
[35] qs2_0.1.7                   lubridate_1.9.3            
[37] forcats_1.0.0               stringr_1.5.1              
[39] dplyr_1.1.4                 purrr_1.0.2                
[41] readr_2.1.5                 tidyr_1.3.1                
[43] tibble_3.2.1                ggplot2_3.5.1              
[45] tidyverse_2.0.0             Seurat_5.2.1               
[47] SeuratObject_5.0.2          sp_2.1-4                   

loaded via a namespace (and not attached):
  [1] spatstat.sparse_3.1-0     httr_1.4.7               
  [3] RColorBrewer_1.1-3        tools_4.4.3              
  [5] sctransform_0.4.1         backports_1.5.0          
  [7] utf8_1.2.4                R6_2.5.1                 
  [9] lazyeval_0.2.2            uwot_0.2.2               
 [11] withr_3.0.2               gridExtra_2.3            
 [13] progressr_0.15.1          textshaping_0.4.0        
 [15] cli_3.6.3                 pacman_0.5.1             
 [17] spatstat.explore_3.3-3    fastDummies_1.7.4        
 [19] labeling_0.4.3            robustbase_0.99-4-1      
 [21] spatstat.data_3.1-4       proxy_0.4-27             
 [23] ggridges_0.5.6            pbapply_1.7-2            
 [25] systemfonts_1.1.0         parallelly_1.39.0        
 [27] TTR_0.24.4                generics_0.1.3           
 [29] vroom_1.6.5               ica_1.0-3                
 [31] spatstat.random_3.3-2     car_3.1-3                
 [33] fansi_1.0.6               abind_1.4-8              
 [35] lifecycle_1.0.4           scatterplot3d_0.3-44     
 [37] yaml_2.3.10               carData_3.0-5            
 [39] SparseArray_1.6.0         Rtsne_0.17               
 [41] grid_4.4.3                promises_1.3.2           
 [43] crayon_1.5.3              miniUI_0.1.1.1           
 [45] lattice_0.22-6            pillar_1.9.0             
 [47] knitr_1.49                boot_1.3-31              
 [49] future.apply_1.11.3       codetools_0.2-20         
 [51] glue_1.8.0                spatstat.univar_3.1-1    
 [53] pcaMethods_1.98.0         vcd_1.4-13               
 [55] vctrs_0.6.5               png_0.1-8                
 [57] spam_2.11-0               gtable_0.3.6             
 [59] cachem_1.1.0              xfun_0.49                
 [61] S4Arrays_1.6.0            mime_0.12                
 [63] RcppEigen_0.3.4.0.2       survival_3.8-3           
 [65] ncdf4_1.23                fitdistrplus_1.2-1       
 [67] ROCR_1.0-11               nlme_3.1-166             
 [69] xts_0.14.1                bit64_4.5.2              
 [71] RcppAnnoy_0.0.22          rprojroot_2.0.4          
 [73] irlba_2.3.5.1             KernSmooth_2.23-26       
 [75] colorspace_2.1-1          nnet_7.3-20              
 [77] tidyselect_1.2.1          smoother_1.3             
 [79] bit_4.5.0                 compiler_4.4.3           
 [81] curl_6.0.1                DelayedArray_0.32.0      
 [83] plotly_4.10.4             stringfish_0.18.0        
 [85] scales_1.3.0              DEoptimR_1.1-3-1         
 [87] lmtest_0.9-40             hexbin_1.28.5            
 [89] rappdirs_0.3.3            digest_0.6.37            
 [91] goftest_1.2-3             spatstat.utils_3.1-1     
 [93] rmarkdown_2.29            RhpcBLASctl_0.23-42      
 [95] XVector_0.46.0            htmltools_0.5.8.1        
 [97] pkgconfig_2.0.3           sparseMatrixStats_1.18.0 
 [99] fastmap_1.2.0             rlang_1.1.4              
[101] htmlwidgets_1.6.4         ggthemes_5.1.0           
[103] UCSC.utils_1.2.0          DelayedMatrixStats_1.28.0
[105] shiny_1.9.1               farver_2.1.2             
[107] zoo_1.8-12                jsonlite_1.8.9           
[109] Formula_1.2-5             GenomeInfoDbData_1.2.13  
[111] dotCall64_1.2             patchwork_1.3.0          
[113] munsell_0.5.1             reticulate_1.42.0        
[115] stringi_1.8.4             zlibbioc_1.52.0          
[117] MASS_7.3-65               listenv_0.9.1            
[119] ggrepel_0.9.6             deldir_2.0-4             
[121] splines_4.4.3             tensor_1.5               
[123] hms_1.1.3                 igraph_2.1.2             
[125] ranger_0.17.0             spatstat.geom_3.3-4      
[127] ggsignif_0.6.4            RcppHNSW_0.6.0           
[129] reshape2_1.4.4            evaluate_1.0.1           
[131] laeken_0.5.3              tzdb_0.4.0               
[133] httpuv_1.6.15             VIM_6.2.2                
[135] RANN_2.6.2                polyclip_1.10-7          
[137] knn.covertree_1.0         future_1.34.0            
[139] scattermore_1.2           broom_1.0.7              
[141] xtable_1.8-4              e1071_1.7-16             
[143] rstatix_0.7.2             later_1.4.1              
[145] ragg_1.3.3                viridisLite_0.4.2        
[147] class_7.3-23              memoise_2.0.1            
[149] cluster_2.1.8.1           ggplot.multistats_1.0.1  
[151] timechange_0.3.0          globals_0.16.3           
[153] here_1.0.1