Subtype: HSPC

造血干祖亚型：HSC、MEP、GMP、MDP。看最后的命名和 marker 点图。

Level-1 HSPC only (clusters 21/24/31). Do not pull Pre-B / Pre-Neut into this object. After inspecting, mixed clusters can be invert-subset in round 2; GMP barcodes that belong with neutrophils can be saved later.

Subset 后重新 NormalizeData：细胞集合变了，HVG / PCA / Harmony 必须重做；丢掉全局 scale.data 后也需要从 counts 重建 data。

pacman::p_load(Seurat, tidyverse, qs2, harmony, tictoc, cowplot, presto, SCpubr, scop)
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")

res_path <- "results/subtype/hspc"
obj_path <- "data/rdata/subtype"
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)

alpha <- .5
pt_size <- .5
raster <- FALSE
raster_dpi <- c(3000, 3000)
min_cutoff <- "q10"
max_cutoff <- "q90"
n_hvg <- 3000

test_genes <- qs_read("data/rdata/test_genes.qs2", nthreads = 10)
seurat_object <-
  qs_read("data/rdata/seurat_object_first_annot.qs2", nthreads = 10)

table(seurat_object$annot)

          B         Ery        HSPC Mono/Mac/DC        Neut         pDC 
       5201        6456        3705        5042       20978         570 
     Plasma    Pre-Neut        T/NK 
        780       13822       15792 

seurat_object <-
  subset(seurat_object, subset = annot %in% c("HSPC"))
dim(seurat_object)

[1] 38606  3705

table(seurat_object$annot)

HSPC 
3705 

table(seurat_object$harmony_clusters_2)

   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
   0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 
  17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 
   0    0    0    0 1477    0    0 1279    0    0    0    0    0    0  949    0 
  33   34   35   36   37   38 
   0    0    0    0    0    0 

if (length(Layers(seurat_object, search = "counts")) > 1) {
  seurat_object <- JoinLayers(seurat_object)
}
seurat_object[["RNA"]]$data <- NULL
seurat_object[["RNA"]]$scale.data <- NULL
DefaultLayer(seurat_object@assays$RNA)

[1] "counts"

Layers(seurat_object)

[1] "counts"

1. Normalize, HVG, Scale, PCA

细胞集合变了，重新找高变基因和降维。PCA 肘图是在选后面用几个主成分，不是生物学结论。

seurat_object <- NormalizeData(seurat_object, verbose = TRUE)
seurat_object <- FindVariableFeatures(
  seurat_object,
  nfeatures = n_hvg,
  verbose = TRUE
)
head(VariableFeatures(seurat_object), 20)

 [1] "CLC"    "IGKC"   "PRG2"   "DNTT"   "HBB"    "CAMP"   "IGHA2"  "RBMS3" 
 [9] "HBA2"   "HBA1"   "CXCL12" "PF4"    "HBD"    "GHR"    "TNC"    "NAV2"  
[17] "HDC"    "LTBP1"  "PPBP"   "LEPR"  

LabelPoints(
  plot = VariableFeaturePlot(seurat_object),
  points = head(VariableFeatures(seurat_object), 20),
  repel = TRUE
)

pca_features <- setdiff(
  VariableFeatures(seurat_object),
  test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr]
)
sum(VariableFeatures(seurat_object) %in% test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr])

[1] 83

seurat_object <- ScaleData(
  seurat_object,
  features = pca_features,
  assay = "RNA"
)
seurat_object <- RunPCA(
  seurat_object,
  npcs = 50,
  features = pca_features,
  verbose = TRUE,
  reduction.name = "pca"
)

num_pcs <- list()
num_pcs$pct <-
  seurat_object[["pca"]]@stdev / sum(seurat_object[["pca"]]@stdev) * 100
num_pcs$cumu <- cumsum(num_pcs$pct)
num_pcs$co1 <- which(num_pcs$cumu > 80 & num_pcs$pct < 5)[1]
num_pcs$co2 <- sort(
  which(
    (num_pcs$pct[seq_along(num_pcs$pct)[-length(num_pcs$pct)]] -
      num_pcs$pct[seq_along(num_pcs$pct)[-1]]) > 0.1
  ),
  decreasing = TRUE
)[1] + 1
num_pcs$pcs <- min(num_pcs$co1, num_pcs$co2)
num_pcs$co1

[1] 35

num_pcs$co2

[1] 15

num_pcs$pcs

[1] 15

ElbowPlot(seurat_object, reduction = "pca", ndims = 50)

2. Harmony

