pacman::p_load(
Seurat, tidyverse, qs2, cowplot, 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")Cell composition
各例里各类细胞占多少。占比是相对数。主模型在 scCODA；这里 Wilcoxon 和 miloR 为补充。
样本级细胞占比（描述 + Wilcoxon 敏感性分析），导出给 07_2_sccoda.ipynb 做 scCODA / tascCODA，再用 miloR 做邻域差异丰度。主结论看成分模型；Wilcoxon 和 Milo 是补充。group 以 CTRL 为对照。
先跑本脚本，再在 conda stroke_bm 里跑 codes/07_2_sccoda.ipynb。
res_path <- "results/composition"
tab_path <- "data/rdata/composition"
if (dir.exists(res_path)) {
old_res <- list.files(res_path, full.names = TRUE)
old_res <- old_res[!grepl("^03_s?ccoda|^03_tasc", basename(old_res))]
if (length(old_res) > 0) unlink(old_res, recursive = TRUE)
} else {
dir.create(res_path, recursive = TRUE)
}
if (!dir.exists(tab_path)) dir.create(tab_path, recursive = TRUE)n_dims <- 20seurat_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"))
cell_meta <- FetchData(
seurat_object,
vars = c("sample", "group", "annot_main", "annot_sub")
) |>
rownames_to_column("barcode")
sample_meta <- cell_meta |>
distinct(sample, group) |>
arrange(group, sample)
print(table(cell_meta$group, cell_meta$sample))
ACI-1 ACI-2 ACI-3 ACI-4 ACI-5 CTRL-1 CTRL-2 CTRL-3 CTRL-4 CTRL-5 CTRL-6
CTRL 0 0 0 0 0 1642 4851 3088 5615 1856 5686
ACI 7204 7246 8998 7217 7787 0 0 0 0 0 0print(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 3231. Sample proportions
堆叠图和箱线：每个点是一例里该类占全部细胞的百分之几。
每个点是一例里该类占全部细胞的比例。先看中性粒、Pre-Neut 和 HSPC 是否整体上移；组成推断以 07.2 的模型为准。
make_prop <- function(meta, cell_col) {
totals <- meta |>
dplyr::count(sample, group, name = "n_total")
meta |>
dplyr::count(sample, group, cell = .data[[cell_col]], name = "n") |>
complete(nesting(sample, group), cell, fill = list(n = 0)) |>
left_join(totals, by = c("sample", "group")) |>
mutate(
prop = n / n_total,
cell = factor(cell, levels = levels(meta[[cell_col]]))
)
}
prop_main <- make_prop(cell_meta, "annot_main")
prop_sub <- make_prop(cell_meta, "annot_sub")
write_csv(cell_meta, file.path(tab_path, "cell_meta.csv"))
write_csv(sample_meta, file.path(tab_path, "sample_meta.csv"))
write_csv(prop_main, file.path(tab_path, "prop_main.csv"))
write_csv(prop_sub, file.path(tab_path, "prop_sub.csv"))sample_order <- sample_meta$sample
ggplot(prop_main, aes(sample, prop, fill = cell)) +
geom_col(width = 0.9) +
scale_x_discrete(limits = sample_order) +
scale_y_continuous(expand = c(0, 0)) +
scale_fill_manual(values = annot_main_col) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "right"
) +
labs(x = NULL, y = "Proportion", fill = "annot_main")
ggsave(
file.path(res_path, "01_prop_main_bar.pdf"),
width = 10, height = 5, dpi = 600
)ggplot(prop_sub, aes(sample, prop, fill = cell)) +
geom_col(width = 0.9) +
scale_x_discrete(limits = sample_order) +
scale_y_continuous(expand = c(0, 0)) +
scale_fill_manual(values = annot_sub_col) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "right",
legend.key.size = unit(0.35, "cm"),
legend.text = element_text(size = 7)
) +
labs(x = NULL, y = "Proportion", fill = "annot_sub")
ggsave(
file.path(res_path, "01_prop_sub_bar.pdf"),
width = 12, height = 6, dpi = 600
)ggplot(prop_main, aes(group, prop, fill = group)) +
geom_boxplot(outlier.shape = NA, width = 0.6) +
geom_jitter(width = 0.12, size = 1.5, alpha = 0.8) +
scale_fill_manual(values = group_col) +
facet_wrap(~cell, scales = "free_y") +
theme_bw() +
theme(legend.position = "none") +
labs(x = NULL, y = "Proportion")
ggsave(
file.path(res_path, "01_prop_main_box.pdf"),
width = 12, height = 8, dpi = 600
)ggplot(prop_sub, aes(group, prop, fill = group)) +
geom_boxplot(outlier.shape = NA, width = 0.6) +
geom_jitter(width = 0.12, size = 1.2, alpha = 0.8) +
scale_fill_manual(values = group_col) +
facet_wrap(~cell, scales = "free_y") +
theme_bw() +
theme(legend.position = "none", strip.text = element_text(size = 8)) +
labs(x = NULL, y = "Proportion")
ggsave(
file.path(res_path, "01_prop_sub_box.pdf"),
width = 16, height = 12, dpi = 600
)2. Wilcoxon (sensitivity)
用 6 vs 5 个样本做 Wilcoxon，只是敏感性分析。主模型在 07.2。
6 vs 5 的 Wilcoxon 不考虑「占比加总为 100%」，只作对照。主结果看 07.2。
wilcox_props <- function(prop_df) {
prop_df |>
group_by(cell) |>
group_modify(function(d, key) {
ctrl <- d$prop[d$group == "CTRL"]
aci <- d$prop[d$group == "ACI"]
tibble(
n_ctrl = length(ctrl),
n_aci = length(aci),
median_ctrl = median(ctrl),
median_aci = median(aci),
delta_aci_minus_ctrl = median_aci - median_ctrl,
p = wilcox.test(aci, ctrl, exact = FALSE)$p.value
)
}) |>
ungroup() |>
mutate(p_adj = p.adjust(p, method = "BH")) |>
arrange(p)
}
wilcox_main <- wilcox_props(prop_main)
wilcox_sub <- wilcox_props(prop_sub)
write_csv(wilcox_main, file.path(res_path, "02_wilcox_main.csv"))
write_csv(wilcox_sub, file.path(res_path, "02_wilcox_sub.csv"))
wilcox_main# A tibble: 10 × 8
cell n_ctrl n_aci median_ctrl median_aci delta_aci_minus_ctrl p p_adj
<fct> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 B 6 5 0.128 0.0193 -0.109 0.00811 0.0135
2 T 6 5 0.312 0.0996 -0.212 0.00811 0.0135
3 pDC 6 5 0.0154 0.00236 -0.0131 0.00811 0.0135
4 Pre-… 6 5 0.0485 0.279 0.231 0.00811 0.0135
5 Neut 6 5 0.0305 0.500 0.469 0.00811 0.0135
6 Ery 6 5 0.182 0.00319 -0.178 0.00811 0.0135
7 Plas… 6 5 0.0213 0.00153 -0.0198 0.0137 0.0171
8 Mono… 6 5 0.0976 0.0438 -0.0538 0.0137 0.0171
9 NK 6 5 0.0459 0.0152 -0.0307 0.0828 0.0920
10 HSPC 6 5 0.0642 0.0552 -0.00897 0.315 0.315 wilcox_sub# A tibble: 34 × 8
cell n_ctrl n_aci median_ctrl median_aci delta_aci_minus_ctrl p p_adj
<fct> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 HSC 6 5 0.0104 0.00208 -0.00828 0.00811 0.0145
2 MEP 6 5 0.00644 0.000556 -0.00588 0.00811 0.0145
3 GMP 6 5 0.00772 0.0382 0.0305 0.00811 0.0145
4 Pro-B 6 5 0.0101 0.000831 -0.00931 0.00811 0.0145
5 Larg… 6 5 0.0146 0.00116 -0.0134 0.00811 0.0145
6 Smal… 6 5 0.0268 0.00167 -0.0252 0.00811 0.0145
7 Naiv… 6 5 0.0384 0.00485 -0.0335 0.00811 0.0145
8 Memo… 6 5 0.0234 0.00511 -0.0183 0.00811 0.0145
9 CD4 … 6 5 0.0368 0.0142 -0.0226 0.00811 0.0145
10 CD8 … 6 5 0.0641 0.0151 -0.0490 0.00811 0.0145
# ℹ 24 more rows3. Export note
把计数表交给下一篇 scCODA。
把每例各类细胞计数交给 scCODA / tascCODA。
4. miloR
不按我们起的名字，而按转录组邻居圈一块细胞再比丰度。方向应和占比大致同向，细节可以不完全一样。
邻域丰度：正 logFC（红）= ACI 更富集。SpatialFDR < 0.1 常当作「这片邻居可信」。应和占比图方向大致一致。
在 Harmony 邻域上检验 ACI vs CTRL。不依赖离散 cluster；邻域再用 annot_main / annot_sub 注释。
pacman::p_load(miloR)
DefaultAssay(seurat_object) <- "RNA"
sce <- as.SingleCellExperiment(seurat_object, assay = "RNA")
reducedDim(sce, "HARMONY") <-
Embeddings(seurat_object, "harmony")[, 1:n_dims, drop = FALSE]
reducedDim(sce, "UMAP") <- Embeddings(seurat_object, "umap")
colData(sce)$sample <- as.character(colData(sce)$sample)
colData(sce)$group <- factor(colData(sce)$group, levels = c("CTRL", "ACI"))
milo <- Milo(sce)
milo <- buildGraph(milo, k = 20, d = n_dims, reduced.dim = "HARMONY")
milo <- makeNhoods(
milo, prop = 0.05, k = 20, d = n_dims,
refined = TRUE, reduced_dims = "HARMONY"
)
milo <- countCells(
milo,
meta.data = data.frame(colData(milo)),
samples = "sample"
)
milo <- calcNhoodDistance(milo, d = n_dims, reduced.dim = "HARMONY")
design_df <- data.frame(colData(milo))[, c("sample", "group")] |>
distinct()
rownames(design_df) <- design_df$sample
design_df <- design_df[colnames(nhoodCounts(milo)), , drop = FALSE]
design_df$group <- factor(design_df$group, levels = c("CTRL", "ACI"))
da_results <- testNhoods(
milo,
design = ~group,
design.df = design_df,
reduced.dim = "HARMONY"
)
da_results <- annotateNhoods(milo, da_results, coldata_col = "annot_main")
da_results <- annotateNhoods(milo, da_results, coldata_col = "annot_sub")
milo <- buildNhoodGraph(milo)
write_csv(da_results, file.path(res_path, "04_milo_da.csv"))
table(sign(da_results$logFC), da_results$SpatialFDR < 0.1)
FALSE TRUE
-1 468 697
1 144 1339summary(da_results$SpatialFDR) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000103 0.0002320 0.0031989 0.1071384 0.0712841 0.9986929 da_results <- da_results |>
mutate(
da_dir = case_when(
SpatialFDR >= 0.1 | is.na(SpatialFDR) ~ "ns",
logFC > 0 ~ "ACI",
TRUE ~ "CTRL"
),
da_dir = factor(da_dir, levels = c("CTRL", "ACI", "ns"))
)
ggplot(da_results, aes(logFC, -log10(SpatialFDR))) +
geom_point(aes(color = da_dir), alpha = 0.7, size = 1.2) +
scale_color_manual(
values = c(group_col, ns = "grey75"),
breaks = c("CTRL", "ACI", "ns"),
labels = c("CTRL enriched", "ACI enriched", "ns")
) +
theme_bw() +
labs(color = NULL)
ggsave(
file.path(res_path, "04_milo_volcano.pdf"),
width = 7, height = 5, dpi = 600
)plotDAbeeswarm(da_results, group.by = "annot_main") +
scale_color_gradient2(
low = milo_div_col["low"],
mid = milo_div_col["mid"],
high = milo_div_col["high"]
) +
theme_bw()
ggsave(
file.path(res_path, "04_milo_beeswarm_main.pdf"),
width = 8, height = 6, dpi = 600
)plotDAbeeswarm(da_results, group.by = "annot_sub") +
scale_color_gradient2(
low = milo_div_col["low"],
mid = milo_div_col["mid"],
high = milo_div_col["high"]
) +
theme_bw()
ggsave(
file.path(res_path, "04_milo_beeswarm_sub.pdf"),
width = 8, height = 10, dpi = 600
)plotNhoodGraphDA(
milo, as.data.frame(da_results), alpha = 0.1, layout = "UMAP"
) +
scale_fill_gradient2(
low = milo_div_col["low"],
mid = milo_div_col["mid"],
high = milo_div_col["high"],
name = "logFC"
)
ggsave(
file.path(res_path, "04_milo_graph.pdf"),
width = 8, height = 8, dpi = 600
)qs_save(milo, file.path(tab_path, "milo.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 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] miloR_2.2.0 edgeR_4.4.0
[3] limma_3.62.1 SingleCellExperiment_1.28.1
[5] SummarizedExperiment_1.36.0 Biobase_2.66.0
[7] GenomicRanges_1.58.0 GenomeInfoDb_1.42.0
[9] IRanges_2.40.0 S4Vectors_0.44.0
[11] BiocGenerics_0.52.0 MatrixGenerics_1.18.0
[13] matrixStats_1.4.1 cowplot_1.1.3
[15] qs2_0.1.7 lubridate_1.9.3
[17] forcats_1.0.0 stringr_1.5.1
[19] dplyr_1.1.4 purrr_1.0.2
[21] readr_2.1.5 tidyr_1.3.1
[23] tibble_3.2.1 ggplot2_3.5.1
[25] tidyverse_2.0.0 Seurat_5.2.1
[27] SeuratObject_5.0.2 sp_2.1-4
loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.22 splines_4.4.3 later_1.4.1
[4] polyclip_1.10-7 fastDummies_1.7.4 lifecycle_1.0.4
[7] globals_0.16.3 lattice_0.22-6 vroom_1.6.5
[10] MASS_7.3-65 magrittr_2.0.3 plotly_4.10.4
[13] rmarkdown_2.29 yaml_2.3.10 httpuv_1.6.15
[16] sctransform_0.4.1 spam_2.11-0 spatstat.sparse_3.1-0
[19] reticulate_1.42.0 pbapply_1.7-2 RColorBrewer_1.1-3
[22] abind_1.4-8 zlibbioc_1.52.0 Rtsne_0.17
[25] ggraph_2.2.1 pracma_2.4.4 tweenr_2.0.3
[28] GenomeInfoDbData_1.2.13 ggrepel_0.9.6 irlba_2.3.5.1
[31] listenv_0.9.1 spatstat.utils_3.1-1 goftest_1.2-3
[34] RSpectra_0.16-2 spatstat.random_3.3-2 fitdistrplus_1.2-1
[37] parallelly_1.39.0 codetools_0.2-20 DelayedArray_0.32.0
[40] ggforce_0.4.2 tidyselect_1.2.1 UCSC.utils_1.2.0
[43] farver_2.1.2 ScaledMatrix_1.14.0 viridis_0.6.5
[46] spatstat.explore_3.3-3 jsonlite_1.8.9 BiocNeighbors_2.0.0
[49] tidygraph_1.3.1 progressr_0.15.1 ggridges_0.5.6
[52] survival_3.8-3 systemfonts_1.1.0 tools_4.4.3
[55] ragg_1.3.3 ica_1.0-3 Rcpp_1.0.13-1
[58] glue_1.8.0 gridExtra_2.3 SparseArray_1.6.0
[61] xfun_0.49 numDeriv_2016.8-1.1 withr_3.0.2
[64] fastmap_1.2.0 fansi_1.0.6 rsvd_1.0.5
[67] digest_0.6.37 timechange_0.3.0 R6_2.5.1
[70] mime_0.12 textshaping_0.4.0 colorspace_2.1-1
[73] scattermore_1.2 gtools_3.9.5 tensor_1.5
[76] spatstat.data_3.1-4 utf8_1.2.4 generics_0.1.3
[79] data.table_1.16.2 graphlayouts_1.2.1 httr_1.4.7
[82] htmlwidgets_1.6.4 S4Arrays_1.6.0 uwot_0.2.2
[85] pkgconfig_2.0.3 gtable_0.3.6 lmtest_0.9-40
[88] XVector_0.46.0 htmltools_0.5.8.1 dotCall64_1.2
[91] scales_1.3.0 png_0.1-8 spatstat.univar_3.1-1
[94] knitr_1.49 tzdb_0.4.0 reshape2_1.4.4
[97] nlme_3.1-166 cachem_1.1.0 zoo_1.8-12
[100] KernSmooth_2.23-26 vipor_0.4.7 parallel_4.4.3
[103] miniUI_0.1.1.1 pillar_1.9.0 grid_4.4.3
[106] vctrs_0.6.5 RANN_2.6.2 promises_1.3.2
[109] BiocSingular_1.22.0 stringfish_0.18.0 beachmat_2.22.0
[112] xtable_1.8-4 cluster_2.1.8.1 beeswarm_0.4.0
[115] evaluate_1.0.1 cli_3.6.3 locfit_1.5-9.10
[118] compiler_4.4.3 rlang_1.1.4 crayon_1.5.3
[121] future.apply_1.11.3 labeling_0.4.3 ggbeeswarm_0.7.2
[124] plyr_1.8.9 stringi_1.8.4 BiocParallel_1.40.2
[127] viridisLite_0.4.2 deldir_2.0-4 munsell_0.5.1
[130] lazyeval_0.2.2 spatstat.geom_3.3-4 pacman_0.5.1
[133] Matrix_1.7-3 RcppHNSW_0.6.0 hms_1.1.3
[136] patchwork_1.3.0 bit64_4.5.2 future_1.34.0
[139] statmod_1.5.0 shiny_1.9.1 ROCR_1.0-11
[142] memoise_2.0.1 igraph_2.1.2 RcppParallel_5.1.9
[145] bit_4.5.0