pacman::p_load(
Seurat, tidyverse, qs2, decoupleR, cowplot, tictoc
)
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 lineage TF activity
中性粒谱系（GMP 到成熟粒）的转录因子活性。展示与发育阶段相关最强的 30 个 TF。活性不等于该 TF 自己的 mRNA。
res_path <- "results/neut_tf"
obj_path <- "data/rdata/activity"
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)
n_tf_top <- 30
traj_levels <- c(
"GMP",
"Promyelocyte", "MMP8 Pre-Neut", "MMP9 Pre-Neut",
"S100A12 Neut", "SLC8A1 Neut", "ISG15 Neut"
)
traj_col <- annot_sub_col[traj_levels]seurat_object <- qs_read("data/rdata/seurat_object_final.qs2", nthreads = 10)
stopifnot(all(c("sample", "group", "annot_sub") %in% colnames(seurat_object[[]])))
keep <- as.character(seurat_object$annot_sub) %in% traj_levels
seurat_object <- seurat_object[, keep]
seurat_object$group <- factor(seurat_object$group, levels = c("CTRL", "ACI"))
seurat_object$annot_sub <- factor(
as.character(seurat_object$annot_sub),
levels = traj_levels
)
if (length(Layers(seurat_object, search = "counts")) > 1) {
seurat_object <- JoinLayers(seurat_object)
}
if (!"data" %in% Layers(seurat_object, search = "data")) {
seurat_object <- NormalizeData(seurat_object, verbose = TRUE)
}
cell_meta <- FetchData(
seurat_object,
vars = c("sample", "group", "annot_sub")
) |>
rownames_to_column("barcode")
print(table(cell_meta$annot_sub, cell_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 7481. CollecTRI ulm
每个细胞、每个转录因子一个活性分数。正分表示靶基因整体上调。
tf_qs2 <- file.path(obj_path, "collectri_neut.qs2")
if (cache_ok(tf_qs2, colnames(seurat_object))) {
tf_wide <- qs_read(tf_qs2, nthreads = 10)
} else {
net_tf <- decoupleR::get_collectri(organism = "human", split_complexes = FALSE)
mat <- as.matrix(GetAssayData(seurat_object, assay = "RNA", layer = "data"))
tic("neut CollecTRI run_ulm")
acts_tf <- decoupleR::run_ulm(
mat = mat,
net = net_tf,
.source = "source",
.target = "target",
.mor = "mor",
minsize = 5
)
toc()
tf_wide <- acts_long_to_wide(acts_tf)
qs_save(tf_wide, tf_qs2, nthreads = 10)
}
tfs <- features_from_wide(tf_wide)
dat_tf <- bind_meta_acts(cell_meta, tf_wide)
length(tfs)[1] 7712. 发育阶段 → top TF
先找随 GMP→成熟粒变化最大的 TF，避免一千个 TF 一起报。
只展示和阶段相关最强的 30 个 TF，不是「只有这 30 个重要」。在 Kruskal BH < 0.05 里按 |ρ| 取 30；不够就按 |ρ| 硬取。
tf_stage <- stage_assoc(dat_tf, tfs, "annot_sub")
write_csv(tf_stage, file.path(res_path, "01_tf_stage.csv"))
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 <- tf_keep$feature
print(tf_keep)# A tibble: 30 × 7
feature n rho p_spearman p_kruskal padj_spearman padj_kruskal
<chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 RORB 32699 0.730 0 0 0 0
2 MYC 32699 -0.720 0 0 0 0
3 KAT7 32699 0.698 0 0 0 0
4 CEBPE 32699 -0.693 0 0 0 0
5 ATF1 32699 0.683 0 0 0 0
6 NFKB1 32699 0.671 0 0 0 0
7 ZBED1 32699 -0.670 0 0 0 0
8 NEUROD1 32699 0.669 0 0 0 0
9 IRF3 32699 0.659 0 0 0 0
10 ZBTB4 32699 -0.657 0 0 0 0
# ℹ 20 more rowsmean_tf <- mean_by_scaled(dat_tf, tf_top, "annot_sub")
p_hm_tf <- plot_heatmap_mean(mean_tf, tf_top, row_col = "annot_sub") +
scale_x_discrete(limits = traj_levels)heatmap fill
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.36313 -0.74263 -0.24594 -0.06902 0.56791 2.44061
|x| q80 / q90 / q95: 0.938 / 1.1 / 1.28 ggsave(file.path(res_path, "01_tf_heatmap_stage.pdf"), p_hm_tf, width = 8, height = 9, dpi = 600)
p_hm_tf
long_tf <- dat_tf |>
dplyr::select(annot_sub, dplyr::all_of(tf_top)) |>
tidyr::pivot_longer(dplyr::all_of(tf_top), names_to = "tf", values_to = "score")
long_tf$tf <- factor(long_tf$tf, levels = tf_top)
long_tf$annot_sub <- factor(long_tf$annot_sub, levels = traj_levels)
p_box_tf <- ggplot(long_tf, aes(x = annot_sub, y = score, fill = annot_sub)) +
geom_boxplot(outlier.size = 0.12, width = 0.7) +
facet_wrap(~tf, scales = "free_y", ncol = 6) +
scale_fill_manual(values = traj_col) +
theme_classic() +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1, size = 7)
) +
xlab(NULL) +
ylab("CollecTRI ulm")
ggsave(file.path(res_path, "01_tf_box_stage.pdf"), p_box_tf, width = 16, height = 12, dpi = 600)3. CTRL vs ACI（top TF）
这些 TF 在两组之间有没有差。看样本级点图。
组间比较这 30 个 TF。看每例一个点。
tf_cell_all <- add_fdr(cell_group_tests(dat_tf, tf_top))
tf_med_all <- sample_medians(dat_tf, tf_top)
tf_sample_all <- add_stars(sample_group_tests(tf_med_all, tf_top))
write_csv(tf_cell_all, file.path(res_path, "02_tf_wilcox_cell_overall.csv"))
write_csv(tf_sample_all, file.path(res_path, "02_tf_wilcox_sample_overall.csv"))
tf_cell_sub <- add_fdr(cell_group_tests(dat_tf, tf_top, strata_col = "annot_sub")) |>
dplyr::mutate(annot_sub = factor(.data$stratum, levels = traj_levels))
tf_med_sub <- sample_medians(dat_tf, tf_top, extra_group = "annot_sub")
tf_sample_sub <- add_stars(sample_group_tests(tf_med_sub, tf_top, strata_col = "annot_sub")) |>
dplyr::mutate(annot_sub = factor(.data$stratum, levels = traj_levels))
write_csv(tf_cell_sub, file.path(res_path, "02_tf_wilcox_cell_subtype.csv"))
write_csv(tf_sample_sub, file.path(res_path, "02_tf_wilcox_sample_subtype.csv"))
print(tf_sample_all)# A tibble: 30 × 8
feature n_sample_ctrl n_sample_aci median_ctrl median_aci
<chr> <int> <int> <dbl> <dbl>
1 RORB 6 5 -0.202 4.41
2 MYC 6 5 9.42 8.20
3 KAT7 6 5 4.15 5.48
4 CEBPE 6 5 9.84 9.82
5 ATF1 6 5 4.01 8.35
6 NFKB1 6 5 5.16 7.89
7 ZBED1 6 5 3.67 3.10
8 NEUROD1 6 5 0.831 3.19
9 IRF3 6 5 0.482 2.31
10 ZBTB4 6 5 5.90 3.08
# ℹ 20 more rows
# ℹ 3 more variables: delta_aci_minus_ctrl <dbl>, p <dbl>, star <chr>print(dplyr::filter(tf_sample_sub, .data$p < 0.05))# A tibble: 106 × 10
stratum feature n_sample_ctrl n_sample_aci median_ctrl median_aci
<chr> <chr> <int> <int> <dbl> <dbl>
1 GMP CEBPE 6 5 7.14 11.7
2 GMP CLOCK 6 5 0.403 1.68
3 GMP STAT1 6 5 1.19 2.53
4 GMP WT1 6 5 1.09 1.95
5 GMP DMTF1 6 5 -0.320 1.43
6 GMP RARG 6 5 -0.559 0.0863
7 GMP BMAL1 6 5 0.356 1.54
8 GMP PHOX2B 6 5 -0.269 1.35
9 GMP JUND 6 5 0.425 1.39
10 GMP DACH1 6 5 0 -2.37
# ℹ 96 more rows
# ℹ 4 more variables: delta_aci_minus_ctrl <dbl>, p <dbl>, star <chr>,
# annot_sub <fct>p_tf_d <- ggplot(
dplyr::mutate(tf_sample_sub, feature = factor(feature, levels = tf_top)),
aes(x = annot_sub, y = feature, fill = delta_aci_minus_ctrl)
) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = star), size = 2.6) +
scale_fill_activity(tf_sample_sub$delta_aci_minus_ctrl) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
xlab(NULL) +
ylab(NULL) +
ggtitle("Top TF sample-level ACI − CTRL (Wilcoxon p)")heatmap fill
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.03312 0.02424 0.71868 0.65395 1.33642 7.46858
|x| q80 / q90 / q95: 1.75 / 2.23 / 2.68 ggsave(file.path(res_path, "02_tf_heatmap_sample_delta.pdf"), p_tf_d, width = 8, height = 9, dpi = 600)
p_tf_d
mean_tf_g <- mean_by_scaled(dat_tf, tf_top, c("annot_sub", "group"))
long_g <- tidyr::pivot_longer(
mean_tf_g, dplyr::all_of(tf_top), names_to = "feature", values_to = "mean"
) |>
dplyr::mutate(
feature = factor(feature, levels = tf_top),
annot_sub = factor(annot_sub, levels = traj_levels),
group = factor(group, levels = c("CTRL", "ACI"))
)
p_tf_g <- ggplot(long_g, aes(x = annot_sub, y = feature, fill = mean)) +
geom_tile(color = "white", linewidth = 0.4) +
facet_wrap(~group, ncol = 2) +
scale_fill_activity(long_g$mean) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 9)) +
xlab(NULL) +
ylab(NULL)heatmap fill
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.7899 -0.8488 -0.2830 -0.1872 0.4735 3.2891
|x| q80 / q90 / q95: 1.08 / 1.32 / 1.48 ggsave(file.path(res_path, "02_tf_heatmap_stage_group.pdf"), p_tf_g, width = 12, height = 9, dpi = 600)
p_tf_g
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] tictoc_1.2.1 cowplot_1.1.3 decoupleR_2.12.0 qs2_0.1.7
[5] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[9] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[13] ggplot2_3.5.1 tidyverse_2.0.0 Seurat_5.2.1 SeuratObject_5.0.2
[17] 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 htmltools_0.5.8.1 sctransform_0.4.1
[13] parallelly_1.39.0 KernSmooth_2.23-26 htmlwidgets_1.6.4
[16] ica_1.0-3 plyr_1.8.9 plotly_4.10.4
[19] zoo_1.8-12 igraph_2.1.2 mime_0.12
[22] lifecycle_1.0.4 pkgconfig_2.0.3 Matrix_1.7-3
[25] R6_2.5.1 fastmap_1.2.0 fitdistrplus_1.2-1
[28] future_1.34.0 shiny_1.9.1 digest_0.6.37
[31] colorspace_2.1-1 patchwork_1.3.0 tensor_1.5
[34] RSpectra_0.16-2 irlba_2.3.5.1 textshaping_0.4.0
[37] labeling_0.4.3 progressr_0.15.1 fansi_1.0.6
[40] spatstat.sparse_3.1-0 timechange_0.3.0 httr_1.4.7
[43] polyclip_1.10-7 abind_1.4-8 compiler_4.4.3
[46] bit64_4.5.2 withr_3.0.2 BiocParallel_1.40.2
[49] fastDummies_1.7.4 MASS_7.3-65 tools_4.4.3
[52] lmtest_0.9-40 httpuv_1.6.15 future.apply_1.11.3
[55] goftest_1.2-3 glue_1.8.0 nlme_3.1-166
[58] promises_1.3.2 grid_4.4.3 Rtsne_0.17
[61] cluster_2.1.8.1 reshape2_1.4.4 generics_0.1.3
[64] gtable_0.3.6 spatstat.data_3.1-4 tzdb_0.4.0
[67] data.table_1.16.2 hms_1.1.3 stringfish_0.18.0
[70] utf8_1.2.4 spatstat.geom_3.3-4 RcppAnnoy_0.0.22
[73] ggrepel_0.9.6 RANN_2.6.2 pillar_1.9.0
[76] vroom_1.6.5 spam_2.11-0 RcppHNSW_0.6.0
[79] later_1.4.1 splines_4.4.3 lattice_0.22-6
[82] bit_4.5.0 survival_3.8-3 deldir_2.0-4
[85] tidyselect_1.2.1 miniUI_0.1.1.1 pbapply_1.7-2
[88] knitr_1.49 gridExtra_2.3 scattermore_1.2
[91] xfun_0.49 matrixStats_1.4.1 stringi_1.8.4
[94] lazyeval_0.2.2 yaml_2.3.10 pacman_0.5.1
[97] evaluate_1.0.1 codetools_0.2-20 cli_3.6.3
[100] uwot_0.2.2 RcppParallel_5.1.9 systemfonts_1.1.0
[103] xtable_1.8-4 reticulate_1.42.0 munsell_0.5.1
[106] Rcpp_1.0.13-1 globals_0.16.3 spatstat.random_3.3-2
[109] png_0.1-8 spatstat.univar_3.1-1 parallel_4.4.3
[112] dotCall64_1.2 listenv_0.9.1 viridisLite_0.4.2
[115] scales_1.3.0 ggridges_0.5.6 crayon_1.5.3
[118] rlang_1.1.4