Subtype: Mono / Mac / DC

单核、cDC2、pDC。本数据没有典型巨噬细胞大类。

Level-1 Mono/Mac/DC + pDC。中性粒不放进来。低质量 cluster 按 harmony_clusters_4 在 round 2 丢掉。

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/myeloid"
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("Mono/Mac/DC", "pDC"))
dim(seurat_object)

[1] 38606  5612

table(seurat_object$annot)

Mono/Mac/DC         pDC 
       5042         570 

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 2720    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    0    0    0    0    0    0    0 1169 1153    0    0    0 
  33   34   35   36   37   38 
   0    0  570    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] "IGKV3-15" "IGHV4-61" "PTGDS"    "HBB"      "IGHG2"    "IGLC3"   
 [7] "CAMP"     "DEFA3"    "IGHV4-31" "C1QC"     "DEFA1"    "GZMB"    
[13] "HBA2"     "C1QB"     "MT1G"     "C1QA"     "TCL1A"    "IGLC2"   
[19] "HBA1"     "LTF"     

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] 118

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: 2.564 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] 37

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 = 18, 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(
  monocyte = c("CD14", "FCN1", "LYZ", "VCAN", "S100A8", "S100A9", "FCGR3A"),
  macrophage = c("CD68", "C1QA", "C1QB", "C1QC", "CSF1R"),
  cDC1 = c("CLEC9A", "XCR1", "BATF3"),
  cDC2 = c("CD1C", "FCER1A", "CLEC10A"),
  pDC = c("CLEC4C", "IL3RA", "LILRA4", "GZMB", "IRF8"),
  cycling = c("MKI67", "TOP2A"),
  HSPC = c("CD34", "SPINK2"),
  neut = c("FCGR3B", "CXCR2", "CAMP", "LTF"),
  lymphoid = c("CD3D", "MS4A1"),
  mast = c("CPA3", "TPSAB1")
)

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 = 10, 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 
333 306 282 278 274 259 259 217 214 213 206 206 198 190 189 181 172 171 168 158 
 21  22  23  24  25  26  27  28  29  30  31  32  33 
150 132 131 117 116  96  91  89  70  45  44  43  14 

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

marker_table <-
  myeloid_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 × 34
   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 RHEX  SEMA… CLEC… S100… ENSG… LCN2  RNAS… FCGR… IFIT… ENSG… CST3  PID1 
 2 Top_2 PLD4  ATF3  CD1C  S100… CYP1… LTF   PDE4D CDKN… MMP9  DYSF  RPL41 SASH1
 3 Top_3 JCHA… CD83  ENHO  S100… TNFA… CAMP  LYZ   ENSG… S100… VCAN  CD74  MTSS1
 4 Top_4 SPIB  GABA… FCER… S100… DYSF  DEFA1 LGAL… HES4  S100… PLBD1 HLA-… SPRE…
 5 Top_5 NIBA… CD30… CCSE… FTL   ENSG… DEFA3 CAVI… UICLM G0S2  SLC2… HLA-… CSF1R
 6 Top_6 ITM2C ENSG… CD1E  NDUF… VCAN  PGLY… TMTC2 LYPD2 IFIT… S100… HLA-… HLA-…
 7 Top_7 LILR… ENSG… DEPT… TAX1… ENSG… S100… RNAS… FMNL2 S100… SERP… TMSB… CPVL 
 8 Top_8 MZB1  ENSG… NDRG2 SMDT1 JUN   S100P SERP… SNED1 FCGR… MEGF9 RPL30 EPB4…
 9 Top_9 C12o… ENSG… HLA-… HNRN… FCAR  CEAC… MS4A… MTSS1 S100… MAPK… RPL28 HLA-…
10 Top_… LIME1 ATP2… LRP1B COPS9 F5    TFF3  CLU   RHOC  CMTM2 S100… HLA-… HLA-…
11 Top_… ENSG… NR4A3 HLA-… MAP3… RBP7  S100… SLC4… TCF7… S100… SERP… RPS27 CALH…
12 Top_… COBL… DSC2  FCGR… COX7B ENSG… SLPI  AZU1  VMO1  S100P AGFG1 RPL39 MAFB 
13 Top_… TCF4  TOM1  CLIC2 PARVG SOCS3 IFIT… PADI4 ENSG… CD177 AGTP… RPL32 HLA-…
14 Top_… CLEC… NR4A2 DTNA  PAIP2 MCEM… MMP9  SLCO… CD30… C1or… PADI4 RPS28 CAMK…
15 Top_… PPP1… NEAT1 AFF3  PJA2  SLC2… CLC   MS4A… CKB   FTH1  HP    RPS24 FCGR…
16 Top_… SERP… GK    ADAM… TAOK1 ENSG… S100… CES1  KNDC1 DUSP1 LYST  RPS23 CAMK1
17 Top_… GZMB  FCAR  HLA-… SKP1  FFAR2 MMP8  METT… SYT17 PGLY… MYO10 RPS2… STAR…
18 Top_… CCDC… ZNF3… HLA-… H4C3  TRIB1 MYL6  ENO1  MS4A7 MYL6  PDE4D ENSG… LGAL…
19 Top_… UGCG  CYP1… CIITA HM13  ENSG… OLFM4 RETN  LINC… LTF   ITGAM RPL12 SH3P…
20 Top_… LINC… RFX2  HLA-… NDUF… CYP1… DEFA4 MAPK… HERC… LCN2  NOS1… RPL37 SECT…
# ℹ 21 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>, `30` <chr>, `31` <chr>, `32` <chr>, `33` <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(6, 9, 19, 21, 22, 25, 26, 31, 32, 33)
# 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"
VlnPlot(
  seurat_object,
  features = c("nFeature_RNA", "nCount_RNA", "percent_mt", "doublet_score"),
  ncol = 2, pt.size = 0
)

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: 2.062 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
PDE4D     1.509939e-133  2.0433808 0.697 0.246 5.829272e-129
PADI4     1.062770e-113  1.2957722 0.763 0.302 4.102928e-109
SERPINB10  1.738286e-86  1.6129046 0.328 0.066  6.710826e-82
VCAN       3.433321e-78  0.9071115 0.980 0.678  1.325468e-73
RNASE2     1.967744e-77  1.4847781 0.539 0.210  7.596673e-73
MEGF9      1.687119e-74  1.0151882 0.892 0.560  6.513293e-70
MAPK14     3.427072e-72  1.2262084 0.738 0.427  1.323055e-67
TMTC2      6.892938e-72  1.3942802 0.691 0.347  2.661088e-67
ITGAM      5.452851e-70  1.0705275 0.688 0.334  2.105128e-65
LYST       1.245044e-67  0.9622545 0.876 0.608  4.806618e-63
PLBD1      1.976691e-66  1.0079216 0.809 0.492  7.631212e-62
DYSF       5.624018e-66  1.1506443 0.647 0.302  2.171208e-61
LYZ        2.506844e-63  0.7104906 0.970 0.866  9.677921e-59
S100A9     3.371459e-61  0.7003017 0.970 0.886  1.301586e-56
RETN       6.327301e-61  1.1431991 0.746 0.441  2.442718e-56
IPCEF1     2.457760e-60  1.2349729 0.611 0.296  9.488428e-56
LRMDA      3.312725e-59  1.1897628 0.848 0.535  1.278911e-54
HP         2.397681e-58  1.4594879 0.379 0.127  9.256487e-54
ANXA6      5.770098e-57  1.1291048 0.585 0.290  2.227604e-52
SERPINB1   7.782842e-57  0.7809587 0.873 0.632  3.004644e-52

7. Markers for annotation

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

骨髓单核/DC 常见 CD14 Mono / CD16 Mono / Mac / cDC1 / cDC2 / pDC。中性粒（FCGR3B/CAMP）不要并进来。

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CD14", "FCN1", "LYZ", "VCAN", "S100A8", "S100A9", "FCGR3A"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CD68", "C1QA", "C1QB", "C1QC", "CSF1R"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CLEC9A", "XCR1", "BATF3"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

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

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CLEC4C", "IL3RA", "LILRA4", "GZMB", "IRF8"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

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 = 18, height = 5)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("FCGR3B", "CXCR2", "CAMP", "LTF", "CD3D", "MS4A1", "CD34", "CPA3"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

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    5 
1925 1503  466  277  205 

myeloid_markers <- FindAllMarkers(seurat_object, only.pos = TRUE)
marker_table <-
  myeloid_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 × 6
   rank   `1`      `2`     `3`             `4`      `5`            
   <chr>  <chr>    <chr>   <chr>           <chr>    <chr>          
 1 Top_1  VCAN     S100A9  RHEX            CLEC10A  FCGR3A         
 2 Top_2  DPYD     FTL     AFF3            CD1C     CDKN1C         
 3 Top_3  NEAT1    RPLP1   ENSG00000225885 ENHO     ENSG00000287682
 4 Top_4  ARHGAP26 TMSB10  TCF4            CCSER1   HES4           
 5 Top_5  DYSF     RPL28   PLD4            FCER1A   UICLM          
 6 Top_6  ARHGAP24 RPL41   NIBAN3          CD1E     FMNL2          
 7 Top_7  ACSL1    RPL39   UGCG            DEPTOR   LYPD2          
 8 Top_8  FNDC3B   RBM8A   MZB1            NDRG2    SNED1          
 9 Top_9  FCAR     HNRNPH1 JCHAIN          HLA-DQA1 MTSS1          
10 Top_10 SLC2A3   RPS28   SPIB            CLIC2    RHOC           
11 Top_11 GAS7     IGKC    COBLL1          ADAM28   CD300H         
12 Top_12 CD14     CLTC    ITM2C           HLA-DPB1 TCF7L2         
13 Top_13 FMN1     CAPZA1  GZMB            AFF3     VMO1           
14 Top_14 WDFY3    TAPBP   PTPRS           HLA-DQB1 CKB            
15 Top_15 SLC11A1  RPS15   LINC01374       DTNA     ENSG00000289263
16 Top_16 ZDHHC20  TMSB4X  FHIP1A          FCGR2B   C1QA           
17 Top_17 AGFG1    OSBPL11 C12orf75        HLA-DPA1 KNDC1          
18 Top_18 MEGF9    LAMTOR5 LINC01478       LRP1B    ABI3           
19 Top_19 IRS2     RBM4    PPP1R14B        HLA-DRA  IL21R          
20 Top_20 CYP1B1   S100A8  CARD11          CD74     HERC2P3        

cluster_map 的 key 必须对应 annot_cluster。反复核对 FeaturePlot / DotPlot / marker 表后再填。

cluster_map <- c(
  "1" = "CD14 Mono",
  "2" = "CD14 Mono",
  "3" = "pDC",
  "4" = "cDC2",
  "5" = "CD16 Mono"
)

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])
  annot_levels <- c("CD14 Mono", "CD16 Mono", "cDC2", "pDC")
  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")
}

CD14 Mono CD16 Mono      cDC2       pDC 
     3428       205       277       466 

if (!"umap" %in% Reductions(seurat_object) &&
    "umap_harmony" %in% Reductions(seurat_object)) {
  seurat_object[["umap"]] <- seurat_object[["umap_harmony"]]
}

DimPlot(
  seurat_object,
  reduction = "umap",
  group.by = "annot_subtype",
  pt.size = pt_size,
  alpha = alpha,
  raster = raster,
  raster.dpi = raster_dpi,
  label = TRUE
) + 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_myeloid.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