Subtype: B / Plasma

B 细胞和浆细胞亚型。PCA 不含 IG 基因，避免克隆型主导分群。

Level-1 B + Plasma。Scale / PCA / clustering 排除 IG 基因，避免克隆型 IG 主导分群。低质量 / doublet 在 subtype 按 harmony_clusters_4 丢掉。

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/b_plasma"
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("B", "Plasma"))
dim(seurat_object)

[1] 38606  5981

table(seurat_object$annot)

     B Plasma 
  5201    780 

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 1970 
  17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 
   0    0    0    0    0    0 1300    0    0    0    0    0    0    0    0  874 
  33   34   35   36   37   38 
 780  636    0    0    0  421 

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-20"  "IGHV3-7"   "IGKV1-39"  "IGHV3-23"  "HBB"       "IGHG2"    
 [7] "IGHA2"     "HBA2"      "IGHG3"     "HBA1"      "IGHV1-69D" "IGKV2D-28"
[13] "IGKV4-1"   "IGHGP"     "IGHG4"     "CAMP"      "IGKV1-5"   "IGHV4-39" 
[19] "IGKV1D-33" "HBD"      

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]
)
sum(VariableFeatures(seurat_object) %in% test_genes$gene[test_genes$unwant])

[1] 19

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

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.626 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] 42

num_pcs$co2

[1] 27

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(
  progenitor = c("CD34", "VPREB1", "IGLL1", "EBF1", "PAX5", "DNTT"),
  mature_B = c("MS4A1", "CD19", "CD79A", "CD79B", "CD22", "CD74"),
  plasma = c("MZB1", "SDC1", "JCHAIN", "XBP1", "PRDM1"),
  cycling = c("MKI67", "TOP2A", "PCNA"),
  HSPC = c("CD34", "SPINK2", "AVP"),
  T_NK = c("CD3D", "NKG7"),
  myeloid = c("CD14", "FCN1", "MPO", "FCGR3B", "CA1")
)

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 
401 378 344 329 323 287 266 263 253 245 230 212 202 197 194 191 188 179 178 171 
 21  22  23  24  25  26  27  28  29  30  31 
142 131 123 114 112 109  61  54  46  45  13 

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

marker_table <-
  b_plasma_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 × 32
   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 LINC… IGF2  DNTT  DTX1  FCER2 FKBP… AIM2  RRM2  CFAP… COL1… LINC… S100…
 2 Top_2 AIM2  ACSM3 LINC… CCDC… COL1… SLAM… SCIMP TOP2A IGLL1 IL4R  TMPR… S100…
 3 Top_3 ARHG… ENSG… ERG   PLD4  IL4R  JSRP1 SSPN  NUSA… RAG2  FCER2 LIMS2 S100…
 4 Top_4 NAPSB DLG2  SOCS2 FCRLA LINC… SDC1  SHIS… MKI67 ROR1  SYN3  MS4A1 TYRO…
 5 Top_5 CCSE… PCDH9 SCHI… SPIB  ADAM… TXND… TEX9  UBE2C SOX4  ADAM… PLD4  S100…
 6 Top_6 BANK1 BMP3  ENSG… NIBA… PLPP5 CHPF  COCH  CENPF VPRE… NETO1 ENSG… MMP9 
 7 Top_7 CTSH  LINC… LCN6  CD72  RPL13 PRDX4 PLEK… TYMS  ENSG… LARG… SUGCT S100…
 8 Top_8 TBC1… STK3… DEPP1 MYO1C RPL1… SPAG4 ENSG… AURKB CMTM8 ADAM… FCRL1 G0S2 
 9 Top_9 RPS27 CDC2… EGFL7 BMP3  BTLA  DERL3 CD86  H3C2  CEMIP LIX1… NCAM2 S100P
10 Top_… LTB   RAG1  CD34  CD79B HVCN1 TNFR… LINC… ASPM  ACSM3 LINC… GNG3  ENSG…
11 Top_… MAST4 SMAD1 LINC… TCL1A RPS27 SEC1… SNED1 PCLAF MYB   HVCN1 PPP1… S100…
12 Top_… RHEX  CCDC… ENSG… KREM… MS4A1 AQP3  ANK3  DIAP… IGLL5 CCR7  CYB5… FCER…
13 Top_… OSBP… LAMP5 SOCS… SAG   CD52  BHLH… BTNL9 H2AX  IRF4  STEA… CD79A CMTM2
14 Top_… COTL1 LHPP  SH2D… KLHL… RPL21 ANKR… BANK1 SPC25 P4HA2 MARC… ENSG… LST1 
15 Top_… RPL41 ENSG… ENSG… PCDH9 BANK1 FKBP2 LTB   KIF11 SNHG7 ZBTB… ALOX5 CDA  
16 Top_… RPL36 LDLR… MIR4… MS4A… LIX1… SELE… RPS14 TPX2  WASF1 PLPP5 RALG… IFIT…
17 Top_… GPR1… ENSG… RBMS3 ENSG… CD74  DNAJ… CRIP2 H1-5  CDHR3 SNTG2 PTPRJ MNDA 
18 Top_… RPL30 YBX3  OR7C1 VPRE… SELL  MEI1  RPS25 CDK1  BACH2 THRB  SPIB  LYZ  
19 Top_… RPS2… ENSG… LEF1… ENSG… CCR7  TENT… RPS29 EZH2  TMEM… ADK   FCRLA PGLY…
20 Top_… RPS21 NSMC… LINC… ENSG… ADAM… FNDC… RPS27 KNL1  ADAM… BANK1 TCL1A CSF3R
# ℹ 19 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>

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(12, 13, 16, 19, 27, 28, 30, 31)
# 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]
  )
  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(24)
# 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]
)
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.215 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
AIM2             0.000000e+00  4.2305499 0.534 0.038  0.000000e+00
COTL1           9.399305e-285  2.6147518 0.786 0.181 3.628696e-280
ARHGAP24        4.249259e-265  2.5206439 0.778 0.186 1.640469e-260
ENSG00000258760 8.908664e-260  5.3913071 0.343 0.013 3.439279e-255
BANK1           2.746154e-254  1.9751628 0.965 0.410 1.060180e-249
GPR183          2.762831e-233  3.4868157 0.441 0.045 1.066619e-228
NAPSB           6.270044e-226  2.3581504 0.657 0.132 2.420613e-221
SCIMP           2.520761e-219  2.8120868 0.498 0.067 9.731651e-215
OSBPL10         1.352712e-217  2.3076359 0.735 0.202 5.222281e-213
LTB             7.364062e-216  1.8094205 0.949 0.410 2.842970e-211
RPS27           3.336190e-205  1.1013035 0.999 0.966 1.287969e-200
RPS29           1.713933e-202  1.4548436 0.978 0.876 6.616809e-198
CTSH            2.097326e-201  2.4816145 0.583 0.117 8.096938e-197
ANK3            3.644819e-200  3.6638078 0.350 0.029 1.407119e-195
RPL41           4.273426e-197  0.8958169 0.996 0.979 1.649799e-192
TNFRSF13B       3.067407e-196  3.3424824 0.449 0.066 1.184203e-191
RPS14           9.808583e-189  0.9508525 0.994 0.950 3.786702e-184
RPL36           4.066437e-188  1.0711884 0.988 0.908 1.569889e-183
LINC01857       1.925334e-187  3.2046947 0.394 0.047 7.432943e-183
B2M             8.117964e-187  0.7567039 0.997 0.985 3.134021e-182

7. Markers for annotation

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

骨髓 B 系常见 Pro-B / Pre-B → naive/memory B → Plasma。Scale/PCA 已排除 IG。看 progenitor vs MS4A1 vs MZB1 决定要不要拆，还是 B + Plasma 两档即可。

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

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("MS4A1", "CD19", "CD79A", "CD79B", "CD22", "CD74"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("MZB1", "SDC1", "JCHAIN", "XBP1", "PRDM1"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("MKI67", "TOP2A", "PCNA"),
  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("CD3D", "NKG7", "CD14", "CD68", "CD34", "CA1", "FCGR3B", "MPO"),
  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

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

用分辨率 0.5 命名。去掉 res 0.1 cluster 1 且 res 0.5 cluster 7 的细胞；res 0.1 cluster 2 记为 Naive B。

cl_01 <- as.character(seurat_object$harmony_clusters_0.1)
cl_05 <- as.character(seurat_object$harmony_clusters_0.5)
drop_mixed <- cl_01 == "1" & cl_05 == "7"
keep_cells <- colnames(seurat_object)[!drop_mixed]
message(sprintf(
  "Drop 0.1=1 & 0.5=7: %s cells; keep %s.",
  sum(drop_mixed),
  length(keep_cells)
))
seurat_object <- subset(seurat_object, cells = keep_cells)

annot_cluster <- "harmony_clusters_0.5"
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    6    7    8 
1059  868  809  539  515  429  269  371 

b_plasma_markers <- FindAllMarkers(seurat_object, only.pos = TRUE)
marker_table <-
  b_plasma_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 × 9
   rank   `1`             `2`       `3`            `4`   `5`   `6`   `7`   `8`  
   <chr>  <chr>           <chr>     <chr>          <chr> <chr> <chr> <chr> <chr>
 1 Top_1  ACSM3           FCER2     COTL1          MKI67 LINC… FKBP… FKBP… DNTT 
 2 Top_2  CCDC191         IL4R      AIM2           TYMS  PLD4  TXND… TXND… LINC…
 3 Top_3  PCDH9           COL19A1   BANK1          NUSA… CCDC… DERL3 SSR4  SOCS2
 4 Top_4  LDLRAD4         LINC00926 ARHGAP24       PCLAF SPIB  SEC1… PRDX4 ERG  
 5 Top_5  RUBCNL          HVCN1     SCIMP          TOP2A FCRLA FNDC… DERL3 SCHI…
 6 Top_6  BACH2           PLPP5     GPR183         LMNB1 ENSG… PRDX4 HBB   ENSG…
 7 Top_7  MME             ADAM28    NAPSB          RRM2  MYO1C TENT… IGHA1 LCN6 
 8 Top_8  IRF4            RPL18A    ENSG000002587… CENPF NIBA… SLAM… JCHA… EGFL7
 9 Top_9  AFF3            SELL      ANK3           SMC4  PTPRJ SDC1  NDUF… SOCS…
10 Top_10 ENSG00000224610 RPL13     TFEC           IGLL1 TCL1A MEI1  PRDX6 CD34 
11 Top_11 ZEB2            HLA-DRB1  OSBPL10        PTTG1 TMPR… JSRP1 SEC1… LINC…
12 Top_12 ADAM23          MS4A1     B2M            UBE2C MS4A… ANKR… SOD1  RBMS3
13 Top_13 WASF1           CD52      TNFRSF13B      H2AX  CD79B SPAG4 COX5B MIR4…
14 Top_14 CD38            BANK1     RPS29          HMGB2 CD72  CHPF  S100… DEPP1
15 Top_15 YBX3            LIX1-AS1  MS4A1          BIRC5 DTX1  AQP3  POLR… ENSG…
16 Top_16 BCL7A           CD74      RPS27          EZH2  PCDH9 BMP6  TUFM  ARL5C
17 Top_17 ACSM1           HLA-DQA1  RPL41          UHRF1 BMP3  FKBP2 TTC1  CYGB 
18 Top_18 PDE4D           ADAM7-AS1 MARCHF1        SMC2  PLEK… BHLH… FIS1  SH2D…
19 Top_19 RAPGEF5         IGHD      CD82           ZWINT CYB5… GAS6  PTK2B VPRE…
20 Top_20 CDC25B          MARCHF1   RPL36          CKS2  ENSG… PRDM1 MED10 MYLK 

cluster_map <- c(
  "1" = "Small Pre-B",
  "2" = "Naive B",
  "3" = "Memory B",
  "4" = "Large Pre-B",
  "5" = "Trans B",
  "6" = "Plasma",
  "7" = "Plasma",
  "8" = "Pro-B"
)

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_01 <- as.character(
    FetchData(seurat_object, vars = "harmony_clusters_0.1")[, 1]
  )
  seurat_object$annot_subtype[cl_01 == "2"] <- "Naive B"
  annot_levels <- c(
    "Pro-B", "Large Pre-B", "Small Pre-B", "Trans B", "Naive B",
    "Memory B", "Plasma"
  )
  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")
  table(
    FetchData(seurat_object, vars = "harmony_clusters_0.1")[, 1],
    seurat_object$annot_subtype,
    useNA = "ifany"
  )
}

   
    Pro-B Large Pre-B Small Pre-B Trans B Naive B Memory B Plasma
  1     0           0        1059     513       0        0      0
  2     0           0           0       0    1014        0      0
  3     0           0           0       0       0      809      0
  4     0           0           0       0       0        0    554
  5     0         539           0       0       0        0      0
  6   371           0           0       0       0        0      0

标注已经定在 annot_subtype 上。丢掉刚才那批细胞后再 Normalize → HVG → PCA → Harmony → UMAP（dims = 1:20，reduction 名 umap），不再聚类。

Idents(seurat_object) <- "sample"
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]
)
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.057 sec elapsed

seurat_object <- RunUMAP(
  seurat_object, dims = 1:20, reduction = "harmony",
  reduction.name = "umap"
)

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_b_plasma.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