一、ScaleData()简介
单细胞基因表达counts矩阵数据经过NormalizeData()归一化处理后，还需要进行scale标准化。ScaleData()函数将归一化的基因表达转换为Z分数（值以 0 为中心，方差为 1）。 它存储在 seurat_obj[['RNA']]@scale.data，用于下游的PCA降维。 默认是仅在高可变基因上运行标准化。
DefaultAssay(seurat_obj) <- 'RNA'
seurat_obj <- NormalizeData(seurat_obj) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()
从帮助文档中可以看出，ScaleData()实际上是对数据进行了scale和center两个步骤；
ScaleData(
object,
features = NULL,
assay = NULL,
vars.to.regress = NULL,
split.by = NULL,
model.use = "linear",
use.umi = FALSE,
do.scale = TRUE,
do.center = TRUE,
scale.max = 10,
block.size = 1000,
min.cells.to.block = 3000,
verbose = TRUE,
...
)
- PCA聚类前对所有细胞的每个基因（或高可变基因）进行标准化基因表达值：
- 这在下游分析中给予同等的权重，因此高表达的基因不会占主导地位；标准化的目的是为了避免那些表达很高的基因的变化，掩盖了表达丰富较低但也有意义的基因，所以把这些基因缩放在合适的范围，以便能看到所有有意义基因的差异。
-
ScaleData函数对数据进行Z-score标准化（Z-score normalization ）处理：
- 使每个基因在所有细胞的表达量均值为 0
- 使每个基因在所有细胞中的表达量方差为 1
-
ScaleData可以选择回归掉不需要的变异来源（regress out unwanted variation）
- 细胞的聚类结果可能受细胞周期状态的影响，可以回归掉细胞周期的影响；
- 常见的几种不感兴趣的变异源：技术噪音（每个细胞检测到的转录本数量、线粒体转录本百分比），批次效应，细胞周期状态等；
- 剔除这种变异可以改善下游分析；
- 在vars.to.regress参数中添加要回归的变量，它们会针对每个特征单独回归，然后对结果残差进行标准化和居中；
主要参数有：
| 参数 | 说明 |
|---|---|
| features | 要标准化/居中的特征基因，默认是高可变基因 |
| vars.to.regress | 要回归的变量， 例如，nUMI 或percent.mito |
| verbose | 显示ScaleData过程的进度条 |
scale.max是指scale标准化后的数据最大值，默认值为 10。如出现极端值，如scale.data数值大于等于10，一律用10表示；减少仅在极少数细胞中基因表达值的影响。
示例：
seurat_obj <- ScaleData(seurat_obj, features = VariableFeatures(seurat_obj))
seurat_obj <- ScaleData(seurat_obj, features = all.genes)
seurat_obj <- ScaleData(seurat_obj, vars.to.regress = c("nFeature_RNA", "percent.mito","S.Score", "G2M.Score"))
ScaleData()现在合并了以前称为RegressOut()函数的功能（根据所提供的变量进行回归，然后标准化残差）。 要使用回归功能，只需将要剔除的变量传递给vars.to.regress 参数。将center设置为TRUE ，将通过减去该基因的平均表达来使每个细胞的基因表达值居中。 如果center 为TRUE，scale设置为TRUE，将居中的基因表达值除以它们的标准差来标准化每个基因的表达水平；如scale设置为FALSE，则除以它们的均方根。
在解析ScaleData()源码之前先熟悉几个基本概念。
二、 统计学中的Z Score
2.1 什么是Z Score？
基础统计学中有一个基本的问题，即已知一组数据符合正态分布，同时给定算术平均值和标准差的情况下，如何计算该组数据的Z score。
Z score的概念是指原始数据距离均值有多少个标准差。当以标准差为单位进行测量时，Z score衡量的是一个数值偏离总体均值以上或以下多少个标准差。如果原始数值高于均值，则 Z score得分为正，如果低于均值，则Z score为负。
Z score其实是标准正态分布（Standard Normal Distribution），即平均值μ=0，标准差 σ=1 的正态分布。SND标准正态分布的直方图如下所示：
这里的横轴就是 z 分数（Z Score），也叫做标准分数（Standard Score）：z score = (x – μ) / σ
事实上，任何一个正态分布都可以通过标准化（Standardization）变成标准正态分布。z分数主要分布范围从-3个标准差(落在正态分布曲线的最左边)到+3个标准差(落在正态分布曲线的最右边)。而标准化的过程，就是按照上面这个公式把随机变量 x 变为 z 分数。 z 分数的值代表 x 的不同取值偏离平均值 μ 多少个标准差 σ。比如，当 z 分数等于 1 时，说明该值偏离平均值 1 个标准差σ。
2.2 Z Score计算公式
z score计算公式：z score = (x – μ) / σ，其中 x 是原始数值，μ 是总体平均值，σ 是总体标准偏差。
当总体均值和总体标准差未知时，可以使用样本均值 (x̄) 和样本标准差 (s) 作为总体值的估计值来计算标准分数。
2.3 如何解释 Z Score？
Z Score的值告诉你与平均值相差多少标准差。如果Z Score等于 0，则它在平均值上。
Z Score为正表示原始数值高于平均值。例如，如果Z Score等于 +1，则它比平均值高 1 个标准差。
负Z Score表明原始数值低于平均水平。例如，如果Z Score等于 -2，则它比平均值低 2 个标准差。
解释Z Score的另一种方法是创建标准正态分布（也称为 z 分数分布或概率分布）。
图2说明了任何标准正态分布 (SND) 的重要特征：
- SND（即 z 分布）始终与原始数值分布的形状相同。例如，如果原始数值的分布是正态分布的，那么 z 分数的分布也是如此。
- 任何 SND 的平均值总是=0。
-
任何 SND 的标准偏差总是=1。因此，原始数值的一个标准偏差（无论这是什么原始值）转换为1 个Z Score单位。
图2.标准正态分布 (SND) 和正态分布
2.4 为什么Z score很重要？
通过将正态分布的值（原始数值）转换为Z score来标准化它们，是非常有用的，因为：
-
研究人员计算Z score出现在标准正态分布内的概率；
SND允许研究人员计算从分布（即样本）中随机获得Z score的概率。例如，有 68% 的概率从平均值中随机选择一个介于 -1 到 +1 标准差之间的Z score（见图 3）。从平均值中随机选择一个介于 -1.96 和 +1.96 标准差之间的分数的概率为 95%（见图 3）。如果随机选择原始数值的可能性低于 5%，那么这是一个具有统计显著性的结果。
图 3.以百分比表示的标准正态分布 (SND) 的比例 - 并使我们能够比较来自不同样本的两个Z score（可能具有不同的均值和标准差）。
这正是我们对转录组表达矩阵进行scale的原因。z-score是常用的标准正态分布化的方法。
基因表达矩阵中的一个基因在不同细胞中的表达量值，通过求z-score值，转换为标准正态分布，经过处理的数据的均值为0，标准差为1，因此z-score也称为零-均值规范化。而z-score考虑到了不同细胞对表达量的影响，计算z-score时，消除到了表达值的平均水平和偏离度的影响。
每个基因都是标准正态分布，在下游分析中给予同等的权重，高表达的基因不会占主导地位。
2.5 R语言的Z score计算
R语言的Z score计算是通过scale()函数求得，Seurat包中ScaleData()函数也基本参照了scale()函数的功能。
scale方法中的两个参数：center和scale
- center和scale默认为真，即T或者TRUE
- center为真表示数据中心化：数据集中的各项数据减去数据集的均值
如有数据集1, 2, 3, 6, 3，其均值为3
那么中心化之后的数据集为1-3,2-3,3-3,6-3,3-3，即：-2,-1,0,3,0 - scale为真表示数据标准化：指中心化之后的数据在除以数据集的标准差，即数据集中的各项数据减去数据集的均值再除以数据集的标准差。
如有数据集1, 2, 3, 6, 3，其均值为3，其标准差为1.87
那么标准化之后的数据集为(1-3)/1.87,(2-3)/1.87,(3-3)/1.87,(6-3)/1.87,(3-3)/1.87，即：-1.069,-0.535,0,1.604,0
示例：
# 限定输出小数点后数字的位数为3位
options(digits=3)
data <- c(1, 2, 3, 6, 3)
data
# [1] 1 2 3 6 3
# 数据中心化
scale(data, center=T, scale=F)
# [,1]
# [1,] -2
# [2,] -1
# [3,] 0
# [4,] 3
# [5,] 0
# attr(,"scaled:center")
# [1] 3
# 数据标准化
scale(data, center=T, scale=T)
[,1]
# [1,] -1.069
# [2,] -0.535
# [3,] 0.000
# [4,] 1.604
# [5,] 0.000
# attr(,"scaled:center")
# [1] 3
# attr(,"scaled:scale")
# [1] 1.87
三、 ScaleData()源码解析
3.1 查看ScaleData()源码
我们查看下ScaleData()的源码，调用的代码比较多，主代码在preprocessing.R文件中，进行scale标准化有两种处理方式：
1）如果细胞表达矩阵的类型为dgCMatrix和dgTMatrix稀疏矩阵，则调用FastSparseRowScale方法；FastSparseRowScale方法是C++函数，是主程序通过Rcpp方式调用。单细胞10X数据调用此方法。
2）若为其他数据格式 ，则调用FastRowScale方法， 是R代码函数。 另外如果vars.to.regress参数不为空，还会调用RegressOutMatrix相关函数。
FastSparseRowScale方法：
Eigen::MatrixXd FastSparseRowScale(Eigen::SparseMatrix<double> mat, bool scale = true, bool center = true,
double scale_max = 10, bool display_progress = true){
mat = mat.transpose();
Progress p(mat.outerSize(), display_progress);
Eigen::MatrixXd scaled_mat(mat.rows(), mat.cols());
for (int k=0; k<mat.outerSize(); ++k){
p.increment();
double colMean = 0;
double colSdev = 0;
for (Eigen::SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
{
colMean += it.value();
}
colMean = colMean / mat.rows();
if (scale == true){
int nnZero = 0;
if(center == true){
for (Eigen::SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
{
nnZero += 1;
colSdev += pow((it.value() - colMean), 2);
}
colSdev += pow(colMean, 2) * (mat.rows() - nnZero);
}
else{
for (Eigen::SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
{
colSdev += pow(it.value(), 2);
}
}
colSdev = sqrt(colSdev / (mat.rows() - 1));
}
else{
colSdev = 1;
}
if(center == false){
colMean = 0;
}
Eigen::VectorXd col = Eigen::VectorXd(mat.col(k));
scaled_mat.col(k) = (col.array() - colMean) / colSdev;
for(int s=0; s<scaled_mat.col(k).size(); ++s){
if(scaled_mat(s,k) > scale_max){
scaled_mat(s,k) = scale_max;
}
}
}
return scaled_mat.transpose();
}
FastRowScale方法：
FastRowScale <- function(
mat,
center = TRUE,
scale = TRUE,
scale_max = 10
) {
# inspired by https://www.r-bloggers.com/a-faster-scale-function/
if (center) {
rm <- rowMeans2(x = mat, na.rm = TRUE)
}
if (scale) {
if (center) {
rsd <- rowSds(mat, center = rm)
} else {
rsd <- sqrt(x = rowSums2(x = mat^2)/(ncol(x = mat) - 1))
}
}
if (center) {
mat <- mat - rm
}
if (scale) {
mat <- mat / rsd
}
if (scale_max != Inf) {
mat[mat > scale_max] <- scale_max
}
return(mat)
}
3.2 示例演示
通过一个小案例，我们看下FastRowScale()和FastSparseRowScale()具体的运算。
示例1：测试FastRowScale
step1：构造一个细胞基因表达矩阵mat(15*10），列为细胞，行为基因；
step2：通过下面的测试，FastRowScale(mat)完全等价于t(scale(t(mat)))的结果；
问题1：为什么不直接调用scale()函数，而需要重新构造FastRowScale()函数？
在https://www.r-bloggers.com/2016/02/a-faster-scale-function/中详细说明了此问题，因为FastRowScale()调用了matrixStats包的colSds函数，整个运算速度要比scale快很多，运算性能得到优化。江湖武功唯快不破。
# Tests for scaling data
# --------------------------------------------------------------------------------
context("Fast Scale Data Functions")
mat <- matrix(rnorm(n = 10*15), nrow = 10, ncol = 15)
# should be the equivalent of t(scale(t(mat)))
test_that("Fast implementation of row scaling returns expected values", {
expect_equal(t(scale(t(mat)))[1:10, 1:15], FastRowScale(mat))
expect_equal(t(scale(t(mat), center = FALSE))[1:10, 1:15],
FastRowScale(mat, center = FALSE))
expect_equal(t(scale(t(mat), scale = FALSE))[1:10, 1:15],
FastRowScale(mat, scale = FALSE))
expect_equal(t(scale(t(mat), scale = FALSE, center = F))[1:10, 1:15],
FastRowScale(mat, scale = FALSE, center = F))
mat.clipped <- FastRowScale(mat, scale_max = 0.2)
expect_true(max(mat.clipped, na.rm = T) >= 0.2)
})
示例2：测试FastSparseRowScale
FastSparseRowScale(mat)同样完全等价于t(scale(t(mat)))的结果；也是为了追求运算速度，改写成C++代码，解决庞大的稀疏矩阵scale运算速度。
# should be the equivalent of t(scale(t(mat)))
mat <- rsparsematrix(10, 15, 0.1)
test_that("Fast implementation of row scaling returns expected values", {
expect_equal(t(scale(t(as.matrix(mat))))[1:10, 1:15], FastSparseRowScale(mat, display_progress = FALSE),
check.attributes = FALSE)
expect_equal(t(scale(t(as.matrix(mat)), center = FALSE))[1:10, 1:15],
FastSparseRowScale(mat, center = FALSE, display_progress = FALSE),
check.attributes = FALSE)
expect_equal(t(scale(t(as.matrix(mat)), scale = FALSE))[1:10, 1:15],
FastSparseRowScale(mat, scale = FALSE, display_progress = FALSE),
check.attributes = FALSE)
expect_equal(t(scale(t(as.matrix(mat)), scale = FALSE, center = F))[1:10, 1:15],
FastSparseRowScale(mat, scale = FALSE, center = F, display_progress = FALSE),
check.attributes = FALSE)
mat.clipped <- FastSparseRowScale(mat, scale_max = 0.2, display_progress = F)
expect_true(max(mat.clipped, na.rm = T) >= 0.2)
})
四、ScaleData()注意事项
4.1 在seurat下游分析中，什么时候用slot ="data"，什么时候用slot ="scale.data"?
- 对于单个样本，通过NormalizeData()归一化处理，使基因表达矩阵符合正态分布，归一化的数据储存在"RNA" assay的 seurat_obj[['RNA']]@data中;
- 归一化的基因表达矩阵，进一步进行scale标准化，使其符合标准正态分布，scale标准化的数据储存在"RNA" assay的 seurat_obj[['RNA']]@scale.data中。
我们也注意到seurat_obj[['RNA']]@data全是非负数，而且是针对基因矩阵的所有基因；而seurat_obj[['RNA']]@scale.data则有正负数，默认情况，只针对高可变基因进行scale标准化；
那么，我们在seurat下游分析中，什么情况使用data，什么时候使用scale.data：
1）下游分析中的PCA线性降维聚类，umap、tsne聚类均是应用高可变基因的scale.data进行后续分析的；
2）在基因可视化分析中，FeaturePlot、FeatureScatter、VlnPlot、DotPlot等函数默认slot ="data"，只有DoHeatmap()默认使用slot = "scale.data"，多个基因跨细胞比较；
我的理解是：
1）只是比较单个基因在不同样本/处理、细胞类型，使用slot ="data"就足够，而且避免使用"scale.data"，导致关注的基因在scale.data矩阵中未找到，关注的基因不是高可变基因；
另外，data的值域范围比scale.data也大很多，便于通过颜色梯度直观看出基因表达强弱；
2）比较多个基因的表达水平则考虑scale标准化的矩阵，如，DoHeatmap，但是为什么DotPlot函数默认slot ="data"?
DotPlot函数内置有scale参数，且默认scale为True，也就是你输入data矩阵，DotPlot函数还是会进行scale处理；
3）FindAllMarkers()找差异基因是默认slot ="data"，它是针对所有基因找差异基因，而不是高可变基因集；
4.2 多样本整合分析中，进行scale标准化要注意DefaultAssay的类型
单个样本只有RNA模式，即1)原始表达矩阵：seurat_obj@assays$RNA@counts; 2）NormalizeData归一化矩阵：seurat_obj@assays$RNA@data；3）scale标准化矩阵：seurat_obj@assays$RNA@scale.data；
但是在多样本，有RNA模式和integrated模式；下面的实例对比选用DefaultAssay为"RNA"和"integrated"，进行scale标准化，下游分析会有很大的差异。
DefaultAssay(pbmc_seurat) <- "RNA"
pbmc_seurat <- NormalizeData(pbmc_seurat, verbose = F)
pbmc_seurat <- FindVariableFeatures(pbmc_seurat, selection.method = "vst", nfeatures = 2000, verbose = F)
pbmc_seurat <- ScaleData(pbmc_seurat, verbose = F)
pbmc_seurat <- RunPCA(pbmc_seurat, npcs = 30, verbose = F)
pbmc_seurat <- RunUMAP(pbmc_seurat, reduction = "pca", dims = 1:30, verbose = F)
在整合前，数据集的UMAP图显示清晰的分离。
DimPlot(pbmc_seurat,reduction = "umap") + plot_annotation(title = "10k 3' PBMC and 10k 5' PBMC cells, before integration")
现在让我们将assay更改为'integrated'，并在整合分析中做同样的事情（它已经标准化并选择了 HVG）：
DefaultAssay(pbmc_seurat) <- "integrated"
pbmc_seurat <- ScaleData(pbmc_seurat, verbose = F)
pbmc_seurat <- RunPCA(pbmc_seurat, npcs = 30, verbose = F)
pbmc_seurat <- RunUMAP(pbmc_seurat, reduction = "pca", dims = 1:30, verbose = F)
最后，绘制integrated模式下的UMAP图，数据非常好地整合在一起。
DefaultAssay设置为"RNA"和"integrated"，它们的data和scale.data不相同；DefaultAssay为"RNA"，只是两个样本简单的merge，没有去除批次。
# before: 对未整合 (RNA) 检测进行归一化、HVG 发现、缩放、PCA 和 UMAP:
DefaultAssay(pbmc_seurat) <- "RNA"
gene_exp_matrix_data1 <- pbmc_seurat@assays$RNA@data
dim(gene_exp_matrix_data1)
# [1] 36601 19190
gene_exp_matrix_data2 <- pbmc_seurat@assays$RNA@scale.data
dim(gene_exp_matrix_data2)
# [1] 2000 19190
# after: 将分析更改为集成，并在集成分析中做同样的事情（它已经标准化并选择了HVG）：
DefaultAssay(pbmc_seurat) <- "integrated"
gene_exp_matrix_data3 <- pbmc_seurat@assays$integrated@data
dim(gene_exp_matrix_data3)
# [1] 2000 19190
gene_exp_matrix_data4 <- pbmc_seurat@assays$integrated@scale.data
dim(gene_exp_matrix_data4)
# [1] 2000 19190
testthat::expect_equal(gene_exp_matrix_data1[row.names(gene_exp_matrix_data3),],gene_exp_matrix_data3)
# Error: gene_exp_matrix_data1[row.names(gene_exp_matrix_data3), ] not equal to `gene_exp_matrix_data3`.
4.3 问题：ScaleData需要按样本分别进行标准化处理吗
比如，有两个样本，merge在一起，用harmony进行批次处理进行整合分析，但是需要
1）按样本split分别进行scale处理，还是2）不区分样本信息，两个样本的细胞混合在一起，一起进行scale处理？
seurat_harmony <- merge(srat_3p,srat_5p)
seurat_harmony <- NormalizeData(seurat_harmony)
seurat_harmony <- FindVariableFeatures(seurat_harmony)
# 按样本分别进行scale
seurat_harmony1 <- ScaleData(seurat_harmony, split.by = "orig.ident")
# 作为整体一起scale
seurat_harmony2 <- ScaleData(seurat_harmony)
testthat::expect_equal(seurat_harmony1@assays$RNA@data,seurat_harmony2@assays$RNA@data)
# seurat_harmony1和seurat_harmony2的scale.data信息不一致
testthat::expect_equal(seurat_harmony1@assays$RNA@scale.data,seurat_harmony2@assays$RNA@scale.data)
# Error: seurat_harmony1@assays$RNA@scale.data not equal to seurat_harmony2@assays$RNA@scale.data.
另外在https://www.singlecellcourse.org/scrna-seq-dataset-integration.html#liger-3-vs-5-10k-pbmc教程中，liger进行数据整合时，scale步骤是区分样本的（split.by = "orig.ident"）
我的理解是如果样本间存在批次效应，最好是按样本分别进行scale处理，但是很多教程在scale步骤，并没有加split.by = "orig.ident"这个参数，有点不理解，希望能解释这块分析的大佬帮吗解答，万分感激。
后记：在github查看了ScaleData相关代码，有些写有"split.by = batch"；这块应该是要求不严谨，可以写split.by，也可以不加，如果有批次效应，最好split.by = batch处理一下。
# 代码1
seurat_object <- NormalizeData(seurat_object)
seurat_object <- FindVariableFeatures(seurat_object)
seurat_object <- ScaleData(seurat_object,split.by = batch)
seurat_object <- RunPCA(seurat_object, verbose = FALSE)
# merge
seurat_object.integrated <- RunHarmony(seurat_object,batch, assay.use = au)
# 代码2
seurat_object.integrated <- SCTransform(seurat_object, vars.to.regress = batch, verbose = FALSE)
seurat_object <- FindVariableFeatures(seurat_object)
seurat_object <- ScaleData(seurat_object,split.by = batch)
