PROGENy pathway activity

PROGENy 通路活性（JAK-STAT、NFκB、TNFα 等），从靶基因推断，不是磷酸化检测。正分更活跃。组间看每个样本一个点。

全细胞 PROGENy（decoupleR mlm），再对 08 中性粒谱系（GMP → 成熟 Neut）做阶段和分组。跑 mlm 前从网络里去掉 Androgen / Estrogen。mlm 只跑一次。统计、箱线、直方图用 mlm 原始 t 值；类群均值热图按[官方 scRNA 教程](https://saezlab.github.io/decoupleR/articles/pw_sc.html)先对每个通路跨细胞 z-score（ScaleData）再取均值。Wilcoxon Δ 热图仍用原始分数。转录因子见 10_neut_decouple.qmd。细胞级 Wilcoxon 用 BH；样本级用原始 p。group 以 CTRL 为对照。

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::filter

res_path <- "results/progeny"
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)

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_main", "annot_sub") %in%
  colnames(seurat_object[[]])))
seurat_object$group <- factor(seurat_object$group, levels = c("CTRL", "ACI"))
seurat_object$annot_main <- factor(
  as.character(seurat_object$annot_main),
  levels = names(annot_main_col)
)
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_main", "annot_sub")
) |>
  rownames_to_column("barcode")
print(table(cell_meta$annot_main, cell_meta$group))

          
            CTRL   ACI
  HSPC      1390  1976
  B         3641   664
  Plasma     445   109
  T         6894  3560
  NK         954   793
  Mono/DC   2060  1850
  pDC        363   103
  Pre-Neut  1216 12291
  Neut       885 16783
  Ery       4890   323

print(table(cell_meta$annot_sub[cell_meta$annot_sub %in% traj_levels],
  cell_meta$group[cell_meta$annot_sub %in% traj_levels]))

               
                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

1. PROGENy mlm

每个细胞、每条通路一个活性分数。正分更活跃。这是 RNA 推断，不是磷酸化。

PROGENy 取 top 500 靶基因，去掉 Androgen / Estrogen 后做多元回归（mlm）。正分表示该通路靶基因整体上调。

pw_drop <- c("Androgen", "Estrogen")
progeny_qs2 <- file.path(obj_path, "progeny_all.qs2")
reuse <- FALSE
if (cache_ok(progeny_qs2, colnames(seurat_object))) {
  progeny_wide <- qs_read(progeny_qs2, nthreads = 10)
  reuse <- !any(pw_drop %in% colnames(progeny_wide))
}
if (!reuse) {
  net <- decoupleR::get_progeny(organism = "human", top = 500) |>
    dplyr::filter(!.data$source %in% pw_drop)
  stopifnot(!any(pw_drop %in% net$source))
  mat <- as.matrix(GetAssayData(seurat_object, assay = "RNA", layer = "data"))
  tic("PROGENy run_mlm")
  acts <- decoupleR::run_mlm(
    mat = mat,
    net = net,
    .source = "source",
    .target = "target",
    .mor = "weight",
    minsize = 5
  )
  toc()
  progeny_wide <- acts_long_to_wide(acts)
  qs_save(progeny_wide, progeny_qs2, nthreads = 10)
}
pw <- features_from_wide(progeny_wide)
stopifnot(!any(pw_drop %in% pw))
dat <- bind_meta_acts(cell_meta, progeny_wide)
stopifnot(nrow(dat) == nrow(cell_meta))
length(pw)

[1] 12

print(pw)

 [1] "EGFR"     "Hypoxia"  "JAK-STAT" "MAPK"     "NFkB"     "PI3K"    
 [7] "TGFb"     "TNFa"     "Trail"    "VEGF"     "WNT"      "p53"     

2. 大类之间

各条通路在不同大类里相对高低。颜色是标准化后的相对值，用来横着比较通路，不是蛋白量。

热图比较各大类的相对高低：每条通路先在细胞间 z-score，再按大类取均值。

kw_main <- kruskal_by_type(dat, pw, "annot_main")
write_csv(kw_main, file.path(res_path, "01_kruskal_annot_main.csv"))
print(kw_main)

# A tibble: 12 × 5
   feature      n     p  padj star 
   <chr>    <int> <dbl> <dbl> <chr>
 1 EGFR     61190     0     0 ***  
 2 Hypoxia  61190     0     0 ***  
 3 JAK-STAT 61190     0     0 ***  
 4 MAPK     61190     0     0 ***  
 5 NFkB     61190     0     0 ***  
 6 PI3K     61190     0     0 ***  
 7 TGFb     61190     0     0 ***  
 8 TNFa     61190     0     0 ***  
 9 Trail    61190     0     0 ***  
