Subtype: T / NK

T / NK 亚型。PCA 不含 TCR 基因，避免克隆型主导分群。

Level-1 T/NK。Scale / PCA / clustering 排除 TR 基因，避免 TCR 克隆型主导分群。cluster 20 的混合在 subtype 里拆开。CellTypist 只作参考。

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/tnk"
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("T/NK"))
dim(seurat_object)

[1] 38606 15792

table(seurat_object$annot)

 T/NK 
15792 

table(seurat_object$harmony_clusters_2)

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

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

[1] "counts"

Layers(seurat_object)

[1] "counts"

1. Normalize, HVG, Scale, PCA

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

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

 [1] "CAMP"            "S100A8"          "S100A12"         "HBB"            
 [5] "CCL4L2"          "LTF"             "HBA2"            "PTGDS"          
 [9] "HBA1"            "LCN2"            "S100A9"          "DEFA3"          
[13] "DEFA1"           "CCL4"            "TRBV7-2"         "ENSG00000227240"
[17] "MT1G"            "MT1E"            "TRBV10-3"        "HSPA6"          

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

[1] 35

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

num_pcs$pcs

[1] 12

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: 8.55 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] 24

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(
  T_lineage = c("CD3D", "CD3E", "CD3G", "CD4", "CD8A", "CD8B"),
  naive = c("CCR7", "SELL", "LEF1", "TCF7", "S1PR1"),
  Tcm = c("IL7R", "ANXA1", "ANXA2", "AQP3", "GPR183", "NELL2", "PTGER2"),
  Tem = c("GZMK", "CCL5", "EOMES", "KLRG1"),
  CTL = c("GZMB", "PRF1", "GZMH", "CX3CR1"),
  NK = c("KLRF1", "NCAM1", "NCR1", "TYROBP", "FCER1G", "KLRD1"),
  NKT = c("TRAV10"),
  Treg = c("FOXP3", "IL2RA", "CTLA4"),
  gdT = c("TRDC", "TRGC1", "TRGC2", "TRGV9", "TRDV1", "TRDV2", "TRDV3"),
  MAIT = c("KLRB1", "SLC4A10", "TRAV1-2", "ZBTB16"),
  cycling = c("MKI67", "TOP2A"),
  contaminants = c(
    "MS4A1", "CD79A", "CD14", "CD68", "CD34", "MPO",
    "AZU1", "NAMPT", "CSF3R", "S100A8", "S100A9", "CA1", "HBB", "PF4", "PPBP"
  )
)

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 
920 780 775 742 714 711 702 688 607 585 531 530 527 518 504 472 468 458 458 401 
 21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36 
378 343 342 335 313 295 280 262 254 234 185 178 106  90  54  52 

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

marker_table <-
  tnk_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 × 37
   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 FCGR… FGFB… GZMK  GZMK  NELL2 CCR7  MNDA  CD40… FCER… RPL32 GZMK  GZMH 
 2 Top_2 GZMB  ENSG… HLA-… IFNG… LEF1  TSHZ2 NAMPT LINC… KLRF1 RPS12 RPL41 ZNF6…
 3 Top_3 KLRC2 ZNF6… CD74  CD27  CCR7  ANKR… CSF3R IL7R  FCGR… RPL30 RPS27 FGFB…
 4 Top_4 FGFB… GZMH  HLA-… CD8B  BACH2 RPL32 FCGR… CD4   CLIC3 RPL13 CD8A  CD52 
 5 Top_5 GZMH  PRSS… HLA-… MS4A1 LINC… LEF1  FPR1  LTB   SPON2 RPS13 CCL5  B2M  
 6 Top_6 GNLY  GZMB  CCL5  CD74  ACTN1 AK5   ACSL1 ENSG… IGFB… RPL11 CD74  NKG7 
 7 Top_7 KLRC3 ZEB2  HLA-… CRTAM CD248 RPS13 S100… GPR1… SH2D… RPL10 B2M   SH3B…
 8 Top_8 KLRF1 ITGB1 HLA-… CD8A  PTPRK FHIT  IL1R2 VIM   NCR1  EEF1… IL32  GZMB 
 9 Top_9 NKG7  CD8A  CD8A  ST8S… LRRN3 RPL30 RGS2  TIMP1 CHST2 RPS1… CD8B  HLA-C
10 Top_… TRDC  LINC… HLA-… DTHD1 ENSG… RPL11 G0S2  INPP… MYOM2 RPS2… RPS18 CCL5 
11 Top_… KLRD1 LILR… CXCR6 CXCR6 NT5E  RPS8  BASP1 ANK3  AKR1… RPS8  RPS3  S100…
12 Top_… CST7  IL7   IL32  TOX   TRAB… RPS12 CMTM2 TNFR… LING… RPS28 RPS1… HLA-B
13 Top_… PRF1  NKG7  CD8B  DUSP2 CD8B  MAL   LST1  TPT1  L3MB… RPL39 RPL28 CD8A 
14 Top_… ENSG… CX3C… GZMA  CCL5  CD79A RPS23 TREM1 IL2   DAB2  RPS23 EEF1… CX3C…
15 Top_… CX3C… MTSS1 HLA-… RGS1  AIF1  RPS14 RBM47 ADAM… PLXN… RPS5  RPL39 EFHD2
16 Top_… TYRO… ADGR… FABP5 THEM… ARMH1 RPL34 SPI1  AQP3  ZMAT4 RPL34 RPS12 TRGV4
17 Top_… ADGR… CD8B  HLA-… SH2D… REG4  RPL10 AQP9  LINC… LINC… RPL1… RPL13 CD3D 
18 Top_… ENSG… C1or… PDCD1 CMC1  EPHA… RPL9  LRRK2 PLP2  GZMB  RPL19 RPS28 GNLY 
19 Top_… PRSS… CD52  CCL4… CD3E  RPS13 RPL19 S100… CD52  CD38  RPL28 RPS25 CD99 
20 Top_… KIR2… MYO3B HLA-A IL32  RPL32 RPL13 KCNJ… DPP4  GSAP  RPS4X HLA-B HLA-…
# ℹ 24 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>, `34` <chr>,
#   `35` <chr>, `36` <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(7, 14, 17, 31, 32, 33, 34, 35)
# 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$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(27)
# 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$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: 6.873 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
GZMK      0.000000e+00  2.4497524 0.952 0.266  0.000000e+00
HLA-DRB1  0.000000e+00  1.7047339 0.802 0.386  0.000000e+00
HLA-DQA1  0.000000e+00  2.6304833 0.424 0.089  0.000000e+00
CD74      0.000000e+00  1.9465836 0.966 0.724  0.000000e+00
CCL5     1.044007e-307  1.0687538 0.986 0.668 4.030492e-303
CD8A     7.278640e-279  1.3892160 0.769 0.347 2.809992e-274
HLA-DRA  5.398178e-271  2.2209343 0.501 0.157 2.084020e-266
CD8B     4.056125e-255  1.3479243 0.723 0.304 1.565907e-250
CXCR6    1.072029e-249  2.4453465 0.338 0.073 4.138673e-245
RGS1     1.281383e-237  2.7427709 0.355 0.089 4.946906e-233
HLA-DPA1 1.798370e-217  1.3826376 0.739 0.424 6.942788e-213
GZMA     5.835397e-213  0.9394927 0.935 0.621 2.252813e-208
CCL4     4.315507e-202  1.3759568 0.786 0.445 1.666045e-197
HLA-DQB1 1.279999e-197  2.0160291 0.436 0.150 4.941566e-193
IL32     3.261859e-194  0.7535600 0.985 0.881 1.259273e-189
TOX      1.641898e-184  1.3542258 0.685 0.372 6.338713e-180
CCL4L2   1.708478e-171  1.9419823 0.451 0.177 6.595751e-167
HLA-DPB1 1.211048e-162  1.1243969 0.752 0.481 4.675372e-158
CRTAM    7.251603e-162  1.8135436 0.366 0.122 2.799554e-157
HLA-DRB5 2.983583e-149  1.9783226 0.325 0.104 1.151842e-144

7. Markers for annotation

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

先分 lineage，再分状态。Scale/PCA 已排除 TR；TRAV1-2 / TRAV10 只用来看，不进 PCA。

Naive / Tcm / Tem / CTL（先看 CD4 vs CD8，再套这套）：

1. Naive：LEF1 TCF7 CCR7 SELL 高，颗粒酶全无。
2. Tcm：淋巴归巢还在（FeaturePlot 对照 CCR7/SELL），LEF1 不如 naive 齐，有 IL7R ANXA1 ANXA2 AQP3 GPR183；CD8 再看 NELL2。没有 GZMK/GZMB。和 naive 糊在一起就并成 naive/Tcm。
3. Tem：归巢掉，GZMK+ CCL5+，GZMB 还不高。
4. CTL：GZMB PRF1 GZMH CX3CR1 高，GZMK 往往下来。CD8 最清楚；CD4 CTL 少但分子一样。

NK vs NKT / NKT-like：

• NK：CD3-（CD3D/E/G 基本无），KLRF1 NCAM1 NCR1 TYROBP FCER1G 高。
• iNKT（真 NKT）：CD3+，TRAV10（Vα24），可有 ZBTB16；人骨髓里很少。和 MAIT 的差别是 MAIT 用 TRAV1-2/SLC4A10，不是 TRAV10。
• CD8 NKT-like（旧稿里的 Tem/NK-like、CTL/NK-like）：CD3+ CD8+，又带 NK 受体，但 KLRF1/NCAM1/TYROBP 不如真 NK。这是常规 CD8，不是 iNKT；不要只靠 NKG7/GNLY/FGFBP2（CTL 也有）。

DotPlot 的 marker_groups 每个基因只出现一次。FeaturePlot 仍可对照共用分子（如 Tcm 图里的 CCR7）。

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CD3D", "CD3E", "CD3G", "CD4", "CD8A", "CD8B"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CCR7", "SELL", "LEF1", "TCF7", "S1PR1"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c(
    "CCR7", "SELL", "IL7R", "ANXA1", "ANXA2",
    "AQP3", "GPR183", "NELL2", "PTGER2"
  ),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("GZMK", "CCL5", "EOMES", "KLRG1", "IL7R"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("GZMB", "PRF1", "GNLY", "NKG7", "GZMH", "CX3CR1", "FGFBP2"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("KLRF1", "NCAM1", "NCR1", "TYROBP", "FCER1G", "KLRD1"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("CD3D", "TRAV10", "KLRF1", "NCAM1", "TYROBP", "FCER1G"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

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

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("TRDC", "TRGC1", "TRGC2", "TRGV9", "TRDV1", "TRDV2", "TRDV3"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 2,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("KLRB1", "SLC4A10", "TRAV1-2", "ZBTB16"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_MAIT.pdf"), width = 12, 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("MS4A1", "CD79A", "CD14", "CD68", "CD34", "MPO", "CA1", "HBB"),
  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 对照图一致。先填标签，再看 cluster 7 的 QC，最后决定 7 要不要进 drop_pre_annot。14 / 15 确定丢掉。丢掉后按 B/plasma：Normalize → HVG → PCA → Harmony → UMAP（dims = 1:20），不再聚类。

annot_cluster <- "harmony_clusters_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    6    7    8    9   10   11   12   13   14   15 
1804 1741 1503 1492 1458 1450  944  778  543  465  396  316  255  126   86 

tnk_markers <- FindAllMarkers(seurat_object, only.pos = TRUE)
marker_table <-
  tnk_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 × 16
   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 GZMH  GZMK  LTB   FCER… FCGR… LEF1  GZMK  NELL2 IFNG… SPRY2 SLC4… RTKN2
 2 Top_2 FGFB… HLA-… IL7R  TYRO… GZMB  LTB   RPL41 LEF1  IL7R  IKZF2 RORC  FOXP3
 3 Top_3 GZMB  HLA-… CD40… KLRF1 GNLY  CCR7  RPS27 CCR7  GZMK  TRDV3 ME1   RBMS3
 4 Top_4 ZNF6… CD74  ADAM… FCGR… KLRC2 MAL   RPL39 BACH2 MS4A1 ITGAD IL23R IL2RA
 5 Top_5 NKG7  CCL5  CD4   CLIC3 FGFB… TSHZ2 RPS28 LINC… CD8A  GCSAM LTK   ENSG…
 6 Top_6 CD52  CD8A  CCR6  GZMB  GZMH  ANKR… RPS12 ACTN1 CD8B  TRDV1 ADAM… CTLA4
 7 Top_7 ZEB2  HLA-… PTPN… FGFB… KLRD1 RPL22 RPS1… AIF1  FXYD2 CMC1  LINC… F5   
 8 Top_8 ENSG… CD8B  ANK3  SPON2 KLRC3 RPS5  RPS2… CD248 PARP8 ENSG… ENSG… DNAH8
 9 Top_9 ITGB1 CXCR6 GPR1… GNLY  TYRO… RPL5  RPL30 PTPRK NELL2 XCL1  P2RY… IKZF2
10 Top_… SH3B… RGS1  AQP3  KLRD1 KLRF1 RPL9  RPL28 LRRN3 CAMK4 XCL2  CCR6  CDK14
11 Top_… PRSS… HLA-… VIM   SH2D… ZEB2  RPS10 RPL32 CD8B  CASK  DRAIC IL7R  RBMS…
12 Top_… GNLY  GZMA  INPP… CEBPD CX3C… RPL29 RPS3  RPS13 PCAT1 ENSG… DPP4  FANK1
13 Top_… EFHD2 CCL4  TIMP1 KLRB1 PRF1  RPL18 RPL13 CD79A CXCR3 CD160 CEBPD ICA1 
14 Top_… S100… HLA-… S100… CX3C… C1or… RPL1… IL32  RPL32 THEM… AOAH  TRBV… ENSG…
15 Top_… CX3C… IL32  TNFR… PTPN… PRSS… RPL3… RPS25 RPS12 LSR   CD74  TRAV… NEBL 
16 Top_… CD8A  TOX   ENSG… IGFB… GPR1… RPS6  RPL34 TRAB… DUSP2 GZMK  PRSS… ENSG…
17 Top_… C1or… CCL4… SERI… MCTP2 ADGR… RPS25 RPS18 ENSG… MT-N… CCL5  KIF5C TBC1…
18 Top_… HLA-C HLA-… TPT1  PRF1  ENSG… RPS13 RPS15 RPS8  CRYB… HLA-… B3GA… HACD1
19 Top_… HLA-… CRTAM LINC… TTC38 TRDC  RPS23 RPS23 NT5E  NR3C2 CRTAM RUNX2 PLCL1
20 Top_… LINC… HLA-… ADAM… GPR1… ENSG… RPS14 B2M   RPS23 BTG1  EOMES KLRB1 LINC…
# ℹ 3 more variables: `13` <chr>, `14` <chr>, `15` <chr>

cluster_map 的 key 必须对应 annot_cluster。

cluster_map <- c(
  "1" = "CD8 CTL",
  "2" = "CD8 Tem",
  "3" = "CD4 Tcm",
  "4" = "CD56dim NK",
  "5" = "CD8 NKT",
  "6" = "CD4 Tnaive",
  "7" = "Ribo",
  "8" = "CD8 Tnaive",
  "9" = "CD8 Tem",
  "10" = "gdT",
  "11" = "MAIT",
  "12" = "CD4 Treg",
  "13" = "CD56bright NK"
)

ribo_cluster <- "7"
cluster_ids <- as.character(FetchData(seurat_object, vars = annot_cluster)[, 1])
is_ribo <- cluster_ids == ribo_cluster
message(sprintf("Cluster %s: %s cells.", ribo_cluster, sum(is_ribo)))

qc_vars <- c(
  annot_cluster, "nFeature_RNA", "nCount_RNA", "percent_mt",
  "is_doublet", "sample"
)
if ("doublet_score" %in% colnames(seurat_object[[]])) {
  qc_vars <- c(qc_vars, "doublet_score")
}
ribo_qc <- FetchData(seurat_object, vars = qc_vars)[is_ribo, ]

print(summary(ribo_qc$nFeature_RNA))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  303.0   622.2   860.0   949.0  1178.5  2871.0 

print(summary(ribo_qc$nCount_RNA))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    501    1059    1634    2028    2605   12716 

print(summary(ribo_qc$percent_mt))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8667  1.4867  1.6509  2.2590  4.8458 

print(table(ribo_qc$is_doublet, useNA = "ifany"))

singlet doublet 
    938       6 

print(table(ribo_qc$sample, useNA = "ifany"))

 ACI-1  ACI-2  ACI-3  ACI-4  ACI-5 CTRL-1 CTRL-2 CTRL-3 CTRL-4 CTRL-5 CTRL-6 
    30     25     44     15     17    103    330    103    201     14     62 

print(round(prop.table(table(ribo_qc$sample)), 3))

 ACI-1  ACI-2  ACI-3  ACI-4  ACI-5 CTRL-1 CTRL-2 CTRL-3 CTRL-4 CTRL-5 CTRL-6 
 0.032  0.026  0.047  0.016  0.018  0.109  0.350  0.109  0.213  0.015  0.066 

seurat_object$ribo_check <- ifelse(is_ribo, paste0("cluster_", ribo_cluster), "other")
VlnPlot(
  seurat_object,
  features = intersect(
    c("nFeature_RNA", "nCount_RNA", "percent_mt", "doublet_score"),
    colnames(seurat_object[[]])
  ),
  group.by = "ribo_check",
  pt.size = 0, ncol = 2
)

ggsave(file.path(res_path, "05_markers_ribo_qc.pdf"), width = 10, height = 8)

DimPlot(
  seurat_object, reduction = "umap", group.by = "sample",
  cells.highlight = colnames(seurat_object)[is_ribo],
  cols.highlight = "red", sizes.highlight = pt_size,
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
) + theme(aspect.ratio = 1)

ggsave(file.path(res_path, "05_markers_ribo_umap.pdf"), width = 8, height = 8)

# 14 / 15 确定删。看完上面 cluster 7 的 QC 后，若要删 7，改成 c("14", "15", "7")
drop_pre_annot <- c("7", "14", "15")
n_before <- ncol(seurat_object)
keep_cells <- colnames(seurat_object)[
  !(as.character(FetchData(seurat_object, vars = annot_cluster)[, 1]) %in%
    as.character(drop_pre_annot))
]
seurat_object <- subset(seurat_object, cells = keep_cells)
message(sprintf(
  "Dropped clusters %s: %s cells; %s remain.",
  paste(drop_pre_annot, collapse = ", "),
  n_before - ncol(seurat_object),
  ncol(seurat_object)
))

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(
    "CD4 Tnaive", "CD4 Tcm", "CD4 Treg",
    "CD8 Tnaive", "CD8 Tem", "CD8 CTL", "CD8 NKT",
    "MAIT", "gdT", "CD56bright NK", "CD56dim NK"
  )
  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")
}

   CD4 Tnaive       CD4 Tcm      CD4 Treg    CD8 Tnaive       CD8 Tem 
         1450          1503           316           778          2284 
      CD8 CTL       CD8 NKT          MAIT           gdT CD56bright NK 
         1804          1458           396           465           255 
   CD56dim NK 
         1492 

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

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$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: 5.558 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_tnk.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