一、ScaleData()簡介
單細胞基因表達counts矩陣數(shù)據(jù)經(jīng)過NormalizeData()歸一化處理后,還需要進行scale標準化。ScaleData()函數(shù)將歸一化的基因表達轉(zhuǎn)換為Z分數(shù)(值以 0 為中心,方差為 1)。 它存儲在 seurat_obj[['RNA']]@scale.data,用于下游的PCA降維。 默認是僅在高可變基因上運行標準化。
DefaultAssay(seurat_obj) <- 'RNA'
seurat_obj <- NormalizeData(seurat_obj) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()
從幫助文檔中可以看出,ScaleData()實際上是對數(shù)據(jù)進行了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聚類前對所有細胞的每個基因(或高可變基因)進行標準化基因表達值:
- 這在下游分析中給予同等的權(quán)重,因此高表達的基因不會占主導(dǎo)地位;標準化的目的是為了避免那些表達很高的基因的變化,掩蓋了表達豐富較低但也有意義的基因,所以把這些基因縮放在合適的范圍,以便能看到所有有意義基因的差異。
-
ScaleData函數(shù)對數(shù)據(jù)進行Z-score標準化(Z-score normalization )處理:
- 使每個基因在所有細胞的表達量均值為 0
- 使每個基因在所有細胞中的表達量方差為 1
-
ScaleData可以選擇回歸掉不需要的變異來源(regress out unwanted variation)
- 細胞的聚類結(jié)果可能受細胞周期狀態(tài)的影響,可以回歸掉細胞周期的影響;
- 常見的幾種不感興趣的變異源:技術(shù)噪音(每個細胞檢測到的轉(zhuǎn)錄本數(shù)量、線粒體轉(zhuǎn)錄本百分比),批次效應(yīng),細胞周期狀態(tài)等;
- 剔除這種變異可以改善下游分析;
- 在vars.to.regress參數(shù)中添加要回歸的變量,它們會針對每個特征單獨回歸,然后對結(jié)果殘差進行標準化和居中;
主要參數(shù)有:
| 參數(shù) | 說明 |
|---|---|
| features | 要標準化/居中的特征基因,默認是高可變基因 |
| vars.to.regress | 要回歸的變量, 例如,nUMI 或percent.mito |
| verbose | 顯示ScaleData過程的進度條 |
scale.max是指scale標準化后的數(shù)據(jù)最大值,默認值為 10。如出現(xiàn)極端值,如scale.data數(shù)值大于等于10,一律用10表示;減少僅在極少數(shù)細胞中基因表達值的影響。
示例:
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()現(xiàn)在合并了以前稱為RegressOut()函數(shù)的功能(根據(jù)所提供的變量進行回歸,然后標準化殘差)。 要使用回歸功能,只需將要剔除的變量傳遞給vars.to.regress 參數(shù)。將center設(shè)置為TRUE ,將通過減去該基因的平均表達來使每個細胞的基因表達值居中。 如果center 為TRUE,scale設(shè)置為TRUE,將居中的基因表達值除以它們的標準差來標準化每個基因的表達水平;如scale設(shè)置為FALSE,則除以它們的均方根。
在解析ScaleData()源碼之前先熟悉幾個基本概念。
二、 統(tǒng)計學(xué)中的Z Score
2.1 什么是Z Score?
基礎(chǔ)統(tǒng)計學(xué)中有一個基本的問題,即已知一組數(shù)據(jù)符合正態(tài)分布,同時給定算術(shù)平均值和標準差的情況下,如何計算該組數(shù)據(jù)的Z score。
Z score的概念是指原始數(shù)據(jù)距離均值有多少個標準差。當以標準差為單位進行測量時,Z score衡量的是一個數(shù)值偏離總體均值以上或以下多少個標準差。如果原始數(shù)值高于均值,則 Z score得分為正,如果低于均值,則Z score為負。
Z score其實是標準正態(tài)分布(Standard Normal Distribution),即平均值μ=0,標準差 σ=1 的正態(tài)分布。SND標準正態(tài)分布的直方圖如下所示:

這里的橫軸就是 z 分數(shù)(Z Score),也叫做標準分數(shù)(Standard Score):z score = (x – μ) / σ
事實上,任何一個正態(tài)分布都可以通過標準化(Standardization)變成標準正態(tài)分布。z分數(shù)主要分布范圍從-3個標準差(落在正態(tài)分布曲線的最左邊)到+3個標準差(落在正態(tài)分布曲線的最右邊)。而標準化的過程,就是按照上面這個公式把隨機變量 x 變?yōu)?z 分數(shù)。 z 分數(shù)的值代表 x 的不同取值偏離平均值 μ 多少個標準差 σ。比如,當 z 分數(shù)等于 1 時,說明該值偏離平均值 1 個標準差σ。
2.2 Z Score計算公式
z score計算公式:z score = (x – μ) / σ,其中 x 是原始數(shù)值,μ 是總體平均值,σ 是總體標準偏差。
當總體均值和總體標準差未知時,可以使用樣本均值 (x?) 和樣本標準差 (s) 作為總體值的估計值來計算標準分數(shù)。
2.3 如何解釋 Z Score?
Z Score的值告訴你與平均值相差多少標準差。如果Z Score等于 0,則它在平均值上。
Z Score為正表示原始數(shù)值高于平均值。例如,如果Z Score等于 +1,則它比平均值高 1 個標準差。
負Z Score表明原始數(shù)值低于平均水平。例如,如果Z Score等于 -2,則它比平均值低 2 個標準差。
解釋Z Score的另一種方法是創(chuàng)建標準正態(tài)分布(也稱為 z 分數(shù)分布或概率分布)。
圖2說明了任何標準正態(tài)分布 (SND) 的重要特征:
- SND(即 z 分布)始終與原始數(shù)值分布的形狀相同。例如,如果原始數(shù)值的分布是正態(tài)分布的,那么 z 分數(shù)的分布也是如此。
- 任何 SND 的平均值總是=0。
-
任何 SND 的標準偏差總是=1。因此,原始數(shù)值的一個標準偏差(無論這是什么原始值)轉(zhuǎn)換為1 個Z Score單位。
圖2.標準正態(tài)分布 (SND) 和正態(tài)分布
2.4 為什么Z score很重要?
通過將正態(tài)分布的值(原始數(shù)值)轉(zhuǎn)換為Z score來標準化它們,是非常有用的,因為:
-
研究人員計算Z score出現(xiàn)在標準正態(tài)分布內(nèi)的概率;
SND允許研究人員計算從分布(即樣本)中隨機獲得Z score的概率。例如,有 68% 的概率從平均值中隨機選擇一個介于 -1 到 +1 標準差之間的Z score(見圖 3)。從平均值中隨機選擇一個介于 -1.96 和 +1.96 標準差之間的分數(shù)的概率為 95%(見圖 3)。如果隨機選擇原始數(shù)值的可能性低于 5%,那么這是一個具有統(tǒng)計顯著性的結(jié)果。
圖 3.以百分比表示的標準正態(tài)分布 (SND) 的比例 - 并使我們能夠比較來自不同樣本的兩個Z score(可能具有不同的均值和標準差)。
這正是我們對轉(zhuǎn)錄組表達矩陣進行scale的原因。z-score是常用的標準正態(tài)分布化的方法。
基因表達矩陣中的一個基因在不同細胞中的表達量值,通過求z-score值,轉(zhuǎn)換為標準正態(tài)分布,經(jīng)過處理的數(shù)據(jù)的均值為0,標準差為1,因此z-score也稱為零-均值規(guī)范化。而z-score考慮到了不同細胞對表達量的影響,計算z-score時,消除到了表達值的平均水平和偏離度的影響。
每個基因都是標準正態(tài)分布,在下游分析中給予同等的權(quán)重,高表達的基因不會占主導(dǎo)地位。
2.5 R語言的Z score計算
R語言的Z score計算是通過scale()函數(shù)求得,Seurat包中ScaleData()函數(shù)也基本參照了scale()函數(shù)的功能。
scale方法中的兩個參數(shù):center和scale
- center和scale默認為真,即T或者TRUE
- center為真表示數(shù)據(jù)中心化:數(shù)據(jù)集中的各項數(shù)據(jù)減去數(shù)據(jù)集的均值
如有數(shù)據(jù)集1, 2, 3, 6, 3,其均值為3
那么中心化之后的數(shù)據(jù)集為1-3,2-3,3-3,6-3,3-3,即:-2,-1,0,3,0 - scale為真表示數(shù)據(jù)標準化:指中心化之后的數(shù)據(jù)在除以數(shù)據(jù)集的標準差,即數(shù)據(jù)集中的各項數(shù)據(jù)減去數(shù)據(jù)集的均值再除以數(shù)據(jù)集的標準差。
如有數(shù)據(jù)集1, 2, 3, 6, 3,其均值為3,其標準差為1.87
那么標準化之后的數(shù)據(jù)集為(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
示例:
# 限定輸出小數(shù)點后數(shù)字的位數(shù)為3位
options(digits=3)
data <- c(1, 2, 3, 6, 3)
data
# [1] 1 2 3 6 3
# 數(shù)據(jù)中心化
scale(data, center=T, scale=F)
# [,1]
# [1,] -2
# [2,] -1
# [3,] 0
# [4,] 3
# [5,] 0
# attr(,"scaled:center")
# [1] 3
# 數(shù)據(jù)標準化
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()的源碼,調(diào)用的代碼比較多,主代碼在preprocessing.R文件中,進行scale標準化有兩種處理方式:
1)如果細胞表達矩陣的類型為dgCMatrix和dgTMatrix稀疏矩陣,則調(diào)用FastSparseRowScale方法;FastSparseRowScale方法是C++函數(shù),是主程序通過Rcpp方式調(diào)用。單細胞10X數(shù)據(jù)調(diào)用此方法。
2)若為其他數(shù)據(jù)格式 ,則調(diào)用FastRowScale方法, 是R代碼函數(shù)。 另外如果vars.to.regress參數(shù)不為空,還會調(diào)用RegressOutMatrix相關(guān)函數(shù)。
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:構(gòu)造一個細胞基因表達矩陣mat(15*10),列為細胞,行為基因;
step2:通過下面的測試,F(xiàn)astRowScale(mat)完全等價于t(scale(t(mat)))的結(jié)果;
問題1:為什么不直接調(diào)用scale()函數(shù),而需要重新構(gòu)造FastRowScale()函數(shù)?
在https://www.r-bloggers.com/2016/02/a-faster-scale-function/中詳細說明了此問題,因為FastRowScale()調(diào)用了matrixStats包的colSds函數(shù),整個運算速度要比scale快很多,運算性能得到優(yōu)化。江湖武功唯快不破。
# 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)))的結(jié)果;也是為了追求運算速度,改寫成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()歸一化處理,使基因表達矩陣符合正態(tài)分布,歸一化的數(shù)據(jù)儲存在"RNA" assay的 seurat_obj[['RNA']]@data中;
- 歸一化的基因表達矩陣,進一步進行scale標準化,使其符合標準正態(tài)分布,scale標準化的數(shù)據(jù)儲存在"RNA" assay的 seurat_obj[['RNA']]@scale.data中。
我們也注意到seurat_obj[['RNA']]@data全是非負數(shù),而且是針對基因矩陣的所有基因;而seurat_obj[['RNA']]@scale.data則有正負數(shù),默認情況,只針對高可變基因進行scale標準化;
那么,我們在seurat下游分析中,什么情況使用data,什么時候使用scale.data:
1)下游分析中的PCA線性降維聚類,umap、tsne聚類均是應(yīng)用高可變基因的scale.data進行后續(xù)分析的;
2)在基因可視化分析中,F(xiàn)eaturePlot、FeatureScatter、VlnPlot、DotPlot等函數(shù)默認slot ="data",只有DoHeatmap()默認使用slot = "scale.data",多個基因跨細胞比較;
我的理解是:
1)只是比較單個基因在不同樣本/處理、細胞類型,使用slot ="data"就足夠,而且避免使用"scale.data",導(dǎo)致關(guān)注的基因在scale.data矩陣中未找到,關(guān)注的基因不是高可變基因;
另外,data的值域范圍比scale.data也大很多,便于通過顏色梯度直觀看出基因表達強弱;
2)比較多個基因的表達水平則考慮scale標準化的矩陣,如,DoHeatmap,但是為什么DotPlot函數(shù)默認slot ="data"?
DotPlot函數(shù)內(nèi)置有scale參數(shù),且默認scale為True,也就是你輸入data矩陣,DotPlot函數(shù)還是會進行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)
在整合前,數(shù)據(jù)集的UMAP圖顯示清晰的分離。
DimPlot(pbmc_seurat,reduction = "umap") + plot_annotation(title = "10k 3' PBMC and 10k 5' PBMC cells, before integration")

現(xiàn)在讓我們將assay更改為'integrated',并在整合分析中做同樣的事情(它已經(jīng)標準化并選擇了 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圖,數(shù)據(jù)非常好地整合在一起。

DefaultAssay設(shè)置為"RNA"和"integrated",它們的data和scale.data不相同;DefaultAssay為"RNA",只是兩個樣本簡單的merge,沒有去除批次。
# before: 對未整合 (RNA) 檢測進行歸一化、HVG 發(fā)現(xiàn)、縮放、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: 將分析更改為集成,并在集成分析中做同樣的事情(它已經(jīng)標準化并選擇了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)不區(qū)分樣本信息,兩個樣本的細胞混合在一起,一起進行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進行數(shù)據(jù)整合時,scale步驟是區(qū)分樣本的(split.by = "orig.ident")

我的理解是如果樣本間存在批次效應(yīng),最好是按樣本分別進行scale處理,但是很多教程在scale步驟,并沒有加split.by = "orig.ident"這個參數(shù),有點不理解,希望能解釋這塊分析的大佬幫嗎解答,萬分感激。
后記:在github查看了ScaleData相關(guān)代碼,有些寫有"split.by = batch";這塊應(yīng)該是要求不嚴謹,可以寫split.by,也可以不加,如果有批次效應(yīng),最好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)