10 VEGF     61190     0     0 ***  
11 WNT      61190     0     0 ***  
12 p53      61190     0     0 ***  

mean_main <- mean_by_scaled(dat, pw, "annot_main")
p_hm <- plot_heatmap_mean(mean_main, pw, row_col = "annot_main") +
  scale_x_discrete(limits = names(annot_main_col))

heatmap fill 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-1.3093 -0.5471 -0.2220 -0.1838  0.1829  1.1293 
|x| q80 / q90 / q95: 0.642 / 0.852 / 1.01 

ggsave(file.path(res_path, "01_heatmap_annot_main.pdf"), p_hm, width = 8, height = 6, dpi = 600)
p_hm

long_pw <- dat |>
  dplyr::select(barcode, annot_main, dplyr::all_of(pw)) |>
  tidyr::pivot_longer(dplyr::all_of(pw), names_to = "pathway", values_to = "score")
long_pw$pathway <- factor(long_pw$pathway, levels = pw)
long_pw$annot_main <- factor(long_pw$annot_main, levels = names(annot_main_col))
p_box <- ggplot(long_pw, aes(x = annot_main, y = score, fill = annot_main)) +
  geom_boxplot(outlier.size = 0.15, width = 0.7) +
  facet_wrap(~pathway, scales = "free_y", ncol = 4) +
  scale_fill_manual(values = annot_main_col) +
  theme_classic() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1, size = 8)
  ) +
  xlab(NULL) +
  ylab("PROGENy mlm")
ggsave(file.path(res_path, "01_box_annot_main.pdf"), p_box, width = 14, height = 10, dpi = 600)
p_box

3. CTRL vs ACI（每个大类）

每个大类里 ACI 减 CTRL。红 = ACI 更高。组间仍叠着平台差异。

直方图看同一亚型里两组分布是否错开；组间结论看每例一个点。细胞很多时分布错开不必等于每例病人都变了。

cell_main <- cell_group_tests(dat, pw, strata_col = "annot_main") |>
  add_fdr() |>
  dplyr::mutate(annot_main = factor(.data$stratum, levels = names(annot_main_col))) |>
  dplyr::arrange(.data$annot_main, .data$padj)
write_csv(cell_main, file.path(res_path, "02_wilcox_cell_annot_main.csv"))

med_main <- sample_medians(dat, pw, extra_group = "annot_main")
sample_main <- sample_group_tests(med_main, pw, strata_col = "annot_main") |>
  add_stars() |>
  dplyr::mutate(annot_main = factor(.data$stratum, levels = names(annot_main_col))) |>
  dplyr::arrange(.data$annot_main, .data$p)
write_csv(sample_main, file.path(res_path, "02_wilcox_sample_annot_main.csv"))

print(dplyr::filter(sample_main, .data$p < 0.05))

# A tibble: 48 × 10
   stratum feature n_sample_ctrl n_sample_aci median_ctrl median_aci
   <chr>   <chr>           <int>        <int>       <dbl>      <dbl>
 1 HSPC    NFkB                6            5     -2.55       -0.926
 2 HSPC    PI3K                6            5     -2.92       -2.04 
 3 HSPC    Trail               6            5     -2.40       -1.44 
 4 HSPC    TGFb                6            5      0.0988      0.551
 5 B       Trail               6            5     -0.736      -1.64 
 6 B       p53                 6            5     -3.80       -6.16 
 7 B       NFkB                6            5     -1.73       -1.22 
 8 B       EGFR                6            5     -1.11       -2.32 
 9 Plasma  MAPK                6            4     -2.31       -1.13 
10 Plasma  NFkB                6            4     -2.20       -0.865
# ℹ 38 more rows
# ℹ 4 more variables: delta_aci_minus_ctrl <dbl>, p <dbl>, star <chr>,
#   annot_main <fct>

hm_delta <- sample_main |>
  dplyr::mutate(
    feature = factor(.data$feature, levels = pw),
    annot_main = factor(.data$annot_main, levels = names(annot_main_col))
  )
p_d <- ggplot(hm_delta, aes(x = annot_main, y = feature, fill = delta_aci_minus_ctrl)) +
  geom_tile(color = "white", linewidth = 0.4) +
  geom_text(aes(label = star), size = 3) +
  scale_fill_activity(hm_delta$delta_aci_minus_ctrl) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  xlab(NULL) +
  ylab(NULL) +
  ggtitle("Sample-level median: ACI − CTRL (star = Wilcoxon p)")

heatmap fill 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-4.04916 -0.44520  0.18574 -0.02489  0.63424  1.62812 
|x| q80 / q90 / q95: 1.08 / 1.39 / 1.94 

