单细胞分析之标准化处理

书籍原始链接：https://www.bioconductor.org/books/release/OSCA/normalization.html

本篇笔记根据书中介绍的内容，按照自己的实际操作情况，进行复现和重新整理。

为什么要做标准化？要解决什么样的问题？进一步消除细胞之间的系统差异（我们前面做质控分析，其实也涉及到这一点），使得后面在比较细胞之间的表达量的时候，这种差异主要来源于生物条件的不同，而非是一些实验过程中的系统误差。

补充小知识点1：标准化与批次校正的区别

标准化：不管批次结构，仅考虑技术偏差（文库之间cDNA捕获与PCR扩增的技术差异，导致覆盖率的不同）

批次校正：仅在批次之间发生，而且必须同时考虑技术偏差和生物学差异，技术偏差影响基因以相似的方式，或者至少是一种与他们的生物物理特性相关的方式，如长度，GC含量，而生物学差异则高度不可预测。

因此，标准化与批次校正是不同的东西，需要使用不同的计算方法。切忌将两者混为一谈。

在众多的标准化策略中，缩放标准化（scaling normalization）是最常用也是最简单的标准化策略。

补充知识点2：缩放标准化

计算方法：将每一个细胞的所有计数值除以细胞特异性的scaling factor。

基本假设：任何细胞特异性偏倚（例如捕获或扩增效率）都会通过缩放该细胞的预期平均数来同等地影响所有基因。

标准化的结果可以用于下游的聚类与降维的分析。

1、加载测试数据集

library(scRNAseq)
example_sce <- ZeiselBrainData()

example_sce数据集的详细信息：

example_sce

 

class: SingleCellExperiment 
dim: 20006 3005 
metadata(0):
assays(1): counts
rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
rowData names(1): featureType
colnames(3005): 1772071015_C02 1772071017_G12 ... 1772066098_A12 1772058148_F03
colData names(22): tissue group # ... altexps_repeat_percent total
reducedDimNames(0):
altExpNames(2): ERCC repeat

2、方案一：文库大小的标准化

文库大小：每一个细胞所有基因表达量的总和。

每个细胞的“library size factor”直接与其库大小成正比，其中定义了比例常数，以便所有单元的平均大小因子等于1。此种标准化可以确保标准化后的表达值与原始数据的范围相同。

library(scater)
#逐个计算3005个细胞的factor，每个细胞有一个对应的factor
lib.sf.zeisel <- librarySizeFactors(example_sce)  
summary(lib.sf.zeisel)

 

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1721  0.5437  0.8635  1.0000  1.2895  4.2466 

在这个数据集中，最大的factor值（4.2466）与最小的factor值（0.1721）之间差40倍之多。这是scRNA-seq数据覆盖率变异性的典型特征。

绘制直方图，统计3005个细胞的library size factor的基本分布情况。

hist(log10(lib.sf.zeisel), xlab="Log10[Size factor]", col='lightblue')
#注意加log10()处理之后，使图像整体更加美观

每一个细胞的library size factor的具体分布
（我感觉通过上述步骤只是计算了一个factor，并没有把数据矩阵中的每一个数据都最终标准化。是做的还不够的）

这种归一化方法的限制条件：这种方法在使用的过程中，存在一个基本的假设，即两个细胞内的上调的基因集与其下调的基因集幅度相同，可以互相抵消（称之为平衡）。但是实际的情况是，这种情况在scRNA-seq的数据中几乎不存在，因此，这种归一化方法不能为下游的分析提供可靠的归一化数据。

由于，归一化的可靠性，并不是探索scRNA-seq数据的主要考虑的因素（所以，在某些情况下，可以睁一只眼闭一只眼？）。在聚类的过程中，这种方法所引起的归一化偏差不会影响簇的分离，而只会影响细胞类型或者簇之间的log10()值的幅度或者方向（方向？）。

因此，这种归一化方法，适用于下游的聚类分析。

3、方案二：去卷积标准化

我们对上文提到的成分偏见（composition biases）这一现象进行重新的分析。假设有两个细胞，细胞A中的geneX相较于细胞B而言是上调的。

（1）这种上调意味着，在A细胞中更多的资源被分配给了geneX，而其他非差异表达基因的覆盖率则相应的降低（当每个细胞中测序的总counts数是固定的）。

（2）当geneX被分配有更多的reads或者UMIs,A细胞的测序总counts是上升的。增加了上述方法中的library size factor的大小，而使得所有正常表达的非差异表达基因的值减少。

在这以上的两种情况下，A细胞中的非差异表达的基因，相较于B细胞而言，会被错误的认为下调了（也就是存在一个超异常值，从而拉高了整体的平均值，在归一化的过程中，一般值相对平均水平低了，如果除去异常值，则是不一样的效果）。

这时，一种标准化的策略是，汇总了许多细胞中的counts以增加counts的大小（单细胞数据中普遍的0或者很低的值，所以总counts很小？），以进行准确的size factor估算。然后，将基于库的size factor“分解”为cell-based factors，以标准化每个细胞的表达谱。（这里我是没怎么明白的？许多撮面疙瘩，揉揉，变成一个面团，然后再重新搓成许多的面疙瘩？）

BiocManager::install("scran")
library(scran)
set.seed(100)

#我们使用了一个预集群步骤quickCluster()，具体的原理下面会介绍。
clust.zeisel <- quickCluster(example_sce) 

#查看cluster的分群
table(clust.zeisel)

#按照分群，计算factor
deconv.sf.zeisel <- calculateSumFactors(example_sce, cluster=clust.zeisel)
summary(deconv.sf.zeisel)

clust.zeisel的具体的数值：

clust.zeisel
  1   2   3   4   5   6   7   8   9  10  11  12  13  14 
139 164 251 234 372 272 151 428 117 148 289 225 111 104 

deconv.sf.zeisel的具体数值（可以与上面的文库大小的factor进行比较，同一套数据）：

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1209  0.4621  0.8082  1.0000  1.3177  4.8508 

• 使用quickCluster()进行聚类，每一个类进行分组的标准化。重新调整分组因子，以便群集之间具有可比性。

比较文库大小标准化得到的每一个细胞的factor与去卷积化得到的factor之间的关系。

plot(lib.sf.zeisel, deconv.sf.zeisel, xlab="Library size factor",
     ylab="Deconvolution size factor", pch=16,
     col=as.integer(factor(example_sce$tissue)))
abline(a=0, b=1, col="green")

绿色的那条线为y=x的渐近线。这张图展示的是两组标准化得到的factor的比较。

 

标准化的方法，对于差异表达基因的分析影响还是挺大的，所以在操作的时候要尽量回避能够产生composition biases的策略（如文库大小标准化策略）。

4、方案三、通过spike-in标准化（在本节将重点介绍spike-in的用途）

基于的基本假设：每一个细胞中添加有等量的spike-inRNA。spike-in转录本的覆盖率的系统差异（systematic difference）只能归因于细胞特异性偏差，如捕获的效率，测序的深度。

为了移除这些偏差，我们通过“spike-in size factor”来均衡细胞之间的spike-in覆盖率。相较于之前的标准化的方法，spike-in标准化方法对系统的生物学条件无要求（比如，缺少许多差异表达基因）。

在实际的处理条件中，spike-in标准化通常被用在如果每一个细胞的总RNA含量的差异是我们感兴趣的内容，并且这必需保留在下游的分析中。对于给定的细胞，其内源RNA总量的上升不会引起其spike-in size factor的变化。这确保着，总RNA含量对整个细胞群表达的影响，不会在scaling的过程中除去。相比之下，上述其他标准化方法只是将总RNA含量的任何变化解释为偏差的一部分，并将其消除。

#加载scran包,librarySizeFactors()函数位于这个包中，如果之前已经加载过不必重复加载了
library(scran)
#计算spike-in factor
spike.sf.zeisel<-librarySizeFactors(altExp(example_sce,"ERCC"))
summary(spike.sf.zeisel)

关于函数computeSpikeFactors()，librarySizeFactor(altExp())等函数的具体用法，可参考链接。

对于同样一套数据集，spike-in factor的值分布情况。

 Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1565  0.7697  0.9844  1.0000  1.1808  3.4171 

5、应用size factors。

（1）缩放以及对数转换

补充知识点3：为什么要进行log转换？

对数值的变化，代表着表达值的对数倍的变化。这对于下游的基于欧氏距离的处理步骤如聚类和降维，非常必要。通过对对数转换后的数据进行操作，我们确保这些过程是基于细胞表达的对数倍来计算距离。

在进行对数转换时，我们通常会添加一个伪计数以避免未定义的零值。对于低丰度基因，较大的伪计数将有效地将细胞之间的对数倍变化缩小至零，这意味着下游的高维分析将更多地由高丰度基因的表达差异来驱动。相反，较小的伪计数将增加低丰度基因的相对贡献。常见的做法是使用1作为伪计数，原因很简单，即它保留了原始矩阵中的稀疏性（即输入中的零在变换后仍为零）。除大多数病理情况外，此方法在所有情况下均适用。

set.seed(100)
clust.zeisel <- quickCluster(example_sce) 
#这里和前面计算factor的操作不太一样的是，将其赋值给了原始数据集
example_sce<- computeSumFactors(example_sce, cluster=clust.zeisel)
example_sce<- logNormCounts(example_sce)
assayNames(example_sce)

 查看结果。

 colnames(colData(example_sce))

下面是colnames()的结果，可以看到sizeFactor被添加到了数据集中。

[1] "tissue"         "group #"        "total mRNA mol" "well"           "sex"           
 [6] "age"            "diameter"       "cell_id"        "level1class"    "level2class"   
[11] "sizeFactor"   

6、还有一种情况，没怎么看懂（这里不写了）。

 

后记：

其实按照操作流程一步一步将代码跑出来，倒不是难的。

我觉得这里非常难的地方，在于理解为什么要这样去处理？应该如何选择适合的方法？这个方法是否适合于我的数据？处理之后会对我的数据下游的分析产生什么样的影响？这个影响是我们想要的吗？

我们要从中做出选择，并尽可能做出优化。很多函数是直接内嵌在一体化操作的R包里，运行，可以直接得到结果。但是，你需要明白，这个函数内部进行着什么操作？是我们想要的吗？