完整的单细胞分析流程——数据标化（normalization）
动机
通常在单细胞RNA测序数据中观察到文库之间测序覆盖率的系统差异。它们通常是由细胞间的cDNA捕获或PCR扩增效率方面的技术差异引起的，这归因于用最少的起始材料难以实现一致的文库制备。标准化旨在消除这些差异，以使它们不干扰细胞之间表达谱的比较。这样可以确保在细胞群体中观察到的任何异质性或差异表达都是由生物学而不是技术偏倚引起的。
在这一点上，规范化和批次校正之间的区别需要注意。归一化的发生与批次结构无关，并且仅考虑技术偏差，而批次矫正仅在批次之间发生，并且必须同时考虑技术偏差和生物学差异。技术偏倚倾向于以相似的方式或至少以与它们的生物物理特性（例如长度，GC含量）有关的方式影响基因，而批次之间的生物学差异可能是高度不可预测的。这样，这两个任务涉及不同的假设，并且通常涉及不同的计算方法（尽管某些软件包旨在一次执行两个步骤，例如zinbwave）。因此，避免混淆“标准化”和“批次校正”的数据非常重要，因为这些数据通常表示不同的事物。
我们将主要关注缩放标准化，这是最简单和最常用的标准化策略。这涉及将每个细胞的所有计数除以特定于细胞的比例因子，通常称为“大小因子”。这里的假设是，任何细胞特异性偏倚（例如，捕获或扩增效率）均会通过缩放该细胞的预期平均数来同等地影响所有基因。每个细胞的大小因子表示该细胞中相对偏差的估计，因此，将其计数除以其大小因子应消除该偏差。然后可以将所得的“归一化数据”用于下游分析，例如聚类和降维。为了演示，我们将使用来自scRNAseq软件包的数据集。
#--- loading ---#
library(scRNAseq)
sce.zeisel <- ZeiselBrainData()
library(scater)
sce.zeisel <- aggregateAcrossFeatures(sce.zeisel,
id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))
#--- gene-annotation ---#
library(org.Mm.eg.db)
rowData(sce.zeisel)$Ensembl <- mapIds(org.Mm.eg.db,
keys=rownames(sce.zeisel), keytype="SYMBOL", column="ENSEMBL")
#--- quality-control ---#
stats <- perCellQCMetrics(sce.zeisel, subsets=list(
Mt=rowData(sce.zeisel)$featureType=="mito"))
qc <- quickPerCellQC(stats, percent_subsets=c("altexps_ERCC_percent",
"subsets_Mt_percent"))
sce.zeisel <- sce.zeisel[,!qc$discard]
sce.zeisel
## class: SingleCellExperiment
## dim: 19839 2816
## metadata(0):
## assays(1): counts
## rownames(19839): 0610005C13Rik 0610007N19Rik ... mt-Tw mt-Ty
## rowData names(2): featureType Ensembl
## colnames(2816): 1772071015_C02 1772071017_G12 ... 1772063068_D01
## 1772066098_A12
## colData names(10): tissue group # ... level1class level2class
## reducedDimNames(0):
## altExpNames(2): ERCC repeat
文库归一化
文库大小归一化是执行缩放归一化的最简单策略。 我们将文库的大小定义为每个细胞中所有基因的计数总和，假定其预期值随任何细胞特异性偏倚而缩放。 然后，在定义比例常数的情况下，每个细胞的“库大小因子”直接与其库大小成正比，从而使所有细胞的平均大小因子等于1。此定义可确保归一化的表达值与原始计数处于相同规模 ——这对解释很有用——尤其是在处理转换后的数据时。
library(scater)
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
summary(lib.sf.zeisel)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.176 0.568 0.868 1.000 1.278 4.084
在Zeisel脑数据中，文库大小因子在细胞之间的差异最大10倍。 这是scRNA-seq数据覆盖范围变异的典型表现。
hist(log10(lib.sf.zeisel), xlab="Log10[Size factor]", col='grey80')
从Zeisel脑数据集中的库大小得出的大小因子分布。
严格来说，文库大小因子的使用是假设任何一对细胞之间的差异表达（DE）基因中都没有“不平衡”。也就是说，基因的一个子集的任何上调都可以通过不同基因子集中的相同下调幅度来抵消。这样可以通过避免合成效应来确保文库大小是相对于细胞特异性相对偏倚的无偏估计。但是，平衡的DE通常在scRNA-seq应用中不存在，这意味着文库大小归一化可能无法为下游分析产生准确的归一化表达值。
在实践中，标准化的准确性不是探索性scRNA-seq数据分析的主要考虑因素。成分偏差通常不会影响细胞群的分离，而只会影响细胞群或细胞类型之间的对数倍数变化的幅度——向着程度较小的方向。因此，库大小归一化通常在许多应用中都是足够的，这些应用的目的是识别细胞群和定义每个细胞群的top标记。
通过反卷积均一化
如前所述，当样本之间存在任何不平衡的差异表达时，就会出现成分偏差。以两个细胞举例，其中单个基因X与细胞B相比在细胞A中被上调。这种上调意味着（i）更多的测序资源用于A中的X，从而当每个细胞的总文库大小通过实验确定时（例如，由于文库量化）；其他的非差异基因的覆盖率降低，或（ii）当为X分配更多的读数或UMI时，A的文库大小增加，从而增加了文库大小因子，并为所有非DE基因产生了较小的归一化表达值。在这两种情况下，最终结果是，与B相比，A中的非DE基因将被错误地下调。
对于大量RNA测序数据分析，消除成分偏差是一个经过充分研究的问题。可以使用DESeq2包中的estimateSizeFactorsFromMatrix（）函数或edgeR包中的calcNormFactors（）函数来执行规范化。这些假设大多数基因不是细胞之间的DE。假设两个细胞之间多数非DE基因之间的计数大小的任何系统性差异都代表了偏差，该偏差用于计算适当的大小因子以将其去除。
然而，由于存在大量的低计数和零计数，单细胞数据应用这些bulk归一化方法可能会有问题。为了克服这个问题，我们汇总了许多细胞的计数以进行准确的大小因子估算。然后，将基于库的大小因子“分解”为基于细胞的大小因子，以标准化每个细胞的表达谱。如下所示，这是使用来自scran的computeSumFactors（）函数执行的。
library(scran)
set.seed(100)
clust.zeisel <- quickCluster(sce.zeisel)
table(clust.zeisel)
## clust.zeisel
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 170 254 441 178 393 148 219 240 189 123 112 103 135 111
deconv.sf.zeisel <- calculateSumFactors(sce.zeisel, cluster=clust.zeisel)
summary(deconv.sf.zeisel)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.119 0.486 0.831 1.000 1.321 4.509
我们使用带有quickCluster（）的预聚类步骤，其中每个聚类中的细胞分别进行归一化，并且将大小因子重新缩放以在各个聚类中具有可比性。这避免了在整个种群中大多数基因都是非DE的假设-在成对的簇之间仅需要非DE多数，这对于高度异质的种群来说是一个较弱的假设。默认情况下，quickCluster（）将基于irlba软件包中的方法对PCA使用近似算法。近似值依赖于随机初始化，因此我们需要设置随机种子（通过set.seed（））以实现可重现性。
我们看到，解卷积大小因子与图7.2中的库大小因子表现出特定于细胞类型的偏差。这与由细胞类型之间强烈的差异表达引入的成分偏倚的存在是一致的。去卷积大小因子的使用针对这些偏差进行调整，以提高下游应用程序的归一化精度。
plot(lib.sf.zeisel, deconv.sf.zeisel, xlab="Library size factor",
ylab="Deconvolution size factor", log='xy', pch=16,
col=as.integer(factor(sce.zeisel$level1class)))
abline(a=0, b=1, col="red")
Zeisel脑数据集中每个细胞的反卷积大小因子，与从库大小得出的等效大小因子相比。 红线对应于两个大小因子之间的同一性
准确的归一化对于涉及对每个基因统计信息的估计和解释的过程而言最重要。 例如，成分偏倚会通过系统性地将对数倍数变化沿一个方向或另一个方向转移来破坏DE分析。 但是，对于基于细胞的分析（如聚类分析），与简单的库大小归一化相比，它往往提供的好处较少。 成分偏差的存在已经暗示了表达谱的巨大差异，因此更改标准化策略不太可能影响聚类过程的结果。
通过spike-in进行归一化
spike-in归一化基于以下假设：向每个细胞中添加了相同量的spike-in RNA。spike-in转录本覆盖范围的系统差异仅归因于细胞特异性偏差，例如捕获效率或测序深度。为了消除这些偏差，我们通过缩放“ spike-in size factor”来均衡细胞间的spike-in覆盖范围。与以前的方法相比，spike-in归一化不需要系统的生物学假设（即，没有许多DE基因）。取而代之的是，它假定将掺入的spike-in转录本（i）以恒定的水平添加到每个细胞中，并且（ii）以与内源基因相同的相对方式响应偏倚。
实际上，如果需要关注单个细胞的总RNA含量差异，并且必须保留在下游分析中，则应使用加标归一化。对于给定的细胞，内源RNA总量的增加不会增加其spike-in大小因子。这确保了总RNA含量在群体间的表达差异不会在缩放时消除。相比之下，上述其他标准化方法将仅将总RNA含量的任何变化解释为偏差的一部分，并将其消除。
举个例子，在不同亲和力的T细胞受体配体刺激后，在涉及T细胞活化的不同数据集上使用spike-in归一化
library(scRNAseq)
sce.richard <- RichardTCellData()
sce.richard <- sce.richard[,sce.richard$`single cell quality`=="OK"]
sce.richard
## class: SingleCellExperiment
## dim: 46603 528
## metadata(0):
## assays(1): counts
## rownames(46603): ENSMUSG00000102693 ENSMUSG00000064842 ...
## ENSMUSG00000096730 ENSMUSG00000095742
## rowData names(0):
## colnames(528): SLX-12611.N701_S502. SLX-12611.N702_S502. ...
## SLX-12612.i712_i522. SLX-12612.i714_i522.
## colData names(13): age individual ... stimulus time
## reducedDimNames(0):
## altExpNames(1): ERCC
我们应用computeSpikeFactors（）方法来估计所有细胞的spike-in大小因子。 通过使用与librarySizeFactors（）中相同的推理，将每个细胞的总spike-in计数转换为大小因子来定义。 scaling将随后消除细胞间spike-in覆盖率的任何差异。
sce.richard <- computeSpikeFactors(sce.richard, "ERCC")
summary(sizeFactors(sce.richard))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.125 0.428 0.627 1.000 1.070 23.316
我们观察到每种处理条件下spike-in大小因子和解卷积大小因子之间存在正相关关系（图7.3），表明它们在测序深度和捕获效率上捕获了相似的技术偏倚。 但是，我们还观察到，就亲和力或时间的增加而言，对T细胞受体的刺激不断增加，导致spike-in因子相对于文库大小因子而言有所降低。 这与刺激过程中生物合成活性和总RNA含量的增加一致，这减少了每个文库中的相对spike-in覆盖率（从而减少了spike-in大小因子），但增加了内源基因的覆盖率（因此增加了文库大小因子）。
to.plot <- data.frame(
DeconvFactor=calculateSumFactors(sce.richard),
SpikeFactor=sizeFactors(sce.richard),
Stimulus=sce.richard$stimulus,
Time=sce.richard$time
)
ggplot(to.plot, aes(x=DeconvFactor, y=SpikeFactor, color=Time)) +
geom_point() + facet_wrap(~Stimulus) + scale_x_log10() +
scale_y_log10() + geom_abline(intercept=0, slope=1, color="red")
将spike-in归一化的大小因子与T细胞数据集中所有细胞的库大小因子作图。 每个图代表不同的配体处理，每个点是根据刺激时间按颜色着色的细胞。
两组尺寸因子之间的差异对下游解释产生了实际影响。 如果将spike-in 大小因子应用于计数矩阵，则未刺激细胞中的表达值将按比例放大，而受刺激细胞中的表达将按比例缩小。 但是，如果使用反卷积大小因子，则会发生相反的情况。 当我们在标准化策略之间切换时，这可以表现为条件之间DE的大小和方向的变化，如下Malat1所示（图7.4）。
# See below for explanation of logNormCounts().
sce.richard.deconv <- logNormCounts(sce.richard, size_factors=to.plot$DeconvFactor)
sce.richard.spike <- logNormCounts(sce.richard, size_factors=to.plot$SpikeFactor)
gridExtra::grid.arrange(
plotExpression(sce.richard.deconv, x="stimulus",
colour_by="time", features="ENSMUSG00000092341") +
theme(axis.text.x = element_text(angle = 90)) +
ggtitle("After deconvolution"),
plotExpression(sce.richard.spike, x="stimulus",
colour_by="time", features="ENSMUSG00000092341") +
theme(axis.text.x = element_text(angle = 90)) +
ggtitle("After spike-in normalization"),
ncol=2
)
使用反卷积大小因子（左）或spike-in大小因子（右）进行归一化后，Malat1的对数归一化表达值的分布。 细胞通过配体亲和力分层，并在刺激后的时间着色。
总RNA含量是否相关——因此，归一化策略的选择——取决于生物学假设。 在大多数情况下，总RNA含量的变化并不令人感兴趣，可以通过应用文库大小或反卷积因子将其标准化。 但是，如果总RNA的差异与感兴趣的生物学过程（例如细胞周期活性或T细胞活化）有关，这可能并不总是合适的。 Spike-in归一化将保留这些差异，以使生物学分组之间的表达变化具有正确的符号。
然而！ 无论我们是否关心总RNA含量，使用spike-in大小因子对spike-in转录本进行标准化都是至关重要的。 从内源基因计数中计算得出的大小因子不应该应用于spike-in转录本，正是因为前者捕获了后者未能实现的总RNA含量差异。 尝试使用基于基因的大小因子对spike-in计数进行归一化将导致过度归一化和错误的定量。 因此，如果需要归一化spike-in数据，我们必须为spike-in转录本计算一组单独的大小因子。 这是由诸如
modelGeneVarWithSpikes（）之类的函数自动执行的。
Applying the size factors
缩放和对数转化
一旦计算出大小因子，就可以使用scater中的logNormCounts（）函数为每个细胞计算归一化的表达值。 这是通过将每个基因/spike-in转录本的计数除以该细胞的合适大小因子来完成的。 该函数还对归一化后的值进行对数转换，从而创建了一个称为“ logcounts”的新assay。 这些对数值将在以下各章中作为我们下游分析的基础。
set.seed(100)
clust.zeisel <- quickCluster(sce.zeisel)
sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)
sce.zeisel <- logNormCounts(sce.zeisel)
assayNames(sce.zeisel)
## [1] "counts" "logcounts"
对数转换很有用，因为对数值的差异表示基因表达的对数倍变化。这在基于欧几里得距离的下游过程中很重要，下游过程包括许多形式的聚类和降维。通过对对数转换后的数据进行操作，我们确保这些过程基于基因表达的对数倍变化来测量细胞之间的距离。比如，一个在细胞类型A中平均表达量为50，在细胞类型B中表达量为10的基因，或在A中为1100，B中为1000的基因，对数转化可以展现出具有强烈相对差异，因此会关注前者。
在进行对数转换时，我们通常会添加一个伪计数以避免值为零。对于低丰度基因，较大的伪计数将有效地将细胞之间的对数倍变化缩小至零，这意味着下游的高维分析将更多地由高丰度基因的表达差异来驱动。相反，较小的伪计数将增加低丰度基因的相对贡献。常见的做法是使用1的伪计数，原因很简单，即实用的原因是它保留原始矩阵中的稀疏性（即原矩阵中的零在变换后仍为零）。除大多数病理情况外，此方法在所有情况下均有效。
顺便说一句，伪计数的增加是出于将尺寸因子居中统一的动机。这确保了伪计数和规范化的表达式值都在同一范围内。伪计数为1可以解释为每个基因的额外reads或UMI。实际上，居中意味着随着计数深度的提高，伪计数的收缩效果减小。这正确地确保了表达的对数倍变化的估计（例如，根据细胞组之间对数值的差异）随着覆盖范围的扩大而变得越来越准确。相反，如果将恒定的伪计数应用于类似百万分之一的度量，则无论我们执行了多少额外的测序，后续对数倍更改的准确性都将永远不会提高。
Downsampling and log-transforming
在极少数情况下，出于由A.Lun所描述的影响，不适合直接对计数进行缩放。 简而言之，这是由于对数归一化计数的平均值与对数变换后的归一化计数的平均值不同而造成的。 它们之间的差异取决于原始计数的均值和方差，因此相对于计数大小，对数计数的平均值存在系统的趋势。 这通常表现为即使在文库大小归一化之后，轨迹也与文库大小密切相关，如图7.5所示，通过合并和拆分方法生成的合成scRNA-seq数据如图5所示。
# TODO: move to scRNAseq.
library(BiocFileCache)
bfc <- BiocFileCache(ask=FALSE)
qcdata <- bfcrpath(bfc, "https://github.com/LuyiTian/CellBench_data/blob/master/data/mRNAmix_qc.RData?raw=true")
env <- new.env()
load(qcdata, envir=env)
sce.8qc <- env$sce8_qc
# Library size normalization and log-transformation.
sce.8qc <- logNormCounts(sce.8qc)
sce.8qc <- runPCA(sce.8qc)
gridExtra::grid.arrange(
plotPCA(sce.8qc, colour_by=I(factor(sce.8qc$mix))),
plotPCA(sce.8qc, colour_by=I(librarySizeFactors(sce.8qc))),
ncol=2)
SORT-seq CellBench数据中所有合并库和拆分库的PCA图，根据具有库大小派生的大小因子的对数归一化表达值计算得出。 每个点都代表一个库，并通过构造它的混合比例（左）或大小因子（右）来着色
由于问题是由于计数大小的差异而引起的，因此最直接的解决方案是降低取样高覆盖率细胞的以匹配低覆盖率细胞。 这使用大小因子来确定达到大小因子的第1个百分位数所需的每个细胞的减采样。 （只有少数几个具有较小尺寸因子的细胞被简单地按比例放大。我们不会尝试将采样缩减为最小尺寸因子，因为这将导致一个尺寸因子非常低的异常细胞过度丢失信息。）我们可以看到 这消除了前两个PC中与库大小因子相关的轨迹，从而提高了基于混合比的已知差异的分辨率（图7.6）。 对数转换仍然是必需的，但是当细胞之间的计数大小相似时，不再会导致均值变化。
sce.8qc2 <- logNormCounts(sce.8qc, downsample=TRUE)
sce.8qc2 <- runPCA(sce.8qc2)
gridExtra::grid.arrange(
plotPCA(sce.8qc2, colour_by=I(factor(sce.8qc2$mix))),
plotPCA(sce.8qc2, colour_by=I(librarySizeFactors(sce.8qc2))),
ncol=2
)
SORT-seq CellBench数据中库库和拆分库的PCA图，根据对库大小因子成比例缩减采样后从对数转换后的计数计算得出。 每个点都代表一个库，并通过用于构建它的混合比例（左）或大小因子（右）进行着色。
虽然减采样是一种方便的解决方案，但由于需要增加高覆盖率细胞的噪声以避免与低覆盖率细胞之间的差异，因此它在统计上是无效的。 它也比简单缩放慢。 因此，我们只建议在按比例缩放的初始分析显示与大小因子高度相关的可疑轨迹后再使用此方法。 在这种情况下，通过减采样重新确定轨迹是否是对数转换的伪像是一件简单的事情。
