Subtype: Neutrophil

中性粒亚型：早幼粒、Pre-Neut、成熟粒。后续组成和状态分析用这些名字。

第一轮只 subset Neut + Pre-Neut，不自动并入 HSPC cluster 21。GMP / promyelocyte barcode 确认后再另存。低质量 / doublet 留给 round 2。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/neutrophil"
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("Neut", "Pre-Neut"))
dim(seurat_object)

[1] 38606 34800

table(seurat_object$annot)

    Neut Pre-Neut 
   20978    13822 

table(seurat_object$harmony_clusters_2)

   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
   0 4060 3900    0 3382 3210    0 3131    0    0 2664 2557 2501 2318    0    0 
  17   18   19   20   21   22   23   24   25   26   27   28   29   30   31   32 
   0 1583 1551    0    0    0    0    0 1266 1262    0    0    0  989    0    0 
  33   34   35   36   37   38 
   0    0    0    0  426    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] "DEFA1"  "IGKC"   "DEFA3"  "H4C3"   "HBB"    "IGLC3"  "DEFA4"  "IGLC2" 
 [9] "HSPA1B" "HBG2"   "CCL4"   "JCHAIN" "IGHA1"  "MPO"    "CCL4L2" "OLFM4" 
[17] "IGHG1"  "CAMP"   "PRG2"   "CLC"   

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

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

num_pcs$co2

[1] 14

num_pcs$pcs

[1] 14

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: 16.562 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] 44

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(
  promyelocyte = c("AZU1", "MPO", "ELANE", "CTSG", "PRTN3", "CD34"),
  pre_neut = c("LTF", "LCN2", "CAMP", "MNDA", "RETN"),
  mature_neut = c("FCGR3B", "CXCR2", "NAMPT", "CSF3R", "S100A8", "S100A9"),
  cycling = c("MKI67", "TOP2A", "PCNA"),
  HSPC = c("CD34", "SPINK2", "AVP"),
  mono = c("CD14", "FCN1", "CD68"),
  lymphoid = c("CD3D", "MS4A1", "CA1", "HBB"),
  mast = c("CPA3", "TPSAB1")
)

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

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

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

   1   10   11   12   13   14   15   16   17   18   19    2   20   21   22   23 
1990 1266 1222 1189 1089 1062  994  990  954  951  926 1842  865  825  741  719 
  24   25   26   27   28   29    3   30   31   32   33   34   35   36   37    4 
 714  654  646  619  607  579 1736  443  406  397  389  301  267  231  208 1693 
   5    6    7    8    9 
1499 1484 1484 1459 1359 

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

marker_table <-
  neutrophil_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 × 38
   rank  `1`   `10`  `11`  `12`  `13`  `14`  `15`  `16`  `17`  `18`  `19`  `2`  
   <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
 1 Top_1 MMP9  NAMPT NAMPT S100… SLC2… NAMPT DEFA1 MMP8  DEFA3 CXCL8 LTF   ENSG…
 2 Top_2 TENM1 NEAT1 DUSP1 ACSL1 IFIT… G0S2  DEFA3 CRIS… DEFA1 NFKB… ABCA… CXCR2
 3 Top_3 S100… LUCA… S100… IFIT… RBM47 RUBC… DEFA4 S100… RPL10 TNFA… LCN2  TXNIP
 4 Top_4 S100… MALA… IFIT… S100… NIBA… ACSL1 SLPI  S100… RPS12 ICAM1 MMP8  LUCA…
 5 Top_5 S100… IFIT… FOS   CMTM2 NEAT1 SRGN  MS4A3 LCN2  ENSG… LINC… RETN  CAMK…
 6 Top_6 LGAL… TREM1 RGS2  RIPO… ACSL4 RGS2  ABCA… PGLY… RPL32 CD83  CEAC… AIF1 
 7 Top_7 HMGB2 SLC2… CMTM2 S100… SLC8… PROK2 PDE4D LTF   RPS27 RELB  CAMP  SEC1…
 8 Top_8 MIR6… FOS   SRGN  MME   CSF3R IFIT… ATP8… MMP9  RPLP1 NEAT1 S100… SLC2…
 9 Top_9 PADI4 KCNJ… FCGR… RGS2  RIPO… CR1   BPI   HP    RPL28 SOD2  LYZ   NIBA…
10 Top_… SH3R… SRGN  CSF3R FCGR… LYN   BCL2… PTMA  CYBB  RPL41 CAMK… CYBB  CPPE…
11 Top_… S100P RGS2  ACSL1 NAMPT ARHG… OSM   AZU1  PADI4 RPS23 SLC2… CRIS… ADGR…
12 Top_… RBP7  MME   SAT1  H3-3A VMP1  AQP9  LRMDA ENSG… RPS14 MALA… CPNE3 VMP1 
13 Top_… SYNE1 FTH1  MME   CR1   ITM2B PTGS2 RPL22 S100… RPL1… FTH1  LTA4H CLEC…
14 Top_… ALOX… FCGR… S100… CSF3R CAMK… MCTP2 GSTP1 ADSS2 RPS29 NFKB… RPL39 SCLT1
15 Top_… IL1R2 CSF3R NEAT1 FKBP5 FTH1  FOS   CD24  SYNE1 RPL13 NIBA… PGLY… SMAP2
16 Top_… SAMS… SAT1  FTH1  VNN2  CXCR2 KCNJ… CEAC… RETN  RPS21 NAMPT OLR1  ENSG…
17 Top_… TXN   DUSP1 H3-3A XPO6  SMAP2 PLAUR CEAC… SVIP  RPS18 LYN   ADSS2 LITAF
18 Top_… RGL4  G0S2  RIPO… SRGN  KCNJ… S100… RNAS… PLBD1 RPL19 ENSG… CD24  RIPO…
19 Top_… HIPK2 RIPO… S100… LINC… ALPL  MME   NUCB2 RPL39 RPS28 PDE4B ANXA1 PPCDC
20 Top_… STK1… CMTM2 FPR1  SELL  LITAF TREM1 ELANE STOM  RPL11 SAT1  RPS8  NEAT1
# ℹ 25 more variables: `20` <chr>, `21` <chr>, `22` <chr>, `23` <chr>,
#   `24` <chr>, `25` <chr>, `26` <chr>, `27` <chr>, `28` <chr>, `29` <chr>,
#   `3` <chr>, `30` <chr>, `31` <chr>, `32` <chr>, `33` <chr>, `34` <chr>,
#   `35` <chr>, `36` <chr>, `37` <chr>, `4` <chr>, `5` <chr>, `6` <chr>,
#   `7` <chr>, `8` <chr>, `9` <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(17, 30, 35, 36)
# drop_clusters <- c("3", "7")

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

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

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

5b. Check cleanup

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

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

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

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

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

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

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

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

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

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

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

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

6. Recluster for annotation

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

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

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

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

toc()

Harmony integration: 16.174 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
MME        0  1.5195838 0.619 0.247         0
CMTM2      0  1.1172016 0.764 0.401         0
RGS2       0  0.8813532 0.866 0.517         0
AQP9       0  0.7338901 0.691 0.375         0
FCGR3B     0  0.6716172 0.809 0.494         0
CR1        0  1.1442153 0.630 0.320         0
KCNJ15     0  0.9230008 0.647 0.365         0
VNN2       0  1.0294321 0.581 0.303         0
RIPOR2     0  0.6215734 0.923 0.646         0
XPO6       0  0.7139534 0.745 0.469         0
PROK2      0  0.8827371 0.579 0.304         0
G0S2       0  0.6653163 0.745 0.476         0
PLXDC2     0  0.8420498 0.622 0.354         0
NAMPT      0  0.9161506 0.990 0.723         0
IL1R2      0  1.0208401 0.841 0.580         0
ACSL1      0  0.9579170 0.931 0.673         0
SELL       0  0.9122509 0.658 0.402         0
LRRK2      0  0.9606403 0.640 0.390         0
S100A4     0  1.0487408 0.968 0.732         0
FTH1       0  0.2727408 0.991 0.769         0

7. Markers for annotation

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

骨髓中性粒按发育轴标，不必硬套肝癌 11 群：

1. Promyelocyte：DEFA1/3/4 AZU1 MPO 高，还没有 LTF/CAMP。
2. Pre-Neut：二次颗粒 LTF LCN2 CAMP RETN 高，MMP9/S100A12 还低。这批细胞多，就是 NP / 带状核前体。
3. Mature：LTF 下来，开始出现 MMP9 SLC8A1 S100A12 TXNIP ISG15。能拆开就按优势基因标成 MMP9 / SLC8A1 / S100A12 / TXNIP / ISG15 Neut；糊在一起就并成 Mature Neut。

LTF+ 又略有 MMP9 的中间态并进 Pre-Neut。Xue 的 TAN 基因（CCL4 SPP1 CD74）骨髓里可以没有。污染看 NKG7/CD3D。

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

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

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

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("LTF", "LCN2", "CAMP", "MNDA", "RETN"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("FCGR3B", "CXCR2", "NAMPT", "CSF3R", "S100A8", "S100A9"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("MMP8", "S100A12", "ISG15", "TXNIP", "ELL2", "PTGS2"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("APOA2", "CD74", "IFIT1", "SPP1", "CCL4", "CD274"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

FeaturePlot(
  seurat_object, reduction = "umap", ncol = 3,
  min.cutoff = min_cutoff, max.cutoff = max_cutoff,
  features = c("MMP9", "SLC8A1", "S100A12", "TXNIP", "ISG15"),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

ggsave(file.path(res_path, "05_markers_states.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", "GNLY", "KLRF1",
    "MS4A1", "CD79A", "CD14", "FCN1",
    "CD34", "CA1", "HBB", "PF4"
  ),
  pt.size = pt_size, alpha = alpha, raster = raster, raster.dpi = raster_dpi
)

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

8. Annotation

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

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

annot_cluster <- "harmony_clusters_0.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 
9013 8301 5934 5740 2102  911  533  371 

neutrophil_markers <- FindAllMarkers(seurat_object, only.pos = TRUE)
marker_table <-
  neutrophil_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  SLC25A37        MMP9   MMP8    MMP9    DEFA1   IFITM3 LTF        CD247
 2 Top_2  LUCAT1          CMTM2  LTF     PGLYRP1 DEFA3   MX1    MMP8       S100…
 3 Top_3  CXCR2           MME    LCN2    S100P   SLPI    IFIT3  CAMP       HLA-…
 4 Top_4  NIBAN1          CR1    CEACAM8 CRISP3  MS4A3   IFIT2  RETN       FTL  
 5 Top_5  CAMK1D          RGS2   RETN    FCN1    PTMA    ISG15  LCN2       DDIT4
 6 Top_6  AQP9            IL1R2  ABCA13  RGL4    ABCA13  IFIT1  ABCA13     HLA-…
 7 Top_7  SLC8A1          NAMPT  CYBB    PADI4   PDE4D   HERC5  S100A9     GNLY 
 8 Top_8  FCGR3B          RIPOR2 CAMP    GAPDH   BPI     IFI6   CYBB       CCL5 
 9 Top_9  RBM47           SELL   ATP8B4  CLC     HMGN2   PARP14 CRISP3     SAP30
10 Top_10 TXNIP           G0S2   CRISP3  TFF3    ATP8B4  XAF1   CEACAM8    RPS3 
11 Top_11 TAGLN2          PLXDC2 PDE4D   TSPO    RPL22   RSAD2  PGLYRP1    CRIP1
12 Top_12 VMP1            LRRK2  ADSS2   TENM1   CD24    EPSTI1 UGCG       CST3 
13 Top_13 TREM1           VNN2   CPNE3   SYNE1   CEACAM8 LY6E   SVIP       CD74 
14 Top_14 RGS2            S100A4 CD24    S100A12 LRMDA   IFI44  ENSG00000… LTB  
15 Top_15 ENSG00000286076 ACSL1  RPL22   PFN1    RNASE2  OAS3   ADSS2      ENSG…
16 Top_16 XPO6            MAP4K4 LTA4H   VIM     DEFA4   CMPK2  LTA4H      CD99 
17 Top_17 CREB5           FPR1   ANXA1   S100A6  RPL36A  DDX60  CEACAM1    FYN  
18 Top_18 ZFP36L1         RBP7   STOM    ACTB    RPS21   TRIM22 LYZ        LIMD2
19 Top_19 SMAP2           SYTL3  PGLYRP1 S100A9  NUCB2   IFITM1 ATP8B4     CARH…
20 Top_20 RIPOR2          DUSP1  DIAPH2  S100A8  H4C3    RIGI   CD24       G0S2 

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

cluster_map <- c(
  "1" = "SLC8A1 Neut",
  "2" = "S100A12 Neut",
  "3" = "MMP8 Pre-Neut",
  "4" = "MMP9 Pre-Neut",
  "5" = "Promyelocyte",
  "6" = "ISG15 Neut",
  "7" = "MMP8 Pre-Neut",
  "8" = "remove"
)

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])
  table(seurat_object$annot_subtype, useNA = "ifany")
}

   ISG15 Neut MMP8 Pre-Neut MMP9 Pre-Neut  Promyelocyte        remove 
          911          6467          5740          2102           371 
 S100A12 Neut   SLC8A1 Neut 
         8301          9013 

DimPlot(
  seurat_object, reduction = "umap", 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, "06_annot_doublets.pdf"),
  width = 8, height = 8, dpi = 600
)

丢掉 remove 和 doublet，再 Normalize → HVG → PCA → Harmony → UMAP（dims = 1:20），不再聚类。

n_before <- ncol(seurat_object)
keep_cells <- colnames(seurat_object)[
  !(seurat_object$annot_subtype %in% "remove" |
    seurat_object$is_doublet %in% "doublet")
]
seurat_object <- subset(seurat_object, cells = keep_cells)
message(sprintf(
  "Dropped remove/doublet: %s cells; %s remain.",
  n_before - ncol(seurat_object),
  ncol(seurat_object)
))
annot_levels <- c(
  "Promyelocyte", "MMP8 Pre-Neut", "MMP9 Pre-Neut",
  "S100A12 Neut", "SLC8A1 Neut", "ISG15 Neut"
)
seurat_object$annot_subtype <- factor(
  seurat_object$annot_subtype,
  levels = intersect(annot_levels, unique(as.character(seurat_object$annot_subtype)))
)

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 | 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: 15.358 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_neutrophil.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