ggsave(file.path(res_path, "02_heatmap_sample_delta.pdf"), p_d, width = 8, height = 6, dpi = 600)
p_d

hm_cell <- cell_main |>
  dplyr::mutate(
    feature = factor(.data$feature, levels = pw),
    annot_main = factor(.data$annot_main, levels = names(annot_main_col))
  )
p_c <- ggplot(hm_cell, aes(x = annot_main, y = feature, fill = delta_aci_minus_ctrl)) +
  geom_tile(color = "white", linewidth = 0.4) +
  geom_text(aes(label = star), size = 3) +
  scale_fill_activity(hm_cell$delta_aci_minus_ctrl) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  xlab(NULL) +
  ylab(NULL) +
  ggtitle("Cell-level: ACI − CTRL (star = BH on Wilcoxon)")

heatmap fill 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-3.9964 -0.4273  0.1647 -0.0164  0.6209  2.1853 
|x| q80 / q90 / q95: 1.25 / 1.46 / 2.4 

ggsave(file.path(res_path, "02_heatmap_cell_delta.pdf"), p_c, width = 8, height = 6, dpi = 600)
p_c

样本级箱线：画样本级 p < 0.05 的通路；若没有则画 |Δ| 最大的 4 个。

sig_pw <- unique(as.character(sample_main$feature[sample_main$p < 0.05]))
if (length(sig_pw) == 0) {
  sig_pw <- sample_main |>
    dplyr::group_by(feature) |>
    dplyr::summarise(m = max(abs(delta_aci_minus_ctrl), na.rm = TRUE), .groups = "drop") |>
    dplyr::slice_max(.data$m, n = 4) |>
    dplyr::pull(feature)
}
long_med <- med_main |>
  tidyr::pivot_longer(dplyr::all_of(pw), names_to = "pathway", values_to = "score") |>
  dplyr::filter(.data$pathway %in% sig_pw)
long_med$annot_main <- factor(long_med$annot_main, levels = names(annot_main_col))
long_med$group <- factor(long_med$group, levels = c("CTRL", "ACI"))
p_s <- ggplot(long_med, aes(x = annot_main, y = score, fill = group)) +
  geom_boxplot(outlier.size = 0.6, width = 0.7, position = position_dodge(0.8)) +
  geom_point(
    position = position_dodge(0.8), size = 1.4, alpha = 0.8,
    aes(group = group)
  ) +
  facet_wrap(~pathway, scales = "free_y") +
  scale_fill_manual(values = group_col) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  xlab(NULL) +
  ylab("Sample median PROGENy")
ggsave(file.path(res_path, "02_box_sample_sig.pdf"), p_s, width = 12, height = 8, dpi = 600)
p_s

细胞级 / 样本级直方图：每个大类 × 通路，CTRL vs ACI。

long_cell_main <- dat |>
  dplyr::select(annot_main, group, dplyr::all_of(pw)) |>
  tidyr::pivot_longer(dplyr::all_of(pw), names_to = "pathway", values_to = "score") |>
  dplyr::mutate(
    annot_main = factor(.data$annot_main, levels = names(annot_main_col)),
    pathway = factor(.data$pathway, levels = pw),
    group = factor(.data$group, levels = c("CTRL", "ACI"))
  )
p_hist_cell_main <- plot_score_hist(long_cell_main, "annot_main", bins = 40) +
  ggtitle("Cell-level")
ggsave(
  file.path(res_path, "02_hist_cell_annot_main.pdf"),
  p_hist_cell_main, width = 22, height = 16, dpi = 600
)
p_hist_cell_main

long_med_main <- med_main |>
  tidyr::pivot_longer(dplyr::all_of(pw), names_to = "pathway", values_to = "score") |>
  dplyr::mutate(
    annot_main = factor(.data$annot_main, levels = names(annot_main_col)),
    pathway = factor(.data$pathway, levels = pw),
    group = factor(.data$group, levels = c("CTRL", "ACI"))
  )
p_hist_med_main <- plot_score_hist(long_med_main, "annot_main", bins = 10) +
  ggtitle("Patient-level (sample median)")
ggsave(
  file.path(res_path, "02_hist_sample_annot_main.pdf"),
  p_hist_med_main, width = 22, height = 16, dpi = 600
)
p_hist_med_main

4. 中性粒谱系：发育阶段

从 GMP 到成熟粒，通路是否随发育升降。

从左 GMP 到右成熟粒。颜色沿阶段变化反映发育程序；脑梗相关的组间差在下一节。

