Clustering

高变基因降维、Harmony 整合样本、Leiden 聚类。UMAP 只用来展示，邻近不等于发育先后。

pacman::p_load(Seurat, tidyverse, qs2, harmony, tictoc, cowplot)
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/dim_reduction"
obj_path <- "data/rdata/dim_reduction"
if (!dir.exists(res_path)) dir.create(res_path, recursive = TRUE)
if (!dir.exists(obj_path)) dir.create(obj_path, recursive = TRUE)

alpha <- .1
pt_size <- .1
raster <- FALSE
raster_dpi <- c(3000, 3000)
min_cutoff <- "q10"
max_cutoff <- "q90"

0. Load Gene Info

基因注释（线粒体、核糖体、IG/TR 等），后面排除基因和高变基因会用到。

genes_type <- 
  read_tsv("data/tx_to_genes_cr_2024a.txt") |>
  dplyr::select(gene_symbol, gene_type) |>
  distinct() |>
  arrange(desc(gene_type)) |>
  filter(!duplicated(gene_symbol))

dim(genes_type)

[1] 38584     2

table(genes_type$gene_type)

      IG_C_gene IG_C_pseudogene       IG_D_gene       IG_J_gene IG_J_pseudogene 
             14               9              37              18               3 
      IG_V_gene IG_V_pseudogene          lncRNA  protein_coding       TR_C_gene 
            145             187           18532           19405               6 
      TR_D_gene       TR_J_gene TR_J_pseudogene       TR_V_gene TR_V_pseudogene 
              5              79               4             107              33 

1. data Preparation

接上质控后的 counts，准备降维。

seurat_object <- qs_read("data/rdata/preprocess/data_prep.qs2", nthreads = 10)
DefaultAssay(seurat_object)

[1] "RNA"

Layers(seurat_object)

[1] "counts"

定义要排除的模式

1. MT-: 线粒体

2. RP[SL]: 核糖体

3. HSP: 热休克蛋白 (应激)

4. DNAJ: HSP伴侣 (应激)

5. HIST: 组蛋白 (细胞周期)

6. JUN|FOS|ZFP36: 早期即刻基因 (解离应激)

test_genes <-
  tibble(gene = rownames(seurat_object)) |>
  left_join(
    genes_type |>
      dplyr::rename(gene = gene_symbol)
  ) |>
  mutate(
    unwant = grepl("^MT-|^RP[SL]|^HSP|^DNAJ|^HIST|^ZFP36", gene),
    #  unwant = grepl("^MT-|^RP[SL]", gene),
    ig = grepl("IG", gene_type),
    tr = grepl("TR", gene_type)
  )
sum(test_genes$unwant)

[1] 208

sum(test_genes$ig)

[1] 413

sum(test_genes$tr)

[1] 234

qs_save(test_genes, "data/rdata/test_genes.qs2", nthreads = 10)

2. Without Integration

先看不校正批次时样本在图上散不散。如果各占一块，后面才需要 Harmony。

Layers(seurat_object)

[1] "counts"

Layers(seurat_object) |> length()

[1] 1

# seurat_object <- 
#   split(seurat_object, f = seurat_object$sample)

seurat_object <- NormalizeData(seurat_object, verbose = TRUE)
Layers(seurat_object)

[1] "counts" "data"  

Layers(seurat_object) |> length()

[1] 2

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

 [1] "IGHG1"    "IGHA1"    "IGHG2"    "HBG2"     "IGKC"     "IGHG3"   
 [7] "JCHAIN"   "IGLC2"    "IGKV3-20" "IGHV3-7"  "IGHG4"    "IGKV1-39"
[13] "PTGDS"    "IGHA2"    "HBG1"     "IGLV1-44" "IGLL5"    "CCL4L2"  
[19] "IGKV4-1"  "IGHV4-61"

Layers(seurat_object)

[1] "counts" "data"  

Layers(seurat_object) |> length()

[1] 2

LabelPoints(
  plot = VariableFeaturePlot(seurat_object),
  points = head(VariableFeatures(seurat_object), 20),
  repel = TRUE
)

sum(VariableFeatures(seurat_object) %in% test_genes$gene[test_genes$unwant])

[1] 21

sum(
  VariableFeatures(seurat_object) %in%
    test_genes$gene[test_genes$ig | test_genes$tr]
)

[1] 264

seurat_object <-
  ScaleData(
    seurat_object,
    features = setdiff(
      VariableFeatures(seurat_object),
      test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr]
    ),
    assay = "RNA"
  )

# pca without integration
seurat_object <-
  RunPCA(
    seurat_object,
    npcs = 50,
    features = setdiff(
      VariableFeatures(seurat_object),
      test_genes$gene[test_genes$unwant | test_genes$ig | test_genes$tr]
    ),
    verbose = TRUE,
    reduction.name = "pca"
  )

Number of PCs

https://hbctraining.github.io/scRNA-seq/lessons/elbow_plot_metric.html

https://www.jianshu.com/p/feac6263faee

num_pcs <- list()
num_pcs$pct <-
  seurat_object[["pca"]]@stdev / sum(seurat_object[["pca"]]@stdev) * 100
num_pcs$cumu <- cumsum(num_pcs$pct)

# Determine which PC exhibits cumulative percent greater than 80% and % variation associated with the PC as less than 5
num_pcs$co1 <- which(num_pcs$cumu > 80 & num_pcs$pct < 5)[1]
num_pcs$co1

[1] 32

# last point where change of % of variation is more than 0.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$co2

[1] 19

num_pcs$pcs <- min(num_pcs$co1, num_pcs$co2)
num_pcs$pcs

[1] 19

ElbowPlot(seurat_object, reduction = "pca", ndims = 50)

ElbowPlot(seurat_object, ndims = 50)$data |> ggplot() +
  geom_point(aes(x = dims, y = stdev)) +
  geom_vline(xintercept = num_pcs$pcs, color = "red") +
  theme_bw() + labs(title = "Elbow plot: quantitative approach")

num_pcs$plot_df <- data.frame(pct = num_pcs$pct,
                              cumu = num_pcs$cumu,
                              rank = seq_along(num_pcs$pct))

# Elbow plot to visualize
ggplot(num_pcs$plot_df,
       aes(cumu, pct, label = rank, color = rank > num_pcs$pcs)) +
  geom_text() +
  geom_vline(xintercept = 90, color = "grey") +
  geom_hline(yintercept = min(num_pcs$pct[num_pcs$pct > 5]), color = "grey") +
  theme_bw()

3. Dimentionality Reduction without Integration

不整合时的 PCA / UMAP，用来对照整合效果。

seurat_object <-
  RunUMAP(
    seurat_object,
    dims = 1:30,
    reduction = "pca",
    reduction.name = "umap_unint"
  )

plot_grid(
  ncol = 2,
  DimPlot(
    seurat_object,
    reduction = "umap_unint",
    group.by = "sample",
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object,
    reduction = "umap_unint",
    group.by = "group",
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object,
    reduction = "umap_unint",
    group.by = "phase",
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object,
    reduction = "umap_unint",
    group.by = "is_doublet",
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
)

ggsave("results/dim_reduction/umap_unint_30.pdf",
       width = 14, height = 12, dpi = 600)

4. With integration

Harmony 按样本整合。看整合后各样本是否混开；疾病分组仍可能分开，但其中叠着 5′/3′。

tic("Harmony integration")
seurat_object <-
  RunHarmony(
    seurat_object,
    c("sample"),
    ncores = 15,
    plot_convergence = TRUE
  )

toc()

Harmony integration: 42.554 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$co1

[1] 41

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$co2

[1] 36

umap_temp <-
  RunUMAP(
    seurat_object,
    dims = 1:20,
    reduction = "harmony",
    reduction.name = "umap_harmony"
  )
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(
  "results/dim_reduction/umap_harmony_20.pdf",
  width = 14, height = 12, dpi = 600
)

umap_temp <-
  RunUMAP(
    seurat_object,
    dims = 1:25,
    reduction = "harmony",
    reduction.name = "umap_harmony"
  )
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(
  "results/dim_reduction/umap_harmony_25.pdf",
  width = 14, height = 12, dpi = 600
)

umap_temp <-
  RunUMAP(
    seurat_object,
    dims = 1:30,
    reduction = "harmony",
    reduction.name = "umap_harmony"
  )
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(
  "results/dim_reduction/umap_harmony_30.pdf",
  width = 14, height = 12, dpi = 600
)

umap_temp <-
  RunUMAP(
    seurat_object,
    dims = 1:35,
    reduction = "harmony",
    reduction.name = "umap_harmony"
  )
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(
  "results/dim_reduction/umap_harmony_35.pdf",
  width = 14, height = 12, dpi = 600
)

umap_temp <-
  RunUMAP(
    seurat_object,
    dims = 1:40,
    reduction = "harmony",
    reduction.name = "umap_harmony"
  )
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(
  "results/dim_reduction/umap_harmony_40.pdf",
  width = 14, height = 12, dpi = 600
)

qs_save(
  seurat_object,
  file = "data/rdata/dim_reduction/seurat_object_harmony.qs2",
  nthreads = 10
)
rm(seurat_object, umap_temp)

5. Visualization of Integration

核对整合后的 UMAP：样本应尽量混，cluster 应成团。

seurat_object <-
  qs_read(
    "data/rdata/dim_reduction/seurat_object_harmony.qs2",
    nthreads = 10
  )

n_dims <- 30
set.seed(123)
seurat_object <-
  RunUMAP(
    seurat_object,
    dims = 1:n_dims,
    reduction = "harmony",
    reduction.name = "umap_harmony"
  )

plot_grid(
  ncol = 2,
  DimPlot(
    seurat_object,
    reduction = "umap_harmony",
    group.by = "sample",
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object,
    reduction = "umap_harmony",
    group.by = "group",
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1),
  DimPlot(
    seurat_object,
    reduction = "umap_harmony",
    group.by = "phase",
    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(
  paste0("results/dim_reduction/umap_harmony_", n_dims, ".pdf"),
  width = 14,
  height = 12,
  dpi = 600
)

6. Level 1 Clustering

在整合后的图上聚类。分辨率越大簇越碎。这一步给 04 起大类名用。

seurat_object <-
  FindNeighbors(
    object = seurat_object,
    reduction = "harmony",
    dims = 1:n_dims
  )

for (res in c(0.1, 0.2, 0.5, 1, 1.5, 2)) {
  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

[1] 1.5

[1] 2

gc()

            used   (Mb) gc trigger    (Mb)   max used   (Mb)
Ncells   4034512  215.5    6803639   363.4    6803639  363.4
Vcells 502298089 3832.3 1514474081 11554.6 1207900648 9215.6

gc()

            used   (Mb) gc trigger    (Mb)   max used   (Mb)
Ncells   4034425  215.5    6803639   363.4    6803639  363.4
Vcells 502298159 3832.3 1514474081 11554.6 1207900648 9215.6

qs_save(seurat_object,
        file = "data/rdata/dim_reduction/seurat_object_clustered.qs2",
        nthreads = 10)

plot_grid(
  ncol = 3,
  DimPlot(
    seurat_object,
    reduction = "umap_harmony",
    group.by = "harmony_clusters_0.1",
    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_0.2",
    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_0.5",
    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_1",
    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_1.5",
    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_2",
    pt.size = pt_size,
    alpha = alpha,
    raster = raster,
    raster.dpi = raster_dpi
  ) + theme(aspect.ratio = 1)
)

ggsave("results/dim_reduction/umap_harmony_clusters.pdf",
       width = 21, height = 12, dpi = 600)

DimPlot(seurat_object,
        reduction = "umap_harmony",
        group.by = "harmony_clusters_0.1",
        pt.size = pt_size,
        alpha = alpha,
        raster = raster,
        raster.dpi = raster_dpi,
        label = TRUE) +
  theme(aspect.ratio = 1)

ggsave("results/dim_reduction/umap_harmony_clusters_0.1.pdf",
       width = 8, height = 8, dpi = 600)

DimPlot(seurat_object,
        reduction = "umap_harmony",
        group.by = "harmony_clusters_0.2",
        pt.size = pt_size,
        alpha = alpha,
        raster = raster,
        raster.dpi = raster_dpi,
        label = TRUE) +
  theme(aspect.ratio = 1)

ggsave("results/dim_reduction/umap_harmony_clusters_0.2.pdf",
       width = 8, height = 8, dpi = 600)

DimPlot(seurat_object,
        reduction = "umap_harmony",
        group.by = "harmony_clusters_0.5",
        pt.size = pt_size,
        alpha = alpha,
        raster = raster,
        raster.dpi = raster_dpi,
        label = TRUE) +
  theme(aspect.ratio = 1)

ggsave("results/dim_reduction/umap_harmony_clusters_0.5.pdf",
       width = 8, height = 8, dpi = 600)

DimPlot(seurat_object,
        reduction = "umap_harmony",
        group.by = "harmony_clusters_1",
        pt.size = pt_size,
        alpha = alpha,
        raster = raster,
        raster.dpi = raster_dpi,
        label = TRUE) +
  theme(aspect.ratio = 1)

ggsave("results/dim_reduction/umap_harmony_clusters_1.pdf",
       width = 8, height = 8, dpi = 600)

DimPlot(seurat_object,
        reduction = "umap_harmony",
        group.by = "harmony_clusters_1.5",
        pt.size = pt_size,
        alpha = alpha,
        raster = raster,
        raster.dpi = raster_dpi,
        label = TRUE) +
  theme(aspect.ratio = 1)

ggsave("results/dim_reduction/umap_harmony_clusters_1.5.pdf",
       width = 8, height = 8, dpi = 600)

DimPlot(seurat_object,
        reduction = "umap_harmony",
        group.by = "harmony_clusters_2",
        pt.size = pt_size,
        alpha = alpha,
        raster = raster,
        raster.dpi = raster_dpi,
        label = TRUE) +
  theme(aspect.ratio = 1)

ggsave("results/dim_reduction/umap_harmony_clusters_2.pdf",
       width = 8, height = 8, dpi = 600)

DimPlot(seurat_object,
        reduction = "umap_harmony",
        group.by = "is_doublet",
        pt.size = pt_size,
        alpha = alpha,
        raster = raster,
        raster.dpi = raster_dpi,
        label = TRUE) +
  theme(aspect.ratio = 1) +
  FeaturePlot(seurat_object,
              reduction = "umap_harmony",
              min.cutoff = min_cutoff,
              max.cutoff = max_cutoff,
              features = "doublet_score",
              pt.size = pt_size,
              alpha = alpha,
              raster = raster,
              raster.dpi = raster_dpi) +
  theme(aspect.ratio = 1)

ggsave("results/dim_reduction/umap_harmony_doublet.pdf",
       width = 16, height = 8, dpi = 600)

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] cowplot_1.1.3      tictoc_1.2.1       harmony_1.2.3      Rcpp_1.0.13-1     
 [5] qs2_0.1.7          lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1     
 [9] dplyr_1.1.4        purrr_1.0.2        readr_2.1.5        tidyr_1.3.1       
[13] tibble_3.2.1       ggplot2_3.5.1      tidyverse_2.0.0    Seurat_5.2.1      
[17] SeuratObject_5.0.2 sp_2.1-4          

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3     jsonlite_1.8.9         magrittr_2.0.3        
  [4] spatstat.utils_3.1-1   farver_2.1.2           rmarkdown_2.29        
  [7] ragg_1.3.3             vctrs_0.6.5            ROCR_1.0-11           
 [10] spatstat.explore_3.3-3 htmltools_0.5.8.1      sctransform_0.4.1     
 [13] parallelly_1.39.0      KernSmooth_2.23-26     htmlwidgets_1.6.4     
 [16] ica_1.0-3              plyr_1.8.9             plotly_4.10.4         
 [19] zoo_1.8-12             igraph_2.1.2           mime_0.12             
 [22] lifecycle_1.0.4        pkgconfig_2.0.3        Matrix_1.7-3          
 [25] R6_2.5.1               fastmap_1.2.0          fitdistrplus_1.2-1    
 [28] future_1.34.0          shiny_1.9.1            digest_0.6.37         
 [31] colorspace_2.1-1       patchwork_1.3.0        tensor_1.5            
 [34] RSpectra_0.16-2        irlba_2.3.5.1          textshaping_0.4.0     
 [37] labeling_0.4.3         progressr_0.15.1       fansi_1.0.6           
 [40] spatstat.sparse_3.1-0  timechange_0.3.0       httr_1.4.7            
 [43] polyclip_1.10-7        abind_1.4-8            compiler_4.4.3        
 [46] bit64_4.5.2            withr_3.0.2            fastDummies_1.7.4     
 [49] MASS_7.3-65            tools_4.4.3            lmtest_0.9-40         
 [52] httpuv_1.6.15          future.apply_1.11.3    goftest_1.2-3         
 [55] glue_1.8.0             nlme_3.1-166           promises_1.3.2        
 [58] grid_4.4.3             Rtsne_0.17             cluster_2.1.8.1       
 [61] reshape2_1.4.4         generics_0.1.3         gtable_0.3.6          
 [64] spatstat.data_3.1-4    tzdb_0.4.0             data.table_1.16.2     
 [67] hms_1.1.3              stringfish_0.18.0      utf8_1.2.4            
 [70] spatstat.geom_3.3-4    RcppAnnoy_0.0.22       ggrepel_0.9.6         
 [73] RANN_2.6.2             pillar_1.9.0           vroom_1.6.5           
 [76] spam_2.11-0            RcppHNSW_0.6.0         later_1.4.1           
 [79] splines_4.4.3          lattice_0.22-6         bit_4.5.0             
 [82] survival_3.8-3         deldir_2.0-4           tidyselect_1.2.1      
 [85] miniUI_0.1.1.1         pbapply_1.7-2          knitr_1.49            
 [88] gridExtra_2.3          scattermore_1.2        RhpcBLASctl_0.23-42   
 [91] xfun_0.49              matrixStats_1.4.1      leidenbase_0.1.31     
 [94] stringi_1.8.4          lazyeval_0.2.2         yaml_2.3.10           
 [97] pacman_0.5.1           evaluate_1.0.1         codetools_0.2-20      
[100] cli_3.6.3              uwot_0.2.2             RcppParallel_5.1.9    
[103] systemfonts_1.1.0      xtable_1.8-4           reticulate_1.42.0     
[106] munsell_0.5.1          globals_0.16.3         spatstat.random_3.3-2 
[109] png_0.1-8              spatstat.univar_3.1-1  parallel_4.4.3        
[112] dotCall64_1.2          listenv_0.9.1          viridisLite_0.4.2     
[115] scales_1.3.0           ggridges_0.5.6         crayon_1.5.3          
[118] rlang_1.1.4