再按样本整合一次。看各样本是否混开。

tic("Harmony integration")
seurat_object <- RunHarmony(
  seurat_object,
  c("sample"),
  ncores = 15,
  plot_convergence = TRUE
)

toc()

Harmony integration: 1.669 sec elapsed

num_pcs <- list()
num_pcs$pct <-
  seurat_object[["harmony"]]@stdev /
  sum(seurat_object[["harmony"]]@stdev) * 100
num_pcs$cumu <- cumsum(num_pcs$pct)
num_pcs$co1 <- which(num_pcs$cumu > 90 & num_pcs$pct < 5)[1]
num_pcs$co2 <- sort(
  which(
    (num_pcs$pct[seq_along(num_pcs$pct)[-length(num_pcs$pct)]] -
      num_pcs$pct[seq_along(num_pcs$pct)[-1]]) > 0.1
  ),
  decreasing = TRUE
)[1] + 1
num_pcs$co1

[1] 43

num_pcs$co2

[1] 39

plot_harmony_umap <- function(obj, dims, file) {
  umap_temp <- RunUMAP(
    obj,
    dims = 1:dims,
    reduction = "harmony",
    reduction.name = "umap_harmony"
  )
  p <- plot_grid(
    ncol = 2,
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "sample",
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1),
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "group",
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1),
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "phase",
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1),
    DimPlot(
      umap_temp, reduction = "umap_harmony", group.by = "is_doublet",
      pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
    ) + theme(aspect.ratio = 1)
  )
  ggsave(file, plot = p, width = 14, height = 12, dpi = 600)
  p
}

plot_harmony_umap(seurat_object, 20, file.path(res_path, "01_dim_umap20.pdf"))

plot_harmony_umap(seurat_object, 30, file.path(res_path, "01_dim_umap30.pdf"))

plot_harmony_umap(seurat_object, 40, file.path(res_path, "01_dim_umap40.pdf"))

# 看过 20/30/40 后改这里
n_dims <- 40
seurat_object <- RunUMAP(
  seurat_object,
  dims = 1:n_dims,
  reduction = "harmony",
  reduction.name = "umap_harmony"
)

3. Clustering

分辨率扫描。这里主要用于找低质量或明显混合的簇。

seurat_object <- FindNeighbors(
  seurat_object,
  reduction = "harmony",
  dims = 1:n_dims
)
for (res in c(1, 2, 3, 4)) {
  print(res)
  set.seed(123)
  seurat_object <- FindClusters(
    seurat_object,
    cluster.name = paste0("harmony_clusters_", res),
    resolution = res,
    algorithm = 4
  )
  gc()
}

[1] 1

[1] 2

[1] 3

[1] 4

plot_grid(
  ncol = 2,
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_1", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_2", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_3", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_4", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend()
)

ggsave(
  file.path(res_path, "02_qc_clusters.pdf"),
  width = 14, height = 12, dpi = 600
)

