pacman::p_load(Seurat, tidyverse, presto, cowplot, SCpubr, scop, qs2)
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/functions.R")Cell Type Annotation
按 marker 把 cluster 收成大类（HSPC、B、中性粒等）。CellTypist 只作参考。
alpha <- .03
pt_size <- .03
raster <- FALSE
raster_dpi <- c(3000, 3000)
min_cutoff <- "q10"
max_cutoff <- "q90"res_path <- "results/dim_reduction"
obj_path <- "data/rdata/dim_reduction"
if (!dir.exists(res_path)) dir.create(res_path, recursive = TRUE)
if (!dir.exists(obj_path)) dir.create(obj_path, recursive = TRUE)
if (!dir.exists("results/figures")) dir.create("results/figures", recursive = TRUE)seurat_object <-
qs_read(
"data/rdata/dim_reduction/seurat_object_clustered.qs2",
nthreads = 10
)
Layers(seurat_object)[1] "counts" "data" "scale.data"DefaultLayer(seurat_object@assays$RNA)[1] "counts"resolutions <- c(0.1, 0.2, 0.5, 1, 1.5, 2)
sc_markers <- list()
for (res in resolutions) {
col_name <- paste0("harmony_clusters_", res)
Idents(seurat_object) <- col_name
cat(sprintf("\n=== Resolution %s Cell Counts ===\n", res))
print(table(Idents(seurat_object)))
sc_markers[[paste0("all_markers_", res)]] <- FindAllMarkers(
seurat_object,
only.pos = TRUE
)
}
=== Resolution 0.1 Cell Counts ===
1 2 3 4 5 6 7 8 9
21124 15795 11443 10354 7046 4188 1048 812 536
=== Resolution 0.2 Cell Counts ===
1 2 3 4 5 6 7 8 9 10 11 12 13
21071 11521 7303 5363 5022 5001 4795 3463 3024 2259 2176 812 536
=== Resolution 0.5 Cell Counts ===
1 2 3 4 5 6 7 8 9 10 11 12 13
12513 9361 8546 5108 5006 4655 4287 4143 3455 2490 2306 2260 2176
14 15 16 17 18
1971 1636 1051 845 537
=== Resolution 1 Cell Counts ===
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
8065 7791 5126 5003 4985 4652 4146 4040 3455 2733 2534 2373 2330 2306 2176 1973
17 18 19 20 21 22 23 24 25 26
1744 1346 1275 980 844 628 538 514 420 369
=== Resolution 1.5 Cell Counts ===
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
5313 5175 4804 4421 4156 3814 3450 3340 3225 2852 2725 2695 2657 2527 2206 2176
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
1973 1813 1746 1709 1465 1306 1167 1003 944 780 633 570 514 421 397 369
=== Resolution 2 Cell Counts ===
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
4118 4060 3900 3783 3382 3210 3132 3131 2723 2720 2664 2557 2501 2318 2043 1970
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
1746 1583 1551 1520 1477 1466 1300 1279 1266 1262 1196 1169 1153 989 949 874
33 34 35 36 37 38
780 636 570 521 426 421 qs_save(sc_markers, file.path(obj_path, "sc_markers.qs2"), nthreads = 10)plot_grid(
ncol = 3,
DimPlot(
seurat_object,
reduction = "umap_harmony",
group.by = "harmony_clusters_0.1",
pt.size = pt_size,
alpha = alpha,
raster = raster,
raster.dpi = raster_dpi
) + theme(aspect.ratio = 1),
DimPlot(
seurat_object,
reduction = "umap_harmony",
group.by = "harmony_clusters_0.2",
pt.size = pt_size,
alpha = alpha,
raster = raster,
raster.dpi = raster_dpi
) + theme(aspect.ratio = 1),
DimPlot(
seurat_object,
reduction = "umap_harmony",
group.by = "harmony_clusters_0.5",
pt.size = pt_size,
alpha = alpha,
raster = raster,
raster.dpi = raster_dpi
) + theme(aspect.ratio = 1),
DimPlot(
seurat_object,
reduction = "umap_harmony",
group.by = "harmony_clusters_1",
pt.size = pt_size,
alpha = alpha,
raster = raster,
raster.dpi = raster_dpi
) + theme(aspect.ratio = 1),
DimPlot(
seurat_object,
reduction = "umap_harmony",
group.by = "harmony_clusters_1.5",
pt.size = pt_size,
alpha = alpha,
raster = raster,
raster.dpi = raster_dpi
) + theme(aspect.ratio = 1),
DimPlot(
seurat_object,
reduction = "umap_harmony",
group.by = "harmony_clusters_2",
pt.size = pt_size,
alpha = alpha,
raster = raster,
raster.dpi = raster_dpi
) + theme(aspect.ratio = 1)
)
# 与下方 cluster_map 对齐：marker 表和 cluster 映射必须用同一个分辨率（res 2）
marker_table <-
sc_markers[["all_markers_2"]] |>
group_by(cluster) |>
slice_head(n = 20) |>
group_by(cluster) |>
mutate(rank = paste0("Top_", row_number())) |>
pivot_wider(id_cols = rank, names_from = cluster, values_from = gene)
# View(marker_table)cluster_map <- c(
"1" = "T/NK",
"2" = "Neut",
"3" = "Neut",
"4" = "T/NK",
"5" = "Neut",
"6" = "Neut",
"7" = "T/NK",
"8" = "Pre-Neut",
"9" = "Ery",
"10" = "Mono/Mac/DC",
"11" = "Pre-Neut",
"12" = "Neut",
"13" = "Pre-Neut",
"14" = "Neut",
"15" = "T/NK",
"16" = "B",
"17" = "Ery",
"18" = "Pre-Neut",
"19" = "Neut",
"20" = "T/NK",
"21" = "HSPC",
"22" = "Ery",
"23" = "B",
"24" = "HSPC",
"25" = "Pre-Neut",
"26" = "Pre-Neut",
"27" = "T/NK",
"28" = "Mono/Mac/DC",
"29" = "Mono/Mac/DC",
"30" = "Pre-Neut",
"31" = "HSPC",
"32" = "B",
"33" = "Plasma",
"34" = "B",
"35" = "pDC",
"36" = "Ery",
"37" = "Pre-Neut",
"38" = "B"
)
seurat_object$annot <-
unname(cluster_map[as.character(seurat_object$harmony_clusters_2)])DimPlot(
seurat_object,
reduction = "umap_harmony",
group.by = "annot",
pt.size = .1, alpha = .1,
raster = FALSE, label = TRUE
) + theme(aspect.ratio = 1)
ggsave("results/dim_reduction/umap_harmony_first_annot.pdf",
width = 8, height = 8, dpi = 600)qs_save(
seurat_object,
file = "data/rdata/seurat_object_first_annot.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] stats graphics grDevices utils datasets methods base
other attached packages:
[1] qs2_0.1.7 scop_0.9.2 SCpubr_2.0.2.9000 cowplot_1.1.3
[5] presto_1.0.0 data.table_1.16.2 Rcpp_1.0.13-1 lubridate_1.9.3
[9] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
[13] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
[17] tidyverse_2.0.0 Seurat_5.2.1 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] bitops_1.0-9 polyclip_1.10-7 fastDummies_1.7.4
[7] lifecycle_1.0.4 thisplot_0.4.6 doParallel_1.0.17
[10] globals_0.16.3 lattice_0.22-6 MASS_7.3-65
[13] magrittr_2.0.3 plotly_4.10.4 rmarkdown_2.29
[16] yaml_2.3.10 thisutils_0.5.2 httpuv_1.6.15
[19] sctransform_0.4.1 spam_2.11-0 spatstat.sparse_3.1-0
[22] reticulate_1.42.0 pbapply_1.7-2 RColorBrewer_1.1-3
[25] abind_1.4-8 zlibbioc_1.52.0 Rtsne_0.17
[28] GenomicRanges_1.58.0 BiocGenerics_0.52.0 tweenr_2.0.3
[31] circlize_0.4.16 GenomeInfoDbData_1.2.13 IRanges_2.40.0
[34] S4Vectors_0.44.0 ggrepel_0.9.6 irlba_2.3.5.1
[37] listenv_0.9.1 spatstat.utils_3.1-1 goftest_1.2-3
[40] RSpectra_0.16-2 spatstat.random_3.3-2 fitdistrplus_1.2-1
[43] parallelly_1.39.0 codetools_0.2-20 RcppRoll_0.3.1
[46] ggforce_0.4.2 tidyselect_1.2.1 shape_1.4.6.1
[49] UCSC.utils_1.2.0 farver_2.1.2 matrixStats_1.4.1
[52] stats4_4.4.3 spatstat.explore_3.3-3 jsonlite_1.8.9
[55] GetoptLong_1.0.5 progressr_0.15.1 ggridges_0.5.6
[58] survival_3.8-3 iterators_1.0.14 systemfonts_1.1.0
[61] foreach_1.5.2 tools_4.4.3 ragg_1.3.3
[64] ica_1.0-3 glue_1.8.0 gridExtra_2.3
[67] xfun_0.49 GenomeInfoDb_1.42.0 withr_3.0.2
[70] fastmap_1.2.0 fansi_1.0.6 digest_0.6.37
[73] timechange_0.3.0 R6_2.5.1 mime_0.12
[76] textshaping_0.4.0 colorspace_2.1-1 scattermore_1.2
[79] tensor_1.5 spatstat.data_3.1-4 utf8_1.2.4
[82] generics_0.1.3 httr_1.4.7 htmlwidgets_1.6.4
[85] uwot_0.2.2 pkgconfig_2.0.3 gtable_0.3.6
[88] ComplexHeatmap_2.22.0 lmtest_0.9-40 XVector_0.46.0
[91] htmltools_0.5.8.1 dotCall64_1.2 clue_0.3-66
[94] scales_1.3.0 png_0.1-8 spatstat.univar_3.1-1
[97] knitr_1.49 Signac_1.16.0 tzdb_0.4.0
[100] reshape2_1.4.4 rjson_0.2.23 nlme_3.1-166
[103] zoo_1.8-12 GlobalOptions_0.1.2 KernSmooth_2.23-26
[106] parallel_4.4.3 miniUI_0.1.1.1 pillar_1.9.0
[109] grid_4.4.3 vctrs_0.6.5 RANN_2.6.2
[112] promises_1.3.2 stringfish_0.18.0 xtable_1.8-4
[115] cluster_2.1.8.1 evaluate_1.0.1 cli_3.6.3
[118] compiler_4.4.3 Rsamtools_2.22.0 rlang_1.1.4
[121] crayon_1.5.3 future.apply_1.11.3 labeling_0.4.3
[124] plyr_1.8.9 stringi_1.8.4 viridisLite_0.4.2
[127] deldir_2.0-4 BiocParallel_1.40.2 munsell_0.5.1
[130] Biostrings_2.74.0 lazyeval_0.2.2 spatstat.geom_3.3-4
[133] pacman_0.5.1 Matrix_1.7-3 RcppHNSW_0.6.0
[136] hms_1.1.3 patchwork_1.3.0 future_1.34.0
[139] shiny_1.9.1 ROCR_1.0-11 igraph_2.1.2
[142] RcppParallel_5.1.9 fastmatch_1.1-4