5.4 标准化理论
5.4.1 简介
在上一章中，我们确定了重要的混杂因素和解释变量。scater允许在后续的统计模型中考虑这些变量，或者根据需要使用normaliseExprs()将它们标准化。
我们将探索如何通过简单的尺度因子（size-factor）标准化来校正文库大小，以消除一些混杂因素和解释变量的影响。
5.4.2 文库大小
ScRNA-seq数据的文库大小各异，来自每个细胞的总read量可能存在很大差异。一些定量方法（如Cufflinks、RSEM）在确定基因表达量时考虑了文库大小，因此不需要标准化。但是，如果使用其他定量方法，则必须通过将表达矩阵的每一列乘以或除以相对于其他细胞的文库大小估计出的“标准化因子”来校正。Bulk RNA测序已经开发了许多用于校正文库大小的方法，并且可以同样应用于scRNA测序（例如UQ、SF、CPM、RPKM、FPKM、TPM）。
5.4.3 标准化
5.4.3.1 CPM
标准化数据的最简单方法是将其转换为每百万计数（counts per million, CPM），即将每列中的值除以其该列总和，然后乘以1,000,000。请注意，为了校正细胞RNA总量，应将ERCC排除在总表达量的计算之外，因此我们将仅使用内源基因。R中CPM函数的示例：
> function (expr_mat, spikes = NULL)
{
norm_factor <- colSums(expr_mat[-spikes, ])
return(t(t(expr_mat)/norm_factor)) * 10^6
}
CPM的一个潜在缺点是如果样本包含高表达且在细胞间差异表达的基因，细胞中的分子总数可能取决于这些基因在细胞中是表达还是关闭，通过总分子数进行标准化可能会隐藏这些基因的差异表达和/或错误地为其他基因创建差异表达。
RPKM、FPKM和TPM是根据基因/转录本的长度进一步调整CPM的变型。
5.4.3.2 相对对数表达（Relative Log Expression, RLE）
尺度因子由DESeq（Anders and Huber 2010）提出并推广。首先计算每个基因在所有细胞中的几何平均值，每个细胞的尺度因子是基因表达与基因几何平均值之比的中位数。该方法的缺点是，由于它使用几何平均值，因此只能使用所有细胞中具有非零表达值的基因进行计算，因此不适用于大型低深度scRNA-seq。edgeR&scater将此方法称为“相对对数表达”（RLE）。R中的SF函数示例：
> function (expr_mat, spikes = NULL)
{
geomeans <- exp(rowMeans(log(expr_mat[-spikes, ])))
SF <- function(cnts) {
median((cnts/geomeans)[(is.finite(geomeans) & geomeans > 0)])
}
norm_factor <- apply(expr_mat[-spikes, ], 2, SF)
return(t(t(expr_mat)/norm_factor))
}
5.4.3.3 上四分位数（UQ）标准化
上四分位数（upper quartile, UQ）由（Bullard等，2010）提出，为每列除以每个文库count数的75%分位数。通常，计算出的分位数会根据细胞间的中位数进行缩放，以保持表达的绝对水平相对一致。该方法的一个缺点是，对于低深度scRNA-seq，大量未检测到的基因可能导致75%分位数为零（或接近于零）。可以通过扩展该思路并使用更高的分位数（99％分位数是scater中的默认值）或在计算75％分位数之前排除零来克服这一限制。R中的UQ函数示例：
> function (expr_mat, spikes = NULL)
{
UQ <- function(x) {
quantile(x[x > 0], 0.75)
}
uq <- unlist(apply(expr_mat[-spikes, ], 2, UQ))
norm_factor <- uq/median(uq)
return(t(t(expr_mat)/norm_factor))
}
5.4.3.4 M值的修剪平均值 (Trimmed Mean of M-values, TMM)
另一种方法称为TMM，是（Robinson and Oshlack 2010）提出的M值的加权修剪平均值。M值是基因在细胞间的log2倍数变化。以一个细胞作为参考，计算其他细胞与该细胞的M值，然后去掉最高和最低的30%来进行修剪，并对剩余值计算加权平均值，以解释对数变换对方差的影响。每个非参考细胞都乘以计算出的因子。此方法的两个潜在问题是修剪后剩余的非零基因不足，以及大多数基因没有差异表达的假设。
5.4.3.5 scran
scran包是CPM的变型，专门用于单细胞数据（L. Lun、Bach和Marioni，2016）。简言之，该方法通过将细胞合并（pooling）计算标准化因子来解决大量零值的问题。由于每个细胞存在于不同的集合中，因此可以使用线性代数方法从集合特异因子中反卷积出细胞特异因子。
5.4.3.6 下采样（Downsampling）
最后，校正文库大小的一个简单方法是对表达矩阵进行下采样，使每个细胞具有大致相同的分子总数。该方法的好处是，通过下采样将引入零值，从而消除由于检测到的基因数量不同而导致的偏差。但缺点是该过程是不确定的，因此每次运行下采样时，得到的表达矩阵都会略有不同。因此，通常必须对多次下采样进行分析以确保结果的鲁棒性。R中的下采样函数示例：
> function (expr_mat)
{
min_lib_size <- min(colSums(expr_mat))
down_sample <- function(x) {
prob <- min_lib_size/sum(x)
return(unlist(lapply(x, function(y) {rbinom(1, y, prob)})))
}
down_sampled_mat <- apply(expr_mat, 2, down_sample)
return(down_sampled_mat)
}
5.4.4 性能
为了比较不同标准化方法的效果，我们使用PCA图查看并通过scater的plotRLE()函数计算每个细胞的相对对数表达。Read较多（较少）的细胞大多数基因表达量高于（低于）中位数，导致整个细胞的RLE为正（负），而标准化后细胞的RLE接近于零。R中RLE函数的示例：
> function (expr_mat, spikes = NULL)
{
RLE_gene <- function(x) {
if (median(unlist(x)) > 0) {
log((x + 1)/(median(unlist(x)) + 1))/log(2)
}
else {
rep(NA, times = length(x))
}
}
if (!is.null(spikes)) {
RLE_matrix <- t(apply(expr_mat[-spikes, ], 1, RLE_gene))
}
else {
RLE_matrix <- t(apply(expr_mat, 1, RLE_gene))
}
cell_RLE <- apply(RLE_matrix, 2, median, na.rm = T)
return(cell_RLE)
}
注意：<
- RLE、TMM和UQ方法是为Bulk RNA测序数据开发的并依赖实验背景，可能不适用于scRNA测序数据，因为它们的基本假设可能被推翻。
-
scater充当edgeR中calcNormFactors函数的功能，它实现几种文库大小标准化方法，可以轻松地将这些方法中的任何一种应用于我们的数据。 -
edgeR对一些标准化方法进行了调整，这可能会导致与原始方法的结果略有不同，例如，edgeR和scater的“RLE”方法使用基于DESeq的“尺度因子”，可能会产生与DESeq/DESeq2包中的estimateSizeFactorsForMatrix方法不同的结果。此外，除非将所有细胞的lib.size设置为1，否则某些版本的edgeR将无法正确计算标准化因子。 - 对于CPM标准化，我们使用
scater的calculateCPM()函数。对于RLE、UQ和TMM，我们过去使用scater的normaliseExprs()函数（但该函数现已弃用）。对于scran，我们使用scran包来计算尺度因子（也可对SingleCellExperiment进行操作）并使用scater的normalize()来标准化数据。所有这些标准化函数都将结果保存到SCE对象的logcounts槽中。对于下采样，我们使用上方的函数。
往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（2）
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（3）
