初探单细胞分析 — 标准化与降维聚类分群的理解

生信技能树jimmy

发布于 2024-03-04 11:34:24

2190

发布于 2024-03-04 11:34:24

文章被收录于专栏：单细胞天地单细胞天地

最近曾老师讲授的单细胞课程四天课程开课啦，通过该课程，我们学习了很多单细胞相关的知识和实战经验。下面的内容简单通俗的讲解一下关于第一天课程中对于标准化数据和降维聚类分群的理解，特别是对基因进行降维的思路使我们能更方便快速的处理单细胞数据。

标准化数据

类似于bulk RNA-seq，single-cell RNA-seq 的原始count数据也是需要进行标准化的。

首先是细胞测序深度进行校正，通过NormalizeData()函数实现，大致步骤是先进行同一文库大小的校正，然后进行对数化操作。

接着是进行z-score的转换，使其成为正态分布，通过ScaleData()函数进行转化。但是由于大部分基因都是不表达的，那些在大部分细胞中都不表达的基因对于后续分析没有帮助，反而增加噪音，以及增加分析时间。所以，在进行ScaleData之前，进行了高变基因的查找（通过FindVariableFeatures()函数实现，默认是2000个高变基因）。ScaleData默认是对高变基因进行转换。

这里梳理一下当进行两个标准化步骤之后，数据结构的变化（基于seurat V5）：

原始count数据：pbmc@assaysRNA@layerscounts
NormalizeData之后的数据：pbmc@assaysRNA@layersdata
ScaleData之后的数据：pbmc@assaysRNA@layersscale.data

下面是实操代码：

# 读入矩阵
pbmc.data <- Read10X(data.dir = "./filtered_gene_bc_matrices/hg19/")
# 创建 seurat 对象，并进行简单的质控
pbmc <- CreateSeuratObject(counts = pbmc.data, 
                           project = "pbmc3k", 
                           min.cells = 3, min.features = 200)
# NormalizeData目的为消除不同细胞测序深度的影响。
# 首先对基因的reads数进行了同一文库大小的校正，然后进行对数化操作
pbmc <- NormalizeData(pbmc)
# 找高变基因，默认的高变基因个数（nfeatures）是2000
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst",
                                nfeatures = 2000)
# ScaleData默认对高变基因的表达量数值进行z-score的转换
pbmc <- ScaleData(pbmc)

降维聚类分群

假设每个基因都是数据的一个维度，由于基因一般有2万多个，所以就存在2万多个维度，这对于分析是很困难的。所以需要保留对整体数据有影响的基因，去除基因变化小的基因来简化分析。上面提到单细胞数据中大部分基因都是不表达的，在课堂上我们查看了一下所分析的示例只有很少的基因的count数不为0（大约5%），也就是95%的基因的count数都为0。

因此，在标准流程中是对基因进行了三次降维处理，再进行聚类分群后，得到我们人眼能够分析的二维可视化降维聚类图。

下面来看看细节。在上述标准化的分析中，已经进行了一次降维（最初的2万多个基因 -> 2000个高变基因）。第二次降维是主成分分析（PCA），通过RunPCA()函数实现，默认的PCs数目是50个。如果有需要，也可以利用elbow point去选择合适的PCs数目，但一般用默认参数即可。第二次降维的变化为：2000个高变基因 -> 50个PCs特征值。然后用PCA的结果进行聚类分群（简单理解为对每个细胞标上一个分群标签），以及用PCA的结果进行UMAP或者t-SNE方法的二维可视化。可视化是最后一次降维，使数据能以一个二维的方式进行分析展示（50个PCs特征值 -> 二维可视化）。

下面是降维聚类分群的实操代码：

# RunPCA的输入是ScaleData()后的数据，默认的维度（npcs）是50
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), 
                  verbose = FALSE)
# 聚类分群：FindNeighbors 是找到基因之间距离
# 默认取PCA的前10个维度进行分析(dims = 1:10)
pbmc <- FindNeighbors(pbmc, dims = 1:10, verbose = FALSE)
# FindClusters 是聚类。越高的分辨率（resolution）得到的聚类数目会更多
pbmc <- FindClusters(pbmc, resolution = 0.5, verbose = FALSE)
table(pbmc$seurat_clusters)
# 可视化降维
pbmc <- RunUMAP(pbmc, dims = 1:10, umap.method = "uwot", metric = "cosine")

理解这个过程后，相信大家应该都不会忘记seurat单细胞的标准分析步骤了吧，每一步都是有其意义滴。以及也解答一些我刚开始疑惑的问题。比如先进行聚类分群（FindNeighbors；FindClusters），还是先进行RunUMAP。答案是都可以，因为大家是独立的，并且都是用PCA的结果，聚类分群相当于对每个细胞给一个分群标签，RunUMAP相当于得到了可视化的二维坐标。

参数探索

在降维时候，涉及到一些参数，比如FindVariableFeatures的高变基因个数（默认为2000个）；RunPCA的PCs数目（默认是50个）。建议大家用默认的就好，或者参考一些文献。其实没有什么对错，是一个权衡吧：高变基因个数越多，PCs的数目越多，保留的数据信息也就越多，更有可能引入噪音，运行速度也越慢。但是也不能太少，否则会丢失很多数据信息。

聚类分群的时候也有一些参数，比如FindNeighbors的dim参数和FindClusters的resolution参数，是与最后的分群数目有关的。dim决定了取多少个PCA的维度进行分析，resolution是分辨率，都是越高，得到的分群数目会越多。该参数的设置与你的课题需要相关，所以曾老师的代码也是对不同的resolution都进行了计算，我们只需要结合结果图去选择合适的参数。

在这里我对降维时候的高变基因个数和PCs数，这两个参数进行了一个简单的探索，结果从UMAP图来看是差不多的。不过具体还是看大家的分析需要啦，如果分析的结果有问题，再想想是否是参数或者其他的问题了。

左图是高变基因1000；PCs数目为10的UMAP图。(pca_pbmc_v1000.pca_10.pdf)
右图是高变基因1000；PCs数目为30的UMAP图。(pca_pbmc_v1000.pca_30.pdf)

左图是高变基因3000；PCs数目为10的UMAP图。(pca_pbmc_v3000.pca_10.pdf)
右图是高变基因3000；PCs数目为50的UMAP图。(pca_pbmc_v3000.pca_50.pdf)

参数探索的代码如下：

library(Seurat)
library(ggplot2)
library(tidyverse)
pbmc.data <- Read10X(data.dir = "./filtered_gene_bc_matrices/hg19/")
## Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, 
                           project = "pbmc3k", 
                           min.cells = 3, min.features = 200)
## Normalizing the data
pbmc <- NormalizeData(pbmc)

variable_nums <- c(1000,2000,3000)
npcs_nums <- c(10,20,30,40,50)
for (i in variable_nums){
  for (j in npcs_nums){
  tmp_sce <- FindVariableFeatures(pbmc, selection.method = "vst",
                                  nfeatures = i)
  tmp_sce <- ScaleData(tmp_sce)
  tmp_sce <- RunPCA(tmp_sce, features = VariableFeatures(object = tmp_sce), 
                 verbose = FALSE, npcs= j)
  assign(paste("pbmc_v", i,".pca_",j, sep = ""), tmp_sce)
  }
}

for (i in variable_nums){
  for (j in npcs_nums){
    tmp_sce <- get(paste("pbmc_v", i,".pca_",j, sep = ""))
    tmp_sce <- FindNeighbors(tmp_sce, dims = 1:10, verbose = FALSE)
    tmp_sce <- FindClusters(tmp_sce, resolution = 0.5, verbose = FALSE)
    table(tmp_sce$seurat_clusters)
    pdf(paste0("pca_pbmc_v",i,".pca_",j,".pdf"))
    p1 <- DimPlot(tmp_sce, reduction = "pca")
    tmp_sce <- RunUMAP(tmp_sce, dims = 1:10, umap.method = "uwot", metric = "cosine")
    p2 <- DimPlot(tmp_sce,label = T,reduction = "umap")
    print(p1)
    print(p2)
    dev.off()
    assign(paste("pbmc_v", i,".pca_dim",j, sep = ""), tmp_sce)
  }
}