pacman::p_load(Seurat, tidyverse, qs2, cowplot, ggpubr)
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")
source("codes/activity_helpers.R")
summarise <- dplyr::summarise
mutate <- dplyr::mutate
filter <- dplyr::filterNeutrophil trajectory × pathway / TF / CytoSig
把粒系拟时间与通路、转录因子、细胞因子分数对齐。看发育进度，以及同一位置上两组是否已经不同。组间仍看每个样本一个点。
res_path <- "results/traj_activity"
obj_path <- "data/rdata/activity"
traj_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)
}
n_tf_top <- 30
n_ck_top <- 12
n_bins <- 6
traj_levels <- c(
"GMP",
"Promyelocyte", "MMP8 Pre-Neut", "MMP9 Pre-Neut",
"S100A12 Neut", "SLC8A1 Neut", "ISG15 Neut"
)1. 读 08 拟时间和 09–11 缓存
读入 08 的粒系拟时间和 09–11 的活性分数，按细胞 barcode 对齐。
traj <- qs_read(file.path(traj_path, "seurat_object_neut_traj.qs2"), nthreads = 10)
pt_raw <- pick_sling_pt(traj[[]])
message(
"Using ", pt_raw, " (cells with PT: ",
sum(is.finite(traj[[]][[pt_raw]])), ")"
)
traj_meta <- FetchData(
traj,
vars = intersect(
c("sample", "group", "annot_sub", pt_raw, "CytoTRACE2_Score"),
colnames(traj[[]])
)
) |>
rownames_to_column("barcode")
traj_meta <- traj_meta[is.finite(as.numeric(traj_meta[[pt_raw]])), , drop = FALSE]
traj_meta$group <- factor(as.character(traj_meta$group), levels = c("CTRL", "ACI"))
traj_meta$annot_sub <- factor(as.character(traj_meta$annot_sub), levels = traj_levels)
traj_meta$pt <- scale01(as.numeric(traj_meta[[pt_raw]]))
progeny_wide <- qs_read(file.path(obj_path, "progeny_all.qs2"), nthreads = 10)
tf_wide <- qs_read(file.path(obj_path, "collectri_neut.qs2"), nthreads = 10)
cs_wide <- qs_read(file.path(obj_path, "cytosig_all.qs2"), nthreads = 10)
pw <- features_from_wide(progeny_wide)
tfs <- features_from_wide(tf_wide)
ck <- features_from_wide(cs_wide)
dat_pw <- inner_join(traj_meta, progeny_wide, by = "barcode")
dat_tf <- inner_join(traj_meta, tf_wide, by = "barcode")
dat_cs <- inner_join(traj_meta, cs_wide, by = "barcode")
stopifnot(nrow(dat_pw) > 100, nrow(dat_tf) > 100, nrow(dat_cs) > 100)
print(table(traj_meta$annot_sub, traj_meta$group))
CTRL ACI
GMP 204 1320
Promyelocyte 154 1774
MMP8 Pre-Neut 294 5696
MMP9 Pre-Neut 768 4821
S100A12 Neut 486 7570
SLC8A1 Neut 280 8465
ISG15 Neut 119 748print(c(progeny = nrow(dat_pw), tf = nrow(dat_tf), cytosig = nrow(dat_cs)))progeny tf cytosig
32699 32699 32699 2. 拟时间进度（有没有提前或滞后）
密度图看细胞在轴上堆在哪。箱线是每例的拟时间中位数，这才是组间。ACI 更靠左 = 更多早阶段细胞（或整条轴前移），不一定是同一个细胞变幼稚。
p_den <- plot_pt_density(traj_meta, "pt", title = "Slingshot PT")
ggsave(file.path(res_path, "01_pt_density.pdf"), p_den, width = 7, height = 4, dpi = 600)
p_den
cyto_ok <- traj_meta[is.finite(as.numeric(traj_meta$CytoTRACE2_Score)), , drop = FALSE]
p_cyto <- plot_pt_density(
cyto_ok,
"CytoTRACE2_Score",
title = "CytoTRACE2 (high = earlier)"
)
ggsave(file.path(res_path, "01_cytotrace_density.pdf"), p_cyto, width = 7, height = 4, dpi = 600)
p_cyto
p_stage <- ggplot(traj_meta, aes(x = annot_sub, y = pt, fill = annot_sub)) +
geom_boxplot(outlier.size = 0.15, width = 0.7) +
scale_fill_manual(values = annot_sub_col[traj_levels]) +
theme_classic() +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)
) +
xlab(NULL) +
ylab("Pseudotime 0–1")
ggsave(file.path(res_path, "01_pt_by_subtype.pdf"), p_stage, width = 8, height = 5, dpi = 600)
p_stage
sample_pt <- traj_meta |>
dplyr::group_by(sample, group) |>
dplyr::summarise(
n_cell = dplyr::n(),
pt_median = median(pt, na.rm = TRUE),
cyto_median = median(CytoTRACE2_Score, na.rm = TRUE),
frac_early = mean(pt <= 1 / 3, na.rm = TRUE),
.groups = "drop"
) |>
dplyr::filter(n_cell >= min_cells_sample)
write_csv(sample_pt, file.path(res_path, "01_sample_pseudotime.csv"))
pt_test <- dplyr::bind_rows(
wilcox_ctrl_aci(sample_pt$pt_median, sample_pt$group) |>
dplyr::mutate(feature = "pt_median"),
wilcox_ctrl_aci(sample_pt$cyto_median, sample_pt$group) |>
dplyr::mutate(feature = "cyto_median"),
wilcox_ctrl_aci(sample_pt$frac_early, sample_pt$group) |>
dplyr::mutate(feature = "frac_early")
) |>
add_stars()
write_csv(pt_test, file.path(res_path, "01_sample_pseudotime_wilcox.csv"))
print(pt_test)# A tibble: 3 × 8
n_ctrl n_aci median_ctrl median_aci delta_aci_minus_ctrl p feature star
<int> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 6 5 0.592 0.747 0.156 0.171 pt_median ns
2 6 5 0.0948 0.0843 -0.0105 0.235 cyto_med… ns
3 6 5 0.224 0.119 -0.105 0.171 frac_ear… ns p_box_pt <- ggplot(sample_pt, aes(x = group, y = pt_median, fill = group)) +
geom_boxplot(width = 0.55, outlier.shape = NA) +
geom_jitter(width = 0.12, size = 2.2, alpha = 0.9) +
scale_fill_manual(values = group_col) +
theme_classic() +
theme(legend.position = "none") +
xlab(NULL) +
ylab("Sample median PT")
ggsave(file.path(res_path, "01_sample_pt_box.pdf"), p_box_pt, width = 4.5, height = 5, dpi = 600)
p_box_pt
3. PROGENy 沿拟时间
曲线看活性随拟时间的形状（平移还是交叉）。组间看按拟时间分段后的每例中位数。
pw_rho <- pt_assoc(dat_pw, "pt", pw)
write_csv(pw_rho, file.path(res_path, "02_progeny_pt_spearman.csv"))
print(pw_rho)# A tibble: 12 × 5
feature n rho p_spearman padj
<chr> <int> <dbl> <dbl> <dbl>
1 JAK-STAT 32699 0.754 0 0
2 NFkB 32699 0.435 0 0
3 p53 32699 0.399 0 0
4 PI3K 32699 -0.277 0 0
5 TNFa 32699 0.274 0 0
6 WNT 32699 0.272 0 0
7 Trail 32699 0.227 0 0
8 MAPK 32699 0.225 0 0
9 EGFR 32699 -0.146 8.96e-155 1.20e-154
10 VEGF 32699 0.142 1.97e-146 2.21e-146
11 Hypoxia 32699 0.142 2.03e-146 2.21e-146
12 TGFb 32699 0.102 6.86e- 77 6.86e- 77p_pw <- plot_pt_smooth(dat_pw, "pt", pw, ncol = 4, ylab = "PROGENy mlm")
ggsave(file.path(res_path, "02_progeny_pt_smooth.pdf"), p_pw, width = 12, height = 8, dpi = 600)
p_pw
med_pw <- sample_pt_bin_medians(dat_pw, "pt", pw, n_bins = n_bins)
bin_pw <- sample_group_tests(med_pw, pw, strata_col = "pt_bin") |>
add_stars()
write_csv(bin_pw, file.path(res_path, "02_progeny_pt_bin_sample.csv"))
p_pw_d <- ggplot(bin_pw, aes(x = stratum, y = feature, fill = delta_aci_minus_ctrl)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = star), size = 3) +
scale_fill_activity(bin_pw$delta_aci_minus_ctrl) +
theme_classic() +
xlab("PT bin (b1 early → b6 late)") +
ylab(NULL) +
ggtitle("PROGENy sample median: ACI − CTRL")heatmap fill
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.29241 -0.39100 0.24562 -0.01532 0.73294 2.43547
|x| q80 / q90 / q95: 1.32 / 2.31 / 2.93 ggsave(file.path(res_path, "02_progeny_pt_bin_delta.pdf"), p_pw_d, width = 8, height = 6, dpi = 600)
p_pw_d
print(dplyr::filter(bin_pw, .data$p < 0.05))# A tibble: 26 × 9
stratum feature n_sample_ctrl n_sample_aci median_ctrl median_aci
<chr> <chr> <int> <int> <dbl> <dbl>
1 b3 Hypoxia 6 5 2.51 4.17
2 b3 JAK-STAT 6 5 1.69 2.18
3 b3 p53 6 5 -0.0956 -4.32
4 b4 EGFR 6 5 1.33 0.462
5 b4 Hypoxia 6 5 3.13 4.53
6 b4 JAK-STAT 6 5 2.44 3.89
7 b4 MAPK 6 5 -2.44 -1.59
8 b4 NFkB 6 5 -1.70 -0.975
9 b4 TGFb 6 5 0.537 1.54
10 b4 VEGF 6 5 0.475 1.22
# ℹ 16 more rows
# ℹ 3 more variables: delta_aci_minus_ctrl <dbl>, p <dbl>, star <chr>4. 转录因子沿拟时间
沿拟时间画与 10 相同的 30 个发育相关 TF。
tf_csv <- "results/neut_tf/01_tf_stage.csv"
if (file.exists(tf_csv)) {
tf_stage <- read_csv(tf_csv, show_col_types = FALSE)
} else {
tf_stage <- stage_assoc(dat_tf, tfs, "annot_sub")
}
tf_keep <- tf_stage |>
dplyr::filter(.data$padj_kruskal < 0.05) |>
dplyr::slice_max(order_by = abs(.data$rho), n = n_tf_top, with_ties = FALSE)
if (nrow(tf_keep) < n_tf_top) {
tf_keep <- tf_stage |>
dplyr::slice_max(order_by = abs(.data$rho), n = n_tf_top, with_ties = FALSE)
}
tf_top <- intersect(tf_keep$feature, tfs)
print(tf_top) [1] "RORB" "MYC" "KAT7" "CEBPE" "ATF1" "NFKB1" "ZBED1"
[8] "NEUROD1" "IRF3" "ZBTB4" "CLOCK" "NR4A2" "PPARD" "STAT1"
[15] "TFAM" "GLI3" "WT1" "DMTF1" "RARG" "BMAL1" "JUNB"
[22] "PHOX2B" "JUND" "SIN3A" "AIRE" "ELF5" "DACH1" "HOPX"
[29] "MYT1" "MAFG" tf_rho <- pt_assoc(dat_tf, "pt", tf_top)
write_csv(tf_rho, file.path(res_path, "03_tf_pt_spearman.csv"))
p_tf <- plot_pt_smooth(dat_tf, "pt", tf_top, ncol = 5, ylab = "CollecTRI ulm")
ggsave(file.path(res_path, "03_tf_pt_smooth.pdf"), p_tf, width = 16, height = 14, dpi = 600)
p_tf
med_tf <- sample_pt_bin_medians(dat_tf, "pt", tf_top, n_bins = n_bins)
bin_tf <- sample_group_tests(med_tf, tf_top, strata_col = "pt_bin") |>
add_stars()
write_csv(bin_tf, file.path(res_path, "03_tf_pt_bin_sample.csv"))
p_tf_d <- ggplot(bin_tf, aes(x = stratum, y = feature, fill = delta_aci_minus_ctrl)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = star), size = 2.5) +
scale_fill_activity(bin_tf$delta_aci_minus_ctrl) +
theme_classic() +
theme(axis.text.y = element_text(size = 8)) +
xlab("PT bin (b1 early → b6 late)") +
ylab(NULL) +
ggtitle("TF sample median: ACI − CTRL")heatmap fill
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.61221 0.01501 0.81826 0.68392 1.33141 5.66800
|x| q80 / q90 / q95: 1.76 / 2.29 / 2.81 ggsave(file.path(res_path, "03_tf_pt_bin_delta.pdf"), p_tf_d, width = 8, height = 9, dpi = 600)
p_tf_d
print(dplyr::filter(bin_tf, .data$p < 0.05))# A tibble: 76 × 9
stratum feature n_sample_ctrl n_sample_aci median_ctrl median_aci
<chr> <chr> <int> <int> <dbl> <dbl>
1 b1 RORB 6 5 -0.270 1.52
2 b1 CEBPE 6 5 6.86 12.5
3 b1 ATF1 6 5 2.60 3.97
4 b1 NFKB1 6 5 2.79 3.90
5 b1 NR4A2 6 5 1.08 1.53
6 b1 PPARD 6 5 1.36 1.97
7 b1 STAT1 6 5 1.48 3.13
8 b1 GLI3 6 5 0.0225 1.02
9 b1 DMTF1 6 5 -0.340 1.32
10 b1 RARG 6 5 -0.565 -0.0811
# ℹ 66 more rows
# ℹ 3 more variables: delta_aci_minus_ctrl <dbl>, p <dbl>, star <chr>5. CytoSig 沿拟时间
选出随拟时间变化最大的细胞因子，看差出现在轴的哪一段。总体组间差以 11 为准。
ck_rho <- pt_assoc(dat_cs, "pt", ck)
write_csv(ck_rho, file.path(res_path, "04_cytosig_pt_spearman.csv"))
ck_top <- ck_rho |>
dplyr::slice_max(order_by = abs(.data$rho), n = n_ck_top, with_ties = FALSE) |>
dplyr::pull(feature)
print(ck_rho |> dplyr::filter(.data$feature %in% ck_top))# A tibble: 12 × 5
feature n rho p_spearman padj
<chr> <int> <dbl> <dbl> <dbl>
1 IL10 32699 0.819 0 0
2 IL1B 32699 0.790 0 0
3 TGFB3 32699 0.742 0 0
4 TNFA 32699 0.719 0 0
5 IFNG 32699 0.706 0 0
6 IL4 32699 -0.681 0 0
7 GCSF 32699 0.678 0 0
8 NO 32699 0.665 0 0
9 LIF 32699 0.592 0 0
10 BMP2 32699 0.591 0 0
11 IL1A 32699 0.589 0 0
12 IL17A 32699 0.588 0 0p_cs <- plot_pt_smooth(dat_cs, "pt", ck_top, ncol = 4, ylab = "CytoSig z")
ggsave(file.path(res_path, "04_cytosig_pt_smooth.pdf"), p_cs, width = 12, height = 8, dpi = 600)
p_cs
med_cs <- sample_pt_bin_medians(dat_cs, "pt", ck_top, n_bins = n_bins)
bin_cs <- sample_group_tests(med_cs, ck_top, strata_col = "pt_bin") |>
add_stars()
write_csv(bin_cs, file.path(res_path, "04_cytosig_pt_bin_sample.csv"))
p_cs_d <- ggplot(bin_cs, aes(x = stratum, y = feature, fill = delta_aci_minus_ctrl)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = star), size = 3) +
scale_fill_activity(bin_cs$delta_aci_minus_ctrl) +
theme_classic() +
xlab("PT bin (b1 early → b6 late)") +
ylab(NULL) +
ggtitle("CytoSig sample median: ACI − CTRL")heatmap fill
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.9975 0.2229 0.6400 0.7506 1.2080 2.9352
|x| q80 / q90 / q95: 1.42 / 1.86 / 2.4 ggsave(file.path(res_path, "04_cytosig_pt_bin_delta.pdf"), p_cs_d, width = 8, height = 6, dpi = 600)
p_cs_d
print(dplyr::filter(bin_cs, .data$p < 0.05))# A tibble: 35 × 9
stratum feature n_sample_ctrl n_sample_aci median_ctrl median_aci
<chr> <chr> <int> <int> <dbl> <dbl>
1 b1 IL17A 6 5 0.976 1.68
2 b3 IL10 6 5 0.333 1.11
3 b3 GCSF 6 5 2.04 3.61
4 b3 NO 6 5 1.66 2.78
5 b3 BMP2 6 5 0.0824 0.722
6 b3 IL17A 6 5 1.62 2.41
7 b4 IL10 6 5 0.644 2.36
8 b4 IL1B 6 5 -0.304 1.13
9 b4 TGFB3 6 5 -0.117 1.42
10 b4 TNFA 6 5 -0.0499 0.582
# ℹ 25 more rows
# ℹ 3 more variables: delta_aci_minus_ctrl <dbl>, p <dbl>, star <chr>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] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggpubr_0.6.0 cowplot_1.1.3 qs2_0.1.7 lubridate_1.9.3
[5] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
[9] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
[13] tidyverse_2.0.0 Seurat_5.2.1 SeuratObject_5.0.2 sp_2.1-4
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 jsonlite_1.8.9 magrittr_2.0.3
[4] spatstat.utils_3.1-1 farver_2.1.2 rmarkdown_2.29
[7] ragg_1.3.3 vctrs_0.6.5 ROCR_1.0-11
[10] spatstat.explore_3.3-3 rstatix_0.7.2 htmltools_0.5.8.1
[13] broom_1.0.7 Formula_1.2-5 sctransform_0.4.1
[16] parallelly_1.39.0 KernSmooth_2.23-26 htmlwidgets_1.6.4
[19] ica_1.0-3 plyr_1.8.9 plotly_4.10.4
[22] zoo_1.8-12 igraph_2.1.2 mime_0.12
[25] lifecycle_1.0.4 pkgconfig_2.0.3 Matrix_1.7-3
[28] R6_2.5.1 fastmap_1.2.0 fitdistrplus_1.2-1
[31] future_1.34.0 shiny_1.9.1 digest_0.6.37
[34] colorspace_2.1-1 patchwork_1.3.0 tensor_1.5
[37] RSpectra_0.16-2 irlba_2.3.5.1 textshaping_0.4.0
[40] labeling_0.4.3 progressr_0.15.1 fansi_1.0.6
[43] spatstat.sparse_3.1-0 timechange_0.3.0 mgcv_1.9-3
[46] httr_1.4.7 polyclip_1.10-7 abind_1.4-8
[49] compiler_4.4.3 bit64_4.5.2 withr_3.0.2
[52] backports_1.5.0 carData_3.0-5 fastDummies_1.7.4
[55] ggsignif_0.6.4 MASS_7.3-65 tools_4.4.3
[58] lmtest_0.9-40 httpuv_1.6.15 future.apply_1.11.3
[61] goftest_1.2-3 glue_1.8.0 nlme_3.1-166
[64] promises_1.3.2 grid_4.4.3 Rtsne_0.17
[67] cluster_2.1.8.1 reshape2_1.4.4 generics_0.1.3
[70] gtable_0.3.6 spatstat.data_3.1-4 tzdb_0.4.0
[73] data.table_1.16.2 hms_1.1.3 car_3.1-3
[76] stringfish_0.18.0 utf8_1.2.4 spatstat.geom_3.3-4
[79] RcppAnnoy_0.0.22 ggrepel_0.9.6 RANN_2.6.2
[82] pillar_1.9.0 vroom_1.6.5 spam_2.11-0
[85] RcppHNSW_0.6.0 later_1.4.1 splines_4.4.3
[88] lattice_0.22-6 bit_4.5.0 survival_3.8-3
[91] deldir_2.0-4 tidyselect_1.2.1 miniUI_0.1.1.1
[94] pbapply_1.7-2 knitr_1.49 gridExtra_2.3
[97] scattermore_1.2 xfun_0.49 matrixStats_1.4.1
[100] stringi_1.8.4 lazyeval_0.2.2 yaml_2.3.10
[103] pacman_0.5.1 evaluate_1.0.1 codetools_0.2-20
[106] cli_3.6.3 uwot_0.2.2 RcppParallel_5.1.9
[109] systemfonts_1.1.0 xtable_1.8-4 reticulate_1.42.0
[112] munsell_0.5.1 Rcpp_1.0.13-1 globals_0.16.3
[115] spatstat.random_3.3-2 png_0.1-8 spatstat.univar_3.1-1
[118] parallel_4.4.3 dotCall64_1.2 listenv_0.9.1
[121] viridisLite_0.4.2 scales_1.3.0 ggridges_0.5.6
[124] crayon_1.5.3 rlang_1.1.4