Subtype: Erythrocyte

红系亚型。珠蛋白主导高变基因是预期现象。

Level-1 Ery（含偏成熟的 cluster 36）。HVG 会被珠蛋白主导，属预期。同时看 MK（PF4/PPBP）和其它谱系污染。

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/erythrocyte"
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("Ery"))
dim(seurat_object)

[1] 38606  6456

table(seurat_object$annot)

 Ery 
6456 

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 2723    0    0    0    0    0    0    0 
  17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 
1746    0    0    0    0 1466    0    0    0    0    0    0    0    0    0    0 
  33   34   35   36   37   38 
   0    0    0  521    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] "HBG2"     "IGHV4-34" "IGKV3-15" "IGLV3-19" "IGKV1-9"  "IGKV3-20"
 [7] "HBG1"     "S100A8"   "CAMP"     "IGHV4-39" "LCN2"     "DEFA1"   
[13] "S100A12"  "IGLV4-69" "IGHV2-5"  "DEFA3"    "LTF"      "GNLY"    
[19] "LYZ"      "XACT"    

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

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: 3.137 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] 22

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(
  MEP = c("GATA2", "KIT", "CD34", "ITGA2B"),
  early = c("GATA1", "KLF1", "TFRC", "CD36", "EPOR", "MYB"),
  intermediate = c("CA1", "HBD", "AHSP"),
  late = c("GYPA", "HBB", "HBA1", "ALAS2", "SLC4A1"),
  MK = c("GP9", "PPBP", "PF4"),
  cycling = c("MKI67", "TOP2A"),
  HSPC = c("CD34", "SPINK2", "AVP"),
  gran_mono = c("MPO", "ELANE", "CAMP", "FCGR3B", "CD14", "FCN1"),
  lymphoid = c("CD3D", "NKG7", "MS4A1", "CD79A", "JCHAIN")
)

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 
356 317 309 308 293 284 259 249 237 229 228 227 215 209 206 194 188 182 181 174 
 21  22  23  24  25  26  27  28  29  30  31  32  33 
172 171 165 147 132 132 131 131 123  93  92  73  49 

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

marker_table <-
  erythrocyte_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 FAM1… RRM2  CENPA PTTG1 H4C3  PTH2R STAG1 MAN1… TMCC2 IL32  JCHA… CDC20
 2 Top_2 APOC1 MYBL2 UBE2C ENSG… H3C3  FIRRE RABG… SLC2… IFIT… CD3E  IGLV… CCNB1
 3 Top_3 NME1  CARH… UBE2S ARL6… H3C2  XACT  DLEU2 SOX6  ARL4A PTPR… IGKV… HMMR 
 4 Top_4 AKR1… H2AC… PLK1  UBE2S CDKN… CRACD SLC2… GPCP… XPO7  CD52  IGHV… CCNB2
 5 Top_5 SRM   CCNE2 KPNA2 CENPF TUBA… CDK6  JAZF1 TRAK2 TCP1… CCL5  IGHV… PLK1 
 6 Top_6 EPCAM TUBB  ASPM  HMMR  H3C4  LINC… SPTB  XPO7  TRIM… CD48  IGHV… DEPD…
 7 Top_7 NMU   PKMY… KIF14 CDC2… H1-2  RYR3  PBX3  TOP1  UBE2H FYN   IGLL5 GAS2…
 8 Top_8 ALDH… SAMD1 CDKN… ASPM  NUSA… THRB  ANK1  ABCC… MKRN1 CD247 IGHM  PTTG1
 9 Top_9 PEBP1 H3C3  CDK1  CCNB1 H1-5  CCDC… NSD2  SLC4… BNIP… ETS1  IGKV… CCND3
10 Top_… MYC   TYMS  CENPF CKAP2 H2AC… MYB   SOX6  SPTA1 PIM1  IL7R  IGKC  PRDX2
11 Top_… PPP1… BTG3  CENPE PTMS  H2AC… SLC3… FAF1  RAPG… TBCEL CD3D  CD74  NEK2 
12 Top_… LDHB  DYNC… ARL6… NEK2  TUBB  ENSG… CENPP RHAG  UBR2  PLAA… IGKV… LMNA 
13 Top_… SYNG… E2F8  TUBA… GPSM2 HMGB2 CCSE… DIAP… ENSG… LPIN2 GIMA… MZB1  CENPA
14 Top_… GCSH  RBBP8 HMGB2 GYPE  H2AC… LDLR… NCAP… TSPA… MAN1… PRKCH IGHA1 ARL6…
15 Top_… GIHCG H1-3  GTSE1 HEMGN H3C10 FAM1… MYO1D JAZF1 BPGM  TRAF… TMSB… UROS 
16 Top_… NME2  FEN1  KIF1… CDKN3 HBM   HTR1F CIT   TBCEL BCL2… PCED… TMSB… CDKN3
17 Top_… GLRX  H3C7  CDKN3 GYPA  SAMD1 PVT1  MKI67 GYPB  ATG14 CD2   B2M   TUBB…
18 Top_… PVT1  SLBP  CCNB1 GYPB  STAG1 RANB… EZH2  BNIP… DCAF… GZMA  IGKV… LXN  
19 Top_… CNRI… TUBA… PTTG1 TROAP MKI67 RREB1 CR1L  HECT… RBM38 LCK   HBG1  BIRC5
20 Top_… MIF   TMEM… CDCA2 HBM   H2AC… ABO   NFIA  SPEC… BSG   CD96  CYBA  MYL4 
# ℹ 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(10, 11, 18, 20, 26, 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(15)
# 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.485 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
)

7. Markers for annotation

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

人骨髓红系是 MEP → Early Ery → Intermediate Ery → Late Ery 连续轴。对照下面四组 marker：MEP 看 GATA2/KIT/CD34（珠蛋白应低）；Early 看 TFRC/CD36/EPOR/GATA1/KLF1；Intermediate 看 CA1/HBD/AHSP；Late 看 HBB/GYPA/ALAS2/SLC4A1（TFRC/MKI67 下降、UMI 低）。

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

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("TFRC", "CD36", "KIT", "MYB", "EPOR", "MKI67"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CA1", "HBD", "AHSP", "KLF1", "GATA1"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("GYPA", "HBB", "HBA1", "HBA2", "ALAS2", "SLC4A1"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("GATA2", "KIT", "TFRC", "CD36", "CA1", "AHSP", "GYPA", "HBB", "ALAS2"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

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

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

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

ggsave(file.path(res_path, "05_markers_MK.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 = 12, height = 5)

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c(
    "CD34", "SPINK2", "AVP",
    "MPO", "ELANE", "CAMP", "FCGR3B", "CD14", "FCN1",
    "CD3D", "NKG7", "MS4A1", "CD79A", "JCHAIN",
    "CPA3", "TPSAB1"
  ),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

8. Annotation

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

命名用的分辨率须与 cluster 对照图一致。按 MEP / Early Ery / Intermediate Ery / Late Ery 填。

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[c("MEP", "early", "intermediate", "late")],
  legend.length = 10
)

ggsave(
  file.path(res_path, "05_markers_stages.pdf"),
  width = 14, height = 6, dpi = 600
)

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 
1686 1586 1436  505 

erythrocyte_markers <- FindAllMarkers(seurat_object, only.pos = TRUE)
marker_table <-
  erythrocyte_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 × 5
   rank   `1`     `2`    `3`      `4`     
   <chr>  <chr>   <chr>  <chr>    <chr>   
 1 Top_1  S100A4  NUSAP1 XPO7     HBB     
 2 Top_2  SNHG29  MKI67  IFIT1B   UBB     
 3 Top_3  S100A6  UBE2S  ARL4A    DCAF12  
 4 Top_4  NPM1    CENPF  TMCC2    MKRN1   
 5 Top_5  RPL23   DLEU2  UBE2H    R3HDM4  
 6 Top_6  ATP5MC2 LBR    MAN1A1   ADIPOR1 
 7 Top_7  DUT     KNL1   TRIM58   BNIP3L  
 8 Top_8  NME2    SMC4   BNIP3L   TRIM58  
 9 Top_9  HSPD1   HMGB2  SLC25A37 ALAS2   
10 Top_10 HNRNPA1 TUBA1B GPCPD1   BCL2L1  
11 Top_11 EEF1B2  H1-0   TRAK2    NCOA4   
12 Top_12 SYNGR1  HEMGN  TBCEL    BPGM    
13 Top_13 NDUFS5  HBM    HECTD4   IFIT1B  
14 Top_14 RPL10A  CDK1   TOP1     OAZ1    
15 Top_15 GAPDH   TOP2A  TSPAN5   RIOK3   
16 Top_16 HINT1   TUBB4B HBD      RBM38   
17 Top_17 LDHB    H4C3   ABCC13   SLC25A39
18 Top_18 NASP    ANP32E BSG      TENT5C  
19 Top_19 ATP5MC1 CDKN2D RAPGEF2  FKBP8   
20 Top_20 EEF1G   TMPO   ALAS2    RNF10   

cluster_map <- c(
  "1" = "Early Ery",
  "2" = "Early Ery",
  "3" = "Late Ery",
  "4" = "Late Ery"
)

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("Early Ery", "Late Ery")
  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")
}

Early Ery  Late Ery 
     3272      1941 

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