ixxmu
commented
1 year ago Seurat 深度总结1---常规流程+默认参数的探讨 by TOP生物信息
1. 为什么写这篇帖子
1.1 Seurat 的重要性
1.2 网络上教程的局限性
2. Setup the Seurat Object
3. Standard pre-processing workflow
3.1 QC and selecting cells for further analysis
4. Normalizing the data
4.1 标准化的方法
4.2 scale.factor 的选择
5. Identification of highly variable features (feature selection)
5.1 找高变基因的方法
5.2 vst 的参数 loess span 选择
5.3 调整高变基因数据
6. Scaling the data
小结一下
往期回顾
1. 为什么写这篇帖子
1.1 Seurat 的重要性
单细胞分析最流行的 R 包 全网教程最多的单细胞分析的流程
1.2 网络上教程的局限性
Seurat 官网流程的搬运工（即 NormalizeData ---> FindVariableFeatures ---> ScaleData ---> RunPCA ---> FindNeighbors ---> FindClusters ---> RunUMAP/RunTSNE），或者分享 Seurat 的进阶流程 SCTransform 没有真正探讨各个默认参数对于分析的影响，例如： FindVariableFeatures 除了 2000，能否选择 1000?3000?5000？ ScaleData 是否应该做回归？如果回归，要回归什么参数？ PC 数目应该如何选择？ 分辨率选择多少？
本篇推文，笔者将基于
Seurat 的官网流程(https://satijalab.org/seurat/articles/pbmc3k_tutorial.html) 自己对于单细胞分析的经验 数篇参考文献抛砖引玉，讨论一下单细胞 Seuat 的流程中，NormalizeData--->FindVariableFeatures ---> ScaleData 的参数选择情况
Schneider et al. J Transl Genet Genom. 2021;5:37-49
2. Setup the Seurat Object
官网中的第一步，使用的是PBMC 的数据(https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz)，读取标准的 Cellranger 下机的文件，然后建立 Seurat 对象。
library(dplyr)
library(Seurat)
library(patchwork)
# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/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)
pbmc
3. Standard pre-processing workflow
3.1 QC and selecting cells for further analysis
3.2.1 细胞层面过滤
在单细胞分析中，最常见的质控指标包括：
UMI counts/cell: 即 nCount_RNA，每个细胞的 UMI 数目 太大：可能是双细胞或者多细胞的集合 genes/cell: 即 nFeature_RNA，每个细胞所检测到的基因数目 太小：可能是无效细胞 太大：可能是双细胞或者多细胞的集合 percent.mt: 映射到线粒体基因组的读数百分比 低质量/垂死的细胞往往表现出广泛的线粒体污染
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
在这一步中，最为关键的就是如何选择质控的标准。质控的标准一般是结合小提琴图 [图 3] + 散点图 [图 4]+既往文献确定
3.2.2 基因层面过滤 (optional)
在质控这一步，也可以进行基因方便的质控，例如删除0表达值的基因和在少于10个细胞中表达的基因
# Output a logical vector for every gene on whether the more than zero counts per cell
# Extract counts
counts <- GetAssayData(object = pbmc, slot = "counts")
# Output a logical vector for every gene on whether the more than zero counts per cell
nonzero <- counts > 0
# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
keep_genes <- Matrix::rowSums(nonzero) >= 10
# Only keeping those genes expressed in more than 10 cells
pbmc_filtered_counts <- counts[keep_genes, ]
# Reassign to filtered Seurat object
pbmc_filtered_seurat <- CreateSeuratObject(pbmc_filtered_counts, meta.data = pbmc@meta.data)
参考：Single-cell RNA-seq analysis workshop ( https://github.com/hbctraining/scRNA-seq )
4. Normalizing the data
标准化的代码很容易:
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
但是其实NormalizeData这个函数中，normalization.method 和 scale.factor 都是可以设置的参数。
4.1 标准化的方法
normalization.method 有三种方法：
LogNormalize：Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale.factor. This is then natural-log transformed using log1p CLR：Applies a centered log ratio transformation RC: Relative counts. Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale.factor. No log-transformation is applied. For counts per million (CPM) set scale.factor = 1e6
从 [图 5] 可以看出，三种方法中，LogN 相对于 CLR 和 RC 两种方法而言，能够有更多的 cluster [图 5A]。但是三种方法注释出来的 cluster 是相似的 [图 5B]。但是，三种方法标准化后，cluster 会略有不同 [图 5C]。
ARI: ARI 值是衡量两个聚类方案之间的相似性/不相似性的指标，值为 1 代表完全相似，值为 0 代表完全不相似。
备注：BCL_bc: B cells; END_en: endothelial cells; EPI_ep: epithelial cancer cells; FIB_fi: fibroblasts; MAC_ma: macrophages
在降维图上，三种不同方法产生的 cluster 方案之间出现一些差异：
cluster 的识别上：有些 cluster 被一种方法识别出来，而其他方法则没有（[图 6A] 中的 cluster 3 和 cluster 12 以及[图 6C] 中的 cluster 10） 相同 cluster 的细分上：某些方法会将一组细胞定义为一个 cluster，但是另外的方法可能会将这群 cluster 分开，例如，成纤维细胞在 LogNormalize 方法中被分成三组，但在 RC 和 CLR 方法中只有两组，而且 RC 和 CLR 方法产生的两组也不相同 [图 6 A-F]。
4.2 scale.factor 的选择
相较于 10000 的默认值来说，增加 10 倍或者减少 10 倍产生的结果变化不大，主要变化的地方在于细胞亚群的数量上。
在降维图上，不同比例因子产生的聚类方案之间出现一些差异，例如：
当比例因子降低时，产生的额外的 cluster（cluster 14 [图 8B]）是一小群最初在基线分析中被确定为巨噬细胞的细胞（cluster 2 [图 8 A]）。在新的方案下，这些细胞转而注释为 B 细胞 [图 8A-B]； 比例因子增加时丢失的 cluster（cluster 12 [图 8A, D]）被重新归类为内皮细胞（cluster 10 [图 8C, F]）。 这种变化细胞类型名字（如：原来的注释信息为巨噬，后来为 B 细胞）的 cluster 比例很小 [图 7B] 。新的 cluser 方案与基线相比，ARI 降低 [图 7C]。
备注：
[图 8C] 中的红线作者标偏了，应该指向 cluster 10 作者使用的注释方法为，使用了预先设定的 gene list, 并把这些 gene list 存储为 txt 格式。再根据现有 cluster 的基因与设定的基因集的相似性来进行注释。
5. Identification of highly variable features (feature selection)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
5.1 找高变基因的方法
在 Seurat 中，有三种找高变基因的方法，分别是 vst、mean.var.plot 和 dispersion。
三种方法能够找到的高变基因各不相同（vst: 1125、 mean.var.plot: 818、dispersion: 1125），其中三种方法共同找到的基因数为 439 (25%) ，有一半被确定为可变的基因只能被一种方法确定（871/1750=50%） [图 9A]。 即使有这些不同的可变基因集，检测到的 cluster 数目是相似的[图 9B]，只有 2%的细胞被不一致地注释 [图 9C]。 [图 9D]显示，mean.var.plot 和 dispersion 相较于 vst 的 ARI 较低，根据[图 10] 的 TSNE 降维图可以看出，细胞分配到 cluster 的情况确实不同。基于 TSNE 的可视化，我们根据它们之间的视觉接近程度，手动将 cluster 分配到大的分群中[图 10A-C]。超过 96%的细胞被置于一致的大的 cluster 中。
5.2 vst 的参数 loess span 选择
因为 vst 方法使用局部多项式回归（loess），所以需要一个输入参数是 loess span。在这里对比了 loess span 的默认值 0.3，以及 0.1 和 0.5 的区别。同时也比较了不同的分箱策略的影响：
[图 11]中 A，D，G 为 cluster 数；B，E，H 为细胞注释相似度；C，F，I 为 ARI 值
A, B, C 比较 loess span 为 0.3（基线）到 0.1 和 0.5 的差异 D, E, F 比较了 bin 数为 20（基线）到 1，10 和 100 的差异 G, H, I 比较了使用等宽（基线）到等频的分箱方法
结果显示，loess span、bin 数以及分箱方法对于 cluster 数目、细胞注释和 ARI 的值没有明显的影响。
5.3 调整高变基因数据
作者在这里，使用了三种不同的高变基因数目，分别为所有检测到的基因的 6.7%、1%和 20%。原因为，如果检测到15,000个基因，6.7% 就会产生一个约1000个可变基因的列表，但是这种选择没有生物学依据，只是为了对比不同数目的可变基因对于聚类的影响。
[图 12] 条形图显示了使用 1,125 （6.7%）的可变基因列表与 168 个基因（1%）和 3,358 个基因（20%）相比在 cluster 数目、相似度和 ARI 的差异。结果提示：
可变基因列表减少到总基因的1%对聚类有很大影响 [图 12B-C] 将可变基因列表增加到所有基因的20%，产生了与使用6.7%基因的基线分析类似的 cluster 方案 [图 12A-C]
6. Scaling the data
ScaleData 也是一句代码：
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
因为使用 all.genes 来做ScaleData的话，占用内存非常多。ScaleData的主要是为了后续的 PCA 降维和使用DoHeatmap画图，但是画热图一般后续使用 ComplexHeatmap 或者 Pheatmap 会比较多，因此只需要对高变基因进行 Scale 即可：
## 默认就是使用高变基因
pbmc <- ScaleData(pbmc)
在ScaleData这一步，比较困扰的就是 vars.to.regress 这个参数，到底什么样的参数需要回归是困扰新手的一个重要问题，在不同文章中也存在不同的设置，例如：
Jin SZ et al. Cell Res. 2020
Chen YP et al. Cell Res. 2020
Bassez A et al. Nat Med. 2021
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
在文章中，作者比较了不回归参数、回归 total UMI count、线粒体基因以及回归 total UMI count 和线粒体基因的区别，结果显示得到的 cluster 数目非常类似，且97%的细胞都可以被指定都同一个细胞类型[图 10A 和图 10B]。但是，[图 10C]显示，相比于回归 total UMI count、线粒体基因两项来说，回归线粒体和不回归参数的 ARI 低于 0.8，且回归线粒体和不回归的聚类情况、分群数目都比较相似，因此这提示了增加回归线粒体基因并不会对结果有太大的影响。
ScaleData中需要一个 maximum scale value，默认为 10 或 50。作者比较发现，10 和 100 对于聚类没有太大的影响 [图 11]。
小结一下
标准化的方法会显著影响聚类结果，且显著高于改变 scale factor 的结果； 增加高变基因数目并不会显著影响聚类效果（~1000 到~3000）,但是降低到~1000 到~125 会影响聚类效果； 回归 UMI count 对于结果有比较大的影响，但是回归线粒体数的影响较小； 改变 maximum scale value 不会显著影响聚类效果
后续 RunPCA--->FindNeighbors--->FindClusters--->RunUMAP 中，还有几个非常重要的参数会影响实际的聚类效果，包括 PC 数、k-parameter、prune parameter 和 resolution parameter，考虑到这些参数在后续的降维聚类中息息相关，因此，我将在下一篇帖子中统一进行讨论。
https://mp.weixin.qq.com/s/VJ4lHzjXcJwefmVlJAi0ig