dat_neut <- dat |>
  dplyr::filter(as.character(.data$annot_sub) %in% traj_levels) |>
  dplyr::mutate(
    annot_sub = factor(as.character(.data$annot_sub), levels = traj_levels)
  )
stopifnot(nrow(dat_neut) > 0)
pw_stage <- stage_assoc(dat_neut, pw, "annot_sub")
write_csv(pw_stage, file.path(res_path, "03_neut_stage.csv"))
print(pw_stage)

# A tibble: 12 × 7
   feature      n    rho p_spearman p_kruskal padj_spearman padj_kruskal
   <chr>    <int>  <dbl>      <dbl>     <dbl>         <dbl>        <dbl>
 1 JAK-STAT 32699  0.749  0                 0     0                    0
 2 NFkB     32699  0.426  0                 0     0                    0
 3 p53      32699  0.393  0                 0     0                    0
 4 WNT      32699  0.271  0                 0     0                    0
 5 TNFa     32699  0.266  0                 0     0                    0
 6 PI3K     32699 -0.257  0                 0     0                    0
 7 Trail    32699  0.246  0                 0     0                    0
 8 MAPK     32699  0.208  1.06e-315         0     1.59e-315            0
 9 Hypoxia  32699  0.158  5.58e-181         0     7.44e-181            0
10 VEGF     32699  0.150  1.19e-164         0     1.43e-164            0
11 EGFR     32699 -0.139  9.50e-142         0     1.04e-141            0
12 TGFb     32699  0.117  9.64e-101         0     9.64e-101            0

mean_pw <- mean_by_scaled(dat_neut, pw, "annot_sub")
p_hm_pw <- plot_heatmap_mean(mean_pw, pw, row_col = "annot_sub") +
  scale_x_discrete(limits = traj_levels)

heatmap fill 
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-1.433538 -0.464712 -0.007734 -0.079204  0.248222  3.431759 
|x| q80 / q90 / q95: 0.579 / 0.769 / 0.978 

ggsave(file.path(res_path, "03_neut_heatmap_stage.pdf"), p_hm_pw, width = 8, height = 6, dpi = 600)
p_hm_pw

long_n <- dat_neut |>
  dplyr::select(annot_sub, dplyr::all_of(pw)) |>
  tidyr::pivot_longer(dplyr::all_of(pw), names_to = "pathway", values_to = "score")
long_n$annot_sub <- factor(long_n$annot_sub, levels = traj_levels)
p_box_n <- ggplot(long_n, aes(x = annot_sub, y = score, fill = annot_sub)) +
  geom_boxplot(outlier.size = 0.15, width = 0.7) +
  facet_wrap(~pathway, scales = "free_y", ncol = 4) +
  scale_fill_manual(values = traj_col) +
  theme_classic() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1, size = 8)
  ) +
  xlab(NULL) +
  ylab("PROGENy mlm")
ggsave(file.path(res_path, "03_neut_box_stage.pdf"), p_box_n, width = 14, height = 10, dpi = 600)
p_box_n

5. 中性粒谱系：CTRL vs ACI

粒系上 ACI 与 CTRL 的通路差。组间看每例中位数。

谱系合计 + 每个 annot_sub。样本级原始 p。

pw_cell_all <- add_fdr(cell_group_tests(dat_neut, pw))
pw_med_all <- sample_medians(dat_neut, pw)
pw_sample_all <- add_stars(sample_group_tests(pw_med_all, pw))
write_csv(pw_cell_all, file.path(res_path, "04_neut_wilcox_cell_overall.csv"))
write_csv(pw_sample_all, file.path(res_path, "04_neut_wilcox_sample_overall.csv"))

pw_cell_sub <- add_fdr(cell_group_tests(dat_neut, pw, strata_col = "annot_sub")) |>
  dplyr::mutate(annot_sub = factor(.data$stratum, levels = traj_levels))
pw_med_sub <- sample_medians(dat_neut, pw, extra_group = "annot_sub")
pw_sample_sub <- add_stars(sample_group_tests(pw_med_sub, pw, strata_col = "annot_sub")) |>
  dplyr::mutate(annot_sub = factor(.data$stratum, levels = traj_levels))
write_csv(pw_cell_sub, file.path(res_path, "04_neut_wilcox_cell_subtype.csv"))
write_csv(pw_sample_sub, file.path(res_path, "04_neut_wilcox_sample_subtype.csv"))
print(dplyr::filter(pw_sample_all, .data$p < 0.05))

# A tibble: 5 × 8
  feature n_sample_ctrl n_sample_aci median_ctrl median_aci delta_aci_minus_ctrl
  <chr>           <int>        <int>       <dbl>      <dbl>                <dbl>
1 Hypoxia             6            5       2.53      4.28                  1.75 
2 NFkB                6            5      -1.28     -0.0219                1.25 
3 TGFb                6            5       0.146     1.33                  1.19 
4 VEGF                6            5       0.304     1.05                  0.745
5 p53                 6            5      -0.116    -3.31                 -3.20 
# ℹ 2 more variables: p <dbl>, star <chr>

print(dplyr::filter(pw_sample_sub, .data$p < 0.05))

# A tibble: 33 × 10
   stratum       feature  n_sample_ctrl n_sample_aci median_ctrl median_aci
   <chr>         <chr>            <int>        <int>       <dbl>      <dbl>
 1 MMP8 Pre-Neut JAK-STAT             6            5       1.54      2.16  
 2 MMP8 Pre-Neut p53                  6            5      -0.572    -4.32  
 3 MMP9 Pre-Neut EGFR                 6            5       1.55      0.627 
 4 MMP9 Pre-Neut Hypoxia              6            5       2.78      4.39  
 5 MMP9 Pre-Neut JAK-STAT             6            5       2.53      3.90  
 6 MMP9 Pre-Neut TGFb                 6            5       0.392     1.47  
 7 MMP9 Pre-Neut VEGF                 6            5       0.276     1.11  
 8 MMP9 Pre-Neut p53                  6            5       0.167    -3.15  
 9 S100A12 Neut  EGFR                 4            5       0.832     0.0725
10 S100A12 Neut  Hypoxia              4            5       3.49      4.67  
# ℹ 23 more rows
# ℹ 4 more variables: delta_aci_minus_ctrl <dbl>, p <dbl>, star <chr>,
#   annot_sub <fct>

细胞级 / 样本级直方图：每个中性粒亚型 × 通路。

long_cell_neut <- dat_neut |>
  dplyr::select(annot_sub, group, dplyr::all_of(pw)) |>
  tidyr::pivot_longer(dplyr::all_of(pw), names_to = "pathway", values_to = "score") |>
  dplyr::mutate(
    annot_sub = factor(.data$annot_sub, levels = traj_levels),
    pathway = factor(.data$pathway, levels = pw),
    group = factor(.data$group, levels = c("CTRL", "ACI"))
  )
p_hist_cell_neut <- plot_score_hist(long_cell_neut, "annot_sub", bins = 40) +
  ggtitle("Cell-level")
ggsave(
  file.path(res_path, "04_neut_hist_cell_subtype.pdf"),
  p_hist_cell_neut, width = 22, height = 12, dpi = 600
)
p_hist_cell_neut

long_med_neut <- pw_med_sub |>
  tidyr::pivot_longer(dplyr::all_of(pw), names_to = "pathway", values_to = "score") |>
  dplyr::mutate(
    annot_sub = factor(.data$annot_sub, levels = traj_levels),
    pathway = factor(.data$pathway, levels = pw),
    group = factor(.data$group, levels = c("CTRL", "ACI"))
  )
p_hist_med_neut <- plot_score_hist(long_med_neut, "annot_sub", bins = 10) +
  ggtitle("Patient-level (sample median)")
ggsave(
  file.path(res_path, "04_neut_hist_sample_subtype.pdf"),
  p_hist_med_neut, width = 22, height = 12, dpi = 600
)
p_hist_med_neut

p_pw_d <- ggplot(
  dplyr::mutate(pw_sample_sub, feature = factor(feature, levels = pw)),
  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.8) +
  scale_fill_activity(pw_sample_sub$delta_aci_minus_ctrl) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  xlab(NULL) +
  ylab(NULL) +
  ggtitle("Neut lineage sample-level ACI − CTRL (Wilcoxon p)")

heatmap fill 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-3.74374 -0.32002  0.21467  0.02656  0.84169  2.16032 
|x| q80 / q90 / q95: 1.19 / 2.1 / 3.06 

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

mean_pw_g <- mean_by_scaled(dat_neut, pw, c("annot_sub", "group"))
long_g <- tidyr::pivot_longer(
  mean_pw_g, dplyr::all_of(pw), names_to = "feature", values_to = "mean"
) |>
  dplyr::mutate(
    annot_sub = factor(annot_sub, levels = traj_levels),
    group = factor(group, levels = c("CTRL", "ACI"))
  )
p_pw_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.4982 -0.5901 -0.1810 -0.1393  0.2833  4.3365 
|x| q80 / q90 / q95: 0.76 / 0.933 / 1.03 

ggsave(file.path(res_path, "04_neut_heatmap_stage_group.pdf"), p_pw_g, width = 12, height = 6, dpi = 600)
p_pw_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