plot_grid(
  ncol = 2,
  FeaturePlot(
    seurat_object, reduction = "umap_harmony", features = "nFeature_RNA",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  FeaturePlot(
    seurat_object, reduction = "umap_harmony", features = "percent_mt",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object, reduction = "umap_harmony", group.by = "is_doublet",
    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_4", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
)

Idents(seurat_object) <- "harmony_clusters_4"
VlnPlot(
  seurat_object,
  features = c("nFeature_RNA", "nCount_RNA", "percent_mt", "doublet_score"),
  ncol = 2, pt.size = 0
)

4. Markers

点图：点越大、颜色越深，表示该簇里越多细胞、表达越高。用来判断要不要丢掉某些簇。

marker_groups <- list(
  HSC = c("CD34", "PROM1", "SPINK2", "AVP", "CRHBP", "HLF"),
  CMP = c("PRSS57", "FLT3", "RNASE2", "KIT", "MYB"),
  GMP = c("AZU1", "MPO", "ELANE", "CTSG", "PRTN3"),
  MDP = c("CSF1R", "KLF4", "IRF8", "LGALS1", "SAMHD1", "LSP1", "SPI1"),
  MEP_MK = c("GATA1", "GATA2", "ITGA2B", "GP9", "PF4", "PPBP"),
  CLP_preB = c("IL7R", "DNTT", "MME", "VPREB1", "IGLL1", "PAX5"),
  cycling = c("MKI67", "TOP2A"),
  contaminants = c("MS4A1", "CD79A", "CD3D", "NKG7", "CA1", "HBB", "FCGR3B", "CD14")
)

Idents(seurat_object) <- "harmony_clusters_4"
do_DotPlot(
  sample = seurat_object,
  flip = FALSE,
  legend.title = "Average Expression",
  legend.title.face = "plain",
  features = marker_groups,
  legend.length = 10
)

ggsave(
  file.path(res_path, "02_qc_dotplot.pdf"),
  width = 16, height = 8, dpi = 600
)

Idents(seurat_object) <- "harmony_clusters_4"
table(Idents(seurat_object))

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
225 199 192 192 181 179 176 176 172 148 134 127 127 126 122 121 118 116 116 111 
 21  22  23  24  25  26  27  28  29 
107 103  90  77  74  65  61  42  28 

hspc_markers <- FindAllMarkers(seurat_object, only.pos = TRUE)

marker_table <-
  hspc_markers |>
  group_by(cluster) |>
  slice_head(n = 20) |>
  mutate(rank = paste0("Top_", row_number())) |>
  pivot_wider(id_cols = rank, names_from = cluster, values_from = gene)
marker_table

# A tibble: 20 × 30
   rank  `1`   `2`   `3`   `4`   `5`   `6`   `7`   `8`   `9`   `10`  `11`  `12` 
   <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
 1 Top_1 VCAN  DEFA1 ITGA… CRHBP ASPM  DEFA1 PRTN3 IFIT… PRTN3 H2AC… HDC   MMP8 
 2 Top_2 IFI30 DEFA3 XACT  AVP   PLK1  DEFA4 CTSG  S100… CTSG  RRM2  AKAP… LTF  
 3 Top_3 S100… BPI   GATA1 PRKG2 UBE2C DEFA3 PRSS… LTF   ELANE H1-5  MS4A2 CRIS…
 4 Top_4 SAMH… CD24  ENSG… ENSG… DLGA… SRGN  ELANE S100… AZU1  H4C3  EPAS1 PGLY…
 5 Top_5 S100… DEFA4 KLF1  RBPMS TOP2A SLPI  AZU1  C1or… MPO   H3C2  BACE2 CAMP 
 6 Top_6 SLC8… SLPI  ENSG… LINC… GTSE1 BPI   MPO   LCN2  CST7  TOP2A ENSG… LCN2 
 7 Top_7 KCTD… ABCA… CNRI… ADGR… CENPF AZU1  CST7  G0S2  SRGN  H2AC… PRG2  ARG1 
 8 Top_8 ANXA2 PDE4D GMPR  IL12… KIF14 ELANE RPS28 CAMP  RPS28 MYBL2 ENPP3 CEAC…
 9 Top_9 LGAL… DEFA… ANK1  ENSG… MKI67 CD24  RPS17 FTH1  DEFA4 LYZ   CLC   TFF3 
10 Top_… FCN1  CEAC… PPP1… HLF   CCNF  PDE4D RPL3… TYRO… RPS17 NUSA… CLEC… CYP4…
11 Top_… LGAL… CEAC… ZNF3… PREX2 KIF11 CEAC… RPL34 IFIT… RNAS… MS4A… TPSB2 MMP9 
12 Top_… FMN1  ATP8… C2or… PCDH9 TPX2  ATP8… RPS8  FCGR… RPL3… H2AC… CDK15 OLFM4
13 Top_… LYZ   PDZD8 PKIG  INPP… CENPE DEFA… NUCB2 S100… RPL2… MKI67 CIBA… CHI3…
14 Top_… MARC… SRGN  MPIG… CALN1 PRC1  MPO   RPL22 S100… RPL22 AURKB THSD… ABCA…
15 Top_… CYBB  MCU   GP1BB CHRM3 CDK1  ANXA3 SRGN  CMTM2 RPL7  TUBA… RHOX… S100…
16 Top_… FGR   ARHG… SMIM1 ANGP… H3C2  RNAS… RPL37 S100… RPL34 FOXM1 TPSA… S100…
17 Top_… TMTC2 ANXA3 PDLI… HOPX  NUSA… CTSG  RPL9  FCER… RPL37 CDCA5 PTPRG CEAC…
18 Top_… RBM47 HS3S… CRACD MMRN1 CDCA8 MS4A3 RPL39 ENSG… RPS8  DIAP… COL1… PADI4
19 Top_… CST3  SDK1  ST8S… NPR3  ANLN  PCOL… RPLP2 MMP9  RPL39 FAM1… SLC2… ANXA3
20 Top_… S100… ARHG… RYR3  CD34  HMMR  RPL39 RPS29 B2M   RPS29 PRC1  NDFI… FCAR 
# ℹ 17 more variables: `13` <chr>, `14` <chr>, `15` <chr>, `16` <chr>,
#   `17` <chr>, `18` <chr>, `19` <chr>, `20` <chr>, `21` <chr>, `22` <chr>,
#   `23` <chr>, `24` <chr>, `25` <chr>, `26` <chr>, `27` <chr>, `28` <chr>,
#   `29` <chr>

pred_file <- "celltypist/predictions_low.tsv"
if (file.exists(pred_file)) {
  pred <- read_tsv(pred_file, show_col_types = FALSE)
  barcode_col <- names(pred)[1]
  pred_vec <- setNames(pred$majority_voting, pred[[barcode_col]])
  seurat_object$celltypist_low <-
    unname(pred_vec[colnames(seurat_object)])
  DimPlot(
    seurat_object,
    reduction = "umap_harmony",
    group.by = "celltypist_low",
    label = TRUE, repel = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
}

5. Round 2 (optional)

去掉低质量或混合簇后重新降维。

按 QC 和 marker 去掉低质量或混合簇，再从 counts 重新 Normalize、高变基因、PCA、Harmony 和聚类。

drop_clusters <- c(11, 18, 27, 29)
# drop_clusters <- c("3", "7")

if (length(drop_clusters) == 0) {
  message("No clusters listed to drop; skip round 2.")
} else {
  n_before <- ncol(seurat_object)
  keep_cells <- colnames(seurat_object)[
    !(as.character(seurat_object$harmony_clusters_4) %in%
      as.character(drop_clusters))
  ]
  seurat_object <- subset(seurat_object, cells = keep_cells)
  message(
    sprintf("Dropped %s cells; %s remain.", n_before - ncol(seurat_object), ncol(seurat_object))
  )

  Idents(seurat_object) <- "sample"
  drop_cols <- grep(
    "^harmony_clusters_|^seurat_clusters$",
    colnames(seurat_object[[]]),
    value = TRUE
  )
  seurat_object@meta.data <-
    seurat_object@meta.data[, setdiff(colnames(seurat_object[[]]), drop_cols), drop = FALSE]

  seurat_object[["RNA"]]$data <- NULL
  seurat_object[["RNA"]]$scale.data <- NULL
  seurat_object <- NormalizeData(seurat_object, verbose = TRUE)
  seurat_object <- FindVariableFeatures(
    seurat_object, nfeatures = n_hvg, verbose = TRUE
  )
  pca_features <- setdiff(
    VariableFeatures(seurat_object),
    test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr]
  )
  seurat_object <- ScaleData(seurat_object, features = pca_features, assay = "RNA")
  seurat_object <- RunPCA(
    seurat_object, npcs = 50, features = pca_features,
    verbose = TRUE, reduction.name = "pca"
  )
  seurat_object <- RunHarmony(
    seurat_object, c("sample"), ncores = 15, plot_convergence = TRUE
  )
  seurat_object <- RunUMAP(
    seurat_object, dims = 1:n_dims, reduction = "harmony",
    reduction.name = "umap_harmony"
  )
  seurat_object <- FindNeighbors(
    seurat_object, reduction = "harmony", dims = 1:n_dims
  )
  for (res in c(1, 2, 3, 4)) {
    set.seed(123)
    seurat_object <- FindClusters(
      seurat_object,
      cluster.name = paste0("harmony_clusters_", res),
      resolution = res,
      algorithm = 4
    )
  }
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_4", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
}

5b. Check cleanup

检查清理后 UMAP 是否还剩明显的垃圾团或双细胞团。

确认清理后的 UMAP 不再有明显的低质量或双细胞团。

plot_grid(
  ncol = 2,
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_1", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_2", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_3", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap_harmony",
    group.by = "harmony_clusters_4", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend()
)

ggsave(
  file.path(res_path, "03_postqc_clusters.pdf"),
  width = 14, height = 12, dpi = 600
)

plot_grid(
  ncol = 2,
  FeaturePlot(
    seurat_object, reduction = "umap_harmony", features = "nFeature_RNA",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  FeaturePlot(
    seurat_object, reduction = "umap_harmony", features = "nCount_RNA",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  FeaturePlot(
    seurat_object, reduction = "umap_harmony", features = "percent_mt",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object, reduction = "umap_harmony", group.by = "is_doublet",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
)

ggsave(
  file.path(res_path, "03_postqc_metrics.pdf"),
  width = 14, height = 12, dpi = 600
)

Idents(seurat_object) <- "harmony_clusters_4"
do_DotPlot(
  sample = seurat_object,
  flip = FALSE,
  legend.title = "Average Expression",
  legend.title.face = "plain",
  features = marker_groups,
  legend.length = 10
)

ggsave(
  file.path(res_path, "03_postqc_dotplot.pdf"),
  width = 16, height = 8, dpi = 600
)

若仍有低质量簇，记下当前高分辨率 ID，下一节 subset 后重跑降维。

drop_clusters_r3 <- c()
# drop_clusters_r3 <- c("2", "5")

if (length(drop_clusters_r3) == 0) {
  message("No extra clusters to drop.")
} else {
  n_before <- ncol(seurat_object)
  keep_cells <- colnames(seurat_object)[
    !(as.character(seurat_object$harmony_clusters_4) %in%
      as.character(drop_clusters_r3))
  ]
  seurat_object <- subset(seurat_object, cells = keep_cells)
  message(
    sprintf("Dropped %s cells; %s remain.", n_before - ncol(seurat_object), ncol(seurat_object))
  )
}

6. Recluster for annotation

用清理后的细胞重新聚类，供最后命名。

用清理后的细胞做最终降维，并在 0.1 / 0.2 / 0.5 / 1 分辨率下聚类，供命名。

Idents(seurat_object) <- "sample"
drop_cols <- grep(
  "^harmony_clusters_|^seurat_clusters$",
  colnames(seurat_object[[]]),
  value = TRUE
)
if (length(drop_cols) > 0) {
  seurat_object@meta.data <-
    seurat_object@meta.data[, setdiff(colnames(seurat_object[[]]), drop_cols), drop = FALSE]
}

seurat_object[["RNA"]]$data <- NULL
seurat_object[["RNA"]]$scale.data <- NULL
seurat_object <- NormalizeData(seurat_object, verbose = TRUE)
seurat_object <- FindVariableFeatures(
  seurat_object, nfeatures = n_hvg, verbose = TRUE
)
pca_features <- setdiff(
  VariableFeatures(seurat_object),
  test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr]
)
seurat_object <- ScaleData(seurat_object, features = pca_features, assay = "RNA")
seurat_object <- RunPCA(
  seurat_object, npcs = 50, features = pca_features,
  verbose = TRUE, reduction.name = "pca"
)
tic("Harmony integration")
seurat_object <- RunHarmony(
  seurat_object, c("sample"), ncores = 15, plot_convergence = TRUE
)

toc()

Harmony integration: 1.403 sec elapsed

seurat_object <- RunUMAP(
  seurat_object, dims = 1:20, reduction = "harmony",
  reduction.name = "umap"
)
seurat_object <- FindNeighbors(
  seurat_object, reduction = "harmony", dims = 1:20
)
for (res in c(0.1, 0.2, 0.5, 1)) {
  print(res)
  set.seed(123)
  seurat_object <- FindClusters(
    seurat_object,
    cluster.name = paste0("harmony_clusters_", res),
    resolution = res,
    algorithm = 4
  )
  gc()
}

[1] 0.1

[1] 0.2

[1] 0.5

[1] 1

plot_grid(
  ncol = 2,
  DimPlot(
    seurat_object, reduction = "umap",
    group.by = "harmony_clusters_0.1", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap",
    group.by = "harmony_clusters_0.2", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap",
    group.by = "harmony_clusters_0.5", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend(),
  DimPlot(
    seurat_object, reduction = "umap",
    group.by = "harmony_clusters_1", label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1) + NoLegend()
)

ggsave(
  file.path(res_path, "04_recluster.pdf"),
  width = 14, height = 12, dpi = 600
)

Idents(seurat_object) <- "harmony_clusters_1"
do_DotPlot(
  sample = seurat_object,
  flip = FALSE,
  legend.title = "Average Expression",
  legend.title.face = "plain",
  features = marker_groups,
  legend.length = 10
)

# 改 ident.1 为要检查是否该拆开的 cluster
FindMarkers(
  seurat_object,
  group.by = "harmony_clusters_1",
  ident.1 = 2,
  only.pos = TRUE
) |>
  head(n = 20)

               p_val avg_log2FC pct.1 pct.2     p_val_adj
PRTN3  3.591580e-213  2.8391528 0.985 0.466 1.386565e-208
CTSG   1.590241e-209  3.3073055 0.951 0.444 6.139284e-205
ELANE  5.238630e-199  2.6919269 0.987 0.586 2.022425e-194
AZU1   1.700258e-187  2.2154108 0.993 0.681 6.564016e-183
MPO    6.948147e-171  1.9694734 0.996 0.673 2.682401e-166
CST7   1.357134e-162  2.5319611 0.898 0.420 5.239353e-158
RPS28  6.429260e-142  0.9514689 1.000 0.978 2.482080e-137
SRGN   1.136528e-124  1.5413160 0.989 0.938 4.387682e-120
RPS17  2.007103e-114  1.5332362 0.958 0.850 7.748622e-110
RPL37A  1.518745e-97  1.0173430 0.976 0.938  5.863268e-93
RPL34   2.957147e-87  0.7928197 0.996 0.971  1.141636e-82
RPL22   2.608549e-79  0.8408146 0.978 0.933  1.007056e-74
RPL37   3.020923e-77  0.7456606 0.996 0.975  1.166258e-72
NUCB2   4.454765e-77  1.6365659 0.805 0.600  1.719807e-72
RPS8    5.144409e-76  0.7229479 0.998 0.975  1.986050e-71
PRSS57  2.332996e-74  2.2658131 0.633 0.298  9.006766e-70
RNASE3  5.282848e-72  2.0485563 0.664 0.328  2.039496e-67
DEFA4   2.694822e-71  1.7773446 0.692 0.301  1.040363e-66
RPS29   3.792236e-68  0.8457849 0.985 0.943  1.464031e-63
RPL39   7.427143e-66  0.6886934 0.998 0.980  2.867323e-61

7. Markers for annotation

命名用的 marker 点图。对照教科书基因，给每个簇起亚型名。

人骨髓 HSPC 主要看 HSC / CMP / GMP / MDP；MEP、CLP/Pre-B、成熟髓系是混入或下游。对照下面 FeaturePlot，和 backup 用过的同一套基因。

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CD34", "PROM1", "SPINK2", "AVP", "CRHBP", "HLF"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_HSC.pdf"), width = 18, height = 10)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("PRSS57", "FLT3", "RNASE2", "KIT", "MYB", "GATA1"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_CMP.pdf"), width = 18, height = 10)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CSF1R", "KLF4", "IRF8", "LGALS1", "SAMHD1", "LSP1", "SPI1", "FLT3"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_MDP.pdf"), width = 18, height = 15)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("AZU1", "MPO", "ELANE", "CTSG", "PRTN3", "CSF3R"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_GMP.pdf"), width = 18, height = 10)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("FCGR1A", "CD36", "CEBPD", "CEBPB", "VCAN"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_MP.pdf"), width = 18, height = 10)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("HLA-DRA", "HLA-DRB1", "FCER1A"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_preDC.pdf"), width = 18, height = 5)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("GATA2", "ITGA2B", "GATA1", "GP9", "PF4", "PPBP"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_MEP_MK.pdf"), width = 18, height = 10)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("IL7R", "DNTT", "MME", "VPREB1", "IGLL1", "PAX5", "TCF4", "CD79A", "TCL1A"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_CLP_preB.pdf"), width = 18, height = 15)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 2,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("MKI67", "TOP2A"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_cycling.pdf"), width = 12, height = 5)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c(
    "MS4A1", "CD3D", "NKG7",
    "CA1", "HBB", "HBD",
    "FCGR3B", "CD14", "CAMP",
    "LTF", "LCN2", "NAMPT"
  ),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_contam.pdf"), width = 18, height = 20)

8. Annotation

最终亚型名字。后面组成和通路分析用这里的标签。

命名用的分辨率须与 cluster 对照图一致。

annot_cluster <- "harmony_clusters_0.1"
Idents(seurat_object) <- annot_cluster

Idents(seurat_object) <- annot_cluster
do_DotPlot(
  sample = seurat_object,
  flip = FALSE,
  legend.title = "Average Expression",
  legend.title.face = "plain",
  features = marker_groups,
  legend.length = 10
)

ggsave(
  file.path(res_path, "05_markers_dotplot.pdf"),
  width = 16, height = 8, dpi = 600
)

Idents(seurat_object) <- annot_cluster
table(Idents(seurat_object))

   1    2    3    4 
1259  863  661  583 

hspc_markers_annot <- FindAllMarkers(seurat_object, only.pos = TRUE)
marker_table <-
  hspc_markers_annot |>
  group_by(cluster) |>
  slice_head(n = 20) |>
  mutate(rank = paste0("Top_", row_number())) |>
  pivot_wider(id_cols = rank, names_from = cluster, values_from = gene)
marker_table

# A tibble: 20 × 5
   rank   `1`     `2`     `3`     `4`            
   <chr>  <chr>   <chr>   <chr>   <chr>          
 1 Top_1  SAMHD1  ELANE   IFITM2  HTR1F          
 2 Top_2  LGALS1  AZU1    LTF     EGFL7          
 3 Top_3  ANXA2   DEFA4   S100A12 FAM30A         
 4 Top_4  LYZ     SRGN    S100A8  NKAIN2         
 5 Top_5  VIM     CTSG    LCN2    ZNF385D        
 6 Top_6  TMSB10  MPO     S100A9  LAPTM4B        
 7 Top_7  SLC8A1  PRTN3   G0S2    CYTL1          
 8 Top_8  MS4A6A  CST7    CAMP    CD34           
 9 Top_9  S100A4  RNASE3  IFITM1  ST8SIA6        
10 Top_10 CST3    SLPI    MMP9    ANGPT1         
11 Top_11 VCAN    RPS28   PGLYRP1 CALN1          
12 Top_12 S100A10 CEACAM6 FCGR3B  NPR3           
13 Top_13 ANXA5   DEFA1   DEFA3   SPINK2         
14 Top_14 HDAC9   NUCB2   CMTM2   GATA2          
15 Top_15 ITGB2   RPS17   NAMPT   MMRN1          
16 Top_16 IFI30   MS4A3   FTH1    ENSG00000289364
17 Top_17 RBM47   DEFA3   S100A6  CHRM3          
18 Top_18 CRIP1   CD24    RGS2    GUCY1A1        
19 Top_19 MKI67   RPL22   IL1R2   DIPK1B         
20 Top_20 LY86    BPI     DEFA1   MEIS1          

cluster_map 的 key 必须对应 annot_cluster（这里用 0.1）。填完后再把 res 1 cluster 10 全部改成 MEP。

cluster_map <- c(
  "1" = "MDP",
  "2" = "GMP",
  "3" = "GMP",
  "4" = "HSC"
)

if (length(cluster_map) == 0) {
  message("cluster_map empty; skip subtype annotation.")
} else {
  cluster_ids <- as.character(
    FetchData(seurat_object, vars = annot_cluster)[, 1]
  )
  seurat_object$annot_subtype <- unname(cluster_map[cluster_ids])
  cl_1 <- as.character(seurat_object$harmony_clusters_1)
  seurat_object$annot_subtype[cl_1 == "10"] <- "MEP"
  annot_levels <- c("HSC", "MEP", "GMP", "MDP")
  seurat_object$annot_subtype <- factor(
    seurat_object$annot_subtype,
    levels = intersect(annot_levels, unique(seurat_object$annot_subtype))
  )
  table(seurat_object$annot_subtype, useNA = "ifany")
}

 HSC  MEP  GMP  MDP 
 404  179 1524 1259 

DimPlot(
  seurat_object, reduction = "umap",
  group.by = "annot_subtype", label = TRUE,
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
) + theme(aspect.ratio = 1)

ggsave(
  file.path(res_path, "06_annot_umap.pdf"),
  width = 8, height = 8, dpi = 600
)

plot_grid(
  ncol = 2,
  DimPlot(
    seurat_object, reduction = "umap", group.by = "annot_subtype",
    label = TRUE,
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object, reduction = "umap", group.by = "sample",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object, reduction = "umap", group.by = "group",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object, reduction = "umap", group.by = "phase",
    pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
)

ggsave(
  file.path(res_path, "06_annot_overview.pdf"),
  width = 14, height = 12, dpi = 600
)

Idents(seurat_object) <- "annot_subtype"
do_DotPlot(
  sample = seurat_object,
  flip = FALSE,
  legend.title = "Average Expression",
  legend.title.face = "plain",
  features = marker_groups,
  legend.length = 10
)

ggsave(
  file.path(res_path, "06_annot_dotplot.pdf"),
  width = 16, height = 8, dpi = 600
)

qs_save(
  seurat_object,
  file = file.path(obj_path, "seurat_object_hspc.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] scop_0.9.2         SCpubr_2.0.2.9000  presto_1.0.0       data.table_1.16.2 
 [5] cowplot_1.1.3      tictoc_1.2.1       harmony_1.2.3      Rcpp_1.0.13-1     
 [9] qs2_0.1.7          lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1     
[13] dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1       
[17] tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0    Seurat_5.2.1      
[21] 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] vroom_1.6.5             globals_0.16.3          lattice_0.22-6         
 [13] MASS_7.3-65             magrittr_2.0.3          limma_3.62.1           
 [16] plotly_4.10.4           rmarkdown_2.29          yaml_2.3.10            
 [19] thisutils_0.5.2         httpuv_1.6.15           sctransform_0.4.1      
 [22] spam_2.11-0             spatstat.sparse_3.1-0   reticulate_1.42.0      
 [25] pbapply_1.7-2           RColorBrewer_1.1-3      abind_1.4-8            
 [28] zlibbioc_1.52.0         Rtsne_0.17              GenomicRanges_1.58.0   
 [31] BiocGenerics_0.52.0     tweenr_2.0.3            circlize_0.4.16        
 [34] GenomeInfoDbData_1.2.13 IRanges_2.40.0          S4Vectors_0.44.0       
 [37] ggrepel_0.9.6           irlba_2.3.5.1           listenv_0.9.1          
 [40] spatstat.utils_3.1-1    goftest_1.2-3           RSpectra_0.16-2        
 [43] spatstat.random_3.3-2   fitdistrplus_1.2-1      parallelly_1.39.0      
 [46] codetools_0.2-20        RcppRoll_0.3.1          ggforce_0.4.2          
 [49] tidyselect_1.2.1        shape_1.4.6.1           UCSC.utils_1.2.0       
 [52] farver_2.1.2            viridis_0.6.5           matrixStats_1.4.1      
 [55] stats4_4.4.3            spatstat.explore_3.3-3  jsonlite_1.8.9         
 [58] GetoptLong_1.0.5        progressr_0.15.1        ggridges_0.5.6         
 [61] survival_3.8-3          iterators_1.0.14        systemfonts_1.1.0      
 [64] foreach_1.5.2           tools_4.4.3             ragg_1.3.3             
 [67] ica_1.0-3               glue_1.8.0              gridExtra_2.3          
 [70] xfun_0.49               GenomeInfoDb_1.42.0     withr_3.0.2            
 [73] fastmap_1.2.0           fansi_1.0.6             digest_0.6.37          
 [76] timechange_0.3.0        R6_2.5.1                mime_0.12              
 [79] textshaping_0.4.0       colorspace_2.1-1        scattermore_1.2        
 [82] tensor_1.5              spatstat.data_3.1-4     RhpcBLASctl_0.23-42    
 [85] utf8_1.2.4              generics_0.1.3          httr_1.4.7             
 [88] htmlwidgets_1.6.4       uwot_0.2.2              pkgconfig_2.0.3        
 [91] gtable_0.3.6            ComplexHeatmap_2.22.0   lmtest_0.9-40          
 [94] XVector_0.46.0          htmltools_0.5.8.1       dotCall64_1.2          
 [97] clue_0.3-66             scales_1.3.0            png_0.1-8              
[100] spatstat.univar_3.1-1   knitr_1.49              Signac_1.16.0          
[103] tzdb_0.4.0              reshape2_1.4.4          rjson_0.2.23           
[106] nlme_3.1-166            zoo_1.8-12              GlobalOptions_0.1.2    
[109] KernSmooth_2.23-26      vipor_0.4.7             parallel_4.4.3         
[112] miniUI_0.1.1.1          ggrastr_1.0.2           pillar_1.9.0           
[115] grid_4.4.3              vctrs_0.6.5             RANN_2.6.2             
[118] promises_1.3.2          stringfish_0.18.0       xtable_1.8-4           
[121] cluster_2.1.8.1         beeswarm_0.4.0          evaluate_1.0.1         
[124] cli_3.6.3               compiler_4.4.3          Rsamtools_2.22.0       
[127] rlang_1.1.4             crayon_1.5.3            leidenbase_0.1.31      
[130] future.apply_1.11.3     labeling_0.4.3          ggbeeswarm_0.7.2       
[133] plyr_1.8.9              stringi_1.8.4           BiocParallel_1.40.2    
[136] viridisLite_0.4.2       deldir_2.0-4            assertthat_0.2.1       
[139] Biostrings_2.74.0       munsell_0.5.1           lazyeval_0.2.2         
[142] spatstat.geom_3.3-4     pacman_0.5.1            Matrix_1.7-3           
[145] RcppHNSW_0.6.0          hms_1.1.3               patchwork_1.3.0        
[148] bit64_4.5.2             future_1.34.0           statmod_1.5.0          
[151] shiny_1.9.1             ROCR_1.0-11             igraph_2.1.2           
[154] RcppParallel_5.1.9      fastmatch_1.1-4         bit_4.5.0