27  单细胞聚类分析

单细胞聚类分析的目的在于从单细胞转录组数据中识别和分类细胞类型,以揭示细胞种类之间的差异以及它们在不同生理或病理状态下的变化。通过聚类算法将相似基因表达模式的细胞聚集在一起,可以识别出不同的细胞类型,从而有助于理解细胞组成和功能,并为研究细胞发展、疾病发生机制等提供重要线索。

27.1 加载R包

使用rm(list = ls())来清空环境中的所有变量。

library(tidyverse)
library(Seurat)
library(cowplot)

rm(list = ls())
options(stringsAsFactors = F)
options(future.globals.maxSize = 10000 * 1024^2)

grp_names <- c("Early Stage", "Late Stage")
grp_colors <- c("#8AC786", "#B897CA")
grp_shapes <- c(15, 16)

27.2 导入数据

  • seurat_obj数据来自于 章节 26,已经做过PCA降纬处理。

seurat_obj <- readRDS("./data/result/scRNA/Seurat_pca.RDS")

seurat_obj

27.3 确定PC纬度

在单细胞聚类分析过程中,对细胞数据进行PCA(主成分分析)以提取主成分(PC)是常见的预处理步骤。由于PC的选取直接关系到聚类结果的质量,因此选择恰当的PC数量至关重要。理想的PC数量应能充分保留原始数据中的关键信息,同时避免引入过多的冗余信息,以确保聚类结果的有效性和准确性。


# seurat_pca <- JackStraw(seurat_pca, num.replicate = 100)
# seurat_pca <- ScoreJackStraw(seurat_pca, dims = 1:30)
# JackStrawPlot(seurat_pca, dims = 1:20)

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

结果:可以观察到PC10 - 20周围是“肘部”,这表明大部分真实信号是在前20个PC捕获的,为了保留更多信息,我们选择40个PC。

27.4 寻找邻居

Seurat使用基于图的聚类方法,该方法采用了 K-最近邻方法,然后尝试将这个图划分为高度相互连接的“簇”。基于PCA空间中的欧氏距离来构建一个K-最近邻(K-nearest neighbor,KNN)图。


seurat_obj <- FindNeighbors(
  object = seurat_obj,
  dims = 1:40,
  verbose = FALSE)

seurat_obj

27.5 寻找聚类

Seurat将以优化标准模块化函数为目标,迭代地将细胞组合在一起。我们将使用FindClusters()函数来执行基于图的聚类。分辨率是一个重要的参数,它设定了下游聚类的“粒度”,并且需要针对每个独立实验进行优化。对于3000 - 5000个细胞的数据集,分辨率设置在0.4 - 1.4之间通常会产生良好的聚类效果。


seurat_cluster <- FindClusters(
  object = seurat_obj,
  resolution = seq(0.4, 1.4, 0.2),
  verbose = FALSE)

seurat_cluster

结果:不同分辨率的聚类结果可以通过seurat_obj@meta.data查看,会看到以RNA_snn_res开头(如果是整合多个来源数据的数据集可能是integrated_snn_res开头)。

  • 选择聚类分辨率,增加分辨率值会导致聚类数量增加,这通常对于更大的数据集是必要的。我们这里选择最大的1.4分辨率。

Idents(object = seurat_cluster) <- "RNA_snn_res.1.4"

seurat_cluster

27.6 可视化聚类

Seurat提供了几种非线性降维技术,如tSNE和UMAP,以可视化和探索这些数据集。这些算法的目标是学习数据的底层流形,以便在低维空间中将相似的细胞放置在一起。基于图的聚类中确定的细胞应该在这些降维图上共定位。作为UMAP和tSNE的输入,我们建议使用与聚类分析相同的主成分(PCs)作为输入。


seurat_cluster <- RunUMAP(
  object = seurat_cluster,
  reduction = "pca",
  dims = 1:40,
  verbose = FALSE)

seurat_cluster <- RunTSNE(
  object = seurat_cluster,
  reduction = "pca",
  dims = 1:40,
  check_duplicates = FALSE,
  verbose = FALSE)

DimPlot(seurat_cluster,
        reduction = "umap",
        label = TRUE,
        label.size = 6) + NoLegend()

DimPlot(seurat_cluster,
        reduction = "tsne",
        label = TRUE,
        label.size = 6) + NoLegend()

结果:总计识别出52个细胞类别。可以通过观察同一个细胞聚类簇它们是否聚集在一起,若它们是分散在不同区域,说明该分辨率下聚类效果不好请重新选择分辨率,但还是没有改善可能需要重新进行数据预处理等步骤了。比如,假定2640在某次聚类显示是一个细胞类群,但它们却是分散在很远的两端,这很大概率是聚类错误。

library(ggh4x)
library(ggunchull)

umap_pl <- ggplot(umap_data, aes(x = umap_1, y = umap_2, fill = ident, color = ident)) +
    geom_point(size = 0.5, show.legend = FALSE) +
    guides(color = "none", x = umap_axis, y = umap_axis) +
    scale_x_continuous(breaks = NULL) +
    scale_y_continuous(breaks = NULL) +  
    theme(aspect.ratio = 1,
         panel.background = element_blank(),
         panel.grid = element_blank(),
         axis.line = element_line(arrow = arrow(type = "closed")),
         axis.title = element_text(hjust = 0.05, face = "italic"))

umap_pl

27.7 输出结果

if (!dir.exists("./data/result/scRNA/")) {
  dir.create("./data/result/scRNA/", recursive = TRUE)
}

saveRDS(seurat_cluster, file = "./data/result/scRNA/Seurat_cluster.RDS", compress = TRUE)

27.8 总结

我们从降维后的数据中筛选出了能够代表足够多原始信息的主成分。这些被筛选出的主成分不仅保留了原始数据的关键特征,而且减少了数据的冗余和噪声,为后续的聚类分析提供了更加清晰和准确的数据基础。
系统信息
sessionInfo()
R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Asia/Shanghai
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] cowplot_1.1.3      Seurat_5.0.3       SeuratObject_5.0.2 sp_2.1-4          
 [5] lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
 [9] purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[13] ggplot2_3.5.1      tidyverse_2.0.0   

loaded via a namespace (and not attached):
  [1] deldir_2.0-4           pbapply_1.7-2          gridExtra_2.3         
  [4] rlang_1.1.3            magrittr_2.0.3         RcppAnnoy_0.0.22      
  [7] spatstat.geom_3.2-9    matrixStats_1.3.0      ggridges_0.5.6        
 [10] compiler_4.3.3         reshape2_1.4.4         png_0.1-8             
 [13] vctrs_0.6.5            pkgconfig_2.0.3        fastmap_1.1.1         
 [16] utf8_1.2.4             promises_1.3.0         rmarkdown_2.26        
 [19] tzdb_0.4.0             xfun_0.43              jsonlite_1.8.8        
 [22] goftest_1.2-3          later_1.3.2            spatstat.utils_3.0-4  
 [25] irlba_2.3.5.1          parallel_4.3.3         cluster_2.1.6         
 [28] R6_2.5.1               ica_1.0-3              spatstat.data_3.0-4   
 [31] stringi_1.8.4          RColorBrewer_1.1-3     reticulate_1.37.0     
 [34] parallelly_1.37.1      scattermore_1.2        lmtest_0.9-40         
 [37] Rcpp_1.0.12            knitr_1.46             tensor_1.5            
 [40] future.apply_1.11.2    zoo_1.8-12             sctransform_0.4.1     
 [43] httpuv_1.6.15          Matrix_1.6-5           splines_4.3.3         
 [46] igraph_2.0.3           timechange_0.3.0       tidyselect_1.2.1      
 [49] abind_1.4-5            rstudioapi_0.16.0      yaml_2.3.8            
 [52] spatstat.random_3.2-3  spatstat.explore_3.2-7 codetools_0.2-19      
 [55] miniUI_0.1.1.1         listenv_0.9.1          plyr_1.8.9            
 [58] lattice_0.22-6         shiny_1.8.1.1          withr_3.0.0           
 [61] ROCR_1.0-11            evaluate_0.23          Rtsne_0.17            
 [64] future_1.33.2          fastDummies_1.7.3      survival_3.7-0        
 [67] polyclip_1.10-6        fitdistrplus_1.1-11    pillar_1.9.0          
 [70] BiocManager_1.30.23    KernSmooth_2.23-22     renv_1.0.0            
 [73] plotly_4.10.4          generics_0.1.3         RcppHNSW_0.6.0        
 [76] hms_1.1.3              munsell_0.5.1          scales_1.3.0          
 [79] globals_0.16.3         xtable_1.8-4           glue_1.7.0            
 [82] lazyeval_0.2.2         tools_4.3.3            data.table_1.15.4     
 [85] RSpectra_0.16-1        RANN_2.6.1             leiden_0.4.3.1        
 [88] dotCall64_1.1-1        grid_4.3.3             colorspace_2.1-0      
 [91] nlme_3.1-164           patchwork_1.2.0        cli_3.6.2             
 [94] spatstat.sparse_3.0-3  spam_2.10-0            fansi_1.0.6           
 [97] viridisLite_0.4.2      uwot_0.2.2             gtable_0.3.5          
[100] digest_0.6.35          progressr_0.14.0       ggrepel_0.9.5         
[103] htmlwidgets_1.6.4      htmltools_0.5.8.1      lifecycle_1.0.4       
[106] httr_1.4.7             mime_0.12              MASS_7.3-60.0.1