OSCA单细胞数据分析笔记6—Normalization

对应原版教程第7章http://bioconductor.org/books/release/OSCA/overview.html 标准化是在剔除不合格细胞之后，尽可能消除细胞文库间大小的差异性，从而得到准确、有意义的分析结果。一起来学习吧~

笔记要点

• 1、背景知识
  • 1.1 为什么要标准化
  • 1.2 标准化≠批次校正
• 2、标准化最常见流程
  • 2.1 假设
  • 2.2 计算细胞文库size factor
  • 2.3 log转换
• 3、标准化特殊情形
  • 3.1 细胞文库标准化之平衡的差异基因？
  • 3.2 根据外参转录本的标准化

1、背景知识

1.1 为什么要标准化

• 因为在制备单细胞文库时，排除低质量细胞的前提下，最理想的测序结果就是每个细胞的文库大小均相同，这样不同细胞的相同基因水平才具有可比性。
• 但由于客观原因(cDNA捕获、PCR复制)等技术误差的存在，细胞间的文库大小会存在一定差异。造成的后果就是某一基因在细胞间的表达差异就由技术误差和生物水平差异构成。从而可能导致错误的结论。
• 标准化这一步简单来说就是消除细胞文库大小的差异性。

1.2 标准化≠批次校正

• 标准化的前提是认为细胞间的文库差异性完全由测序时的技术误差造成的。而这种技术误差总体是可消除的。因为我们常常假设技术误差对所有的基因表达有相同方向/模式的影响。例如技术误差使基因A的表达量降低10%，那么其它所有基因的表达量都会降低10%。
• 批次效应则不仅仅意味着批次之间有更为明显的技术误差存在，而且各个批次的取样细胞的生理状态不同所造成的额外生物水平误差往往是难以预测的，因此需要采取不同于标准化的思路进行校正。
• 关于批次效应的校正方法，在之后章节会有说明。这里只是想强调下标准化与批次校正是两个不同的概念。

2、标准化最常见流程

2.1 假设

• 每一个细胞测序时受到的的技术误差，对其内的所有基因表达有相同模式的影响。
• 用一个指标(size factor)评价每一个细胞受到的技术误差大小、方向(相对所有细胞文库的均值，而不是绝对意义上的偏离标准文库的距离)。然后分别根据每个细胞的指标的值进行对该细胞的标准化。最后使基因在不同细胞间的表达具有可比性。可能说的比较抽象，下面举例子应该会好理解些。

2.2 计算细胞文库size factor

sce.zeisel
# dim: 19839 2816 
# assays(1): counts

(1)首先计算每个细胞文库的大小

lib.size=apply(counts(sce.zeisel), 2,sum)
head(lib.size)
#1772071015_C02 1772071017_G12 1772071017_A05 1772071014_B06 1772067065_H06 1772071017_E02 
#         22354          22869          32594          33525          21694          25919

(2)计算所有细胞文库均值【参考标准】

mean(lib.size)
#15550.11

(3)计算每个细胞的细胞文库size.factor

lib.sf=lib.size/mean(lib.size)
head(lib.sf)
# 1772071015_C02 1772071017_G12 1772071017_A05 1772071014_B06 1772067065_H06 1772071017_E02 
# 1.437546       1.470664       2.096062       2.155933       1.395102       1.666805

如上结果表示相对每个细胞文库相对于平均细胞文库大小的比例。

library(scater) #一步计算
lib.sf.zeisel <- librarySizeFactors(sce.zeisel)
summary(lib.sf.zeisel)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.1757  0.5680  0.8680  1.0000  1.2783  4.0839
#均值一定为0
#值大于1的需要缩小；小于1的需要放大

(4)然后根据每个细胞的细胞文库size factor对该细胞的所有基因表达进行归一化处理，即基因表达量除以细胞文库因子。这里的“一”就代表平均细胞文库大小。

#以第一个细胞为例
normalized <- counts(sce.zeisel[,1])/lib.sf[1]

2.3 log转换

• log转换最主要的目的是我们在之后计算细胞的距离时主要想依据细胞间基因表达差异倍数变化而不是简单的差值。试想细胞1的基因A表达量为50，B的表达量为1000；细胞2的基因A表达量为10，B的表达量1000。那个基因变化最明显的，显然使用差异倍数更有意义。而取log之后的减法即为除法，可以表示log FoldChange
• 由于单细胞测序结果的稀疏性，很多基因的表达值为0，无法取log。所以会为所有细胞的基因加一个pseudo-count表达(一般取1，log后还是为0，未改变稀疏性)。

#一般取log2
log_normalized=log2(normalized + 1)

如上就得到了最终标准化后的表达矩阵。

2.4 logNormCounts一步完成

sce.zeisel <- logNormCounts(sce.zeisel)
assays(sce.zeisel)
#List of length 2
#names(2): counts logcounts
counts(sce.zeisel)[1:4,1:4]
logcounts(sce.zeisel)[1:4,1:4]

如上结果，标准化后的结果保存在logcounts assay里面

3、标准化特殊情形

3.1 细胞文库标准化之平衡的差异基因？

(1)在使用文库标准化size factor时，有一个潜在的假设就是balanced differentially expressed genes,在某一细胞中若存在某些基因表达上调，那么一定会存在某些基因下调。

• 举例1：假设一个细胞测序表达有100个基因，文库大小为100。基因2-99表达量均为1，而基因1表达量为1.5；而基因100的表达量为0.5(对照细胞100个基因表达全为1)。就可以说是balanced DE，符合我们预期的假设。因为不影响剩下的98个non-DEG

(2)但是单细胞测序中，往往会出现imbalanced DEG，即比较多的出现上调基因或者下调基因，从而使原本的non-DEG变成DEG。术语称之为composition bias

• 举例2：假设一个细胞测序表达有100个基因，文库大小为101。基因2-100表达量均为1，而基因1表达量为2(对照细胞100个基因表达全为1)。
• 举例3：假设一个细胞测序表达有100个基因，文库大小为100。基因1表达量为2，基因2-100表达量为0.99，而(对照细胞100个基因表达全为1)。
• 无论是例2，还是例3，在经过标准化之后的差异分析结果就是基因1真实相对上调；基因2-99表面相对下调，其实本质为non-DEG。

(3) 从对之后的分析影响来看，作者认为composition bias对于单细胞之后的聚类分群、Top marker gene结影响不会很大。但如果想进行单基因水平的分析，还是最好消除这种误差。

(4) 如何最大化避免composition bias

• 对于传统的Bulk RNA-seq数据，DESeq2包的estimateSizeFactorsFromMatrix()函数、edgeR包的calcNormFactors()函数会考虑到这种误差。它们主要是假设大部分的基因均为non-DEG，所以当一个测序结果里的大部分基因出现相对于对照组的系统偏移时，会认为是composition bias，计算合适的size factor予以去除。
• 对于单细胞测序结果来说，因为表达矩阵的稀疏性，上述的函数方法不太适合。可以设置calculateSumFactors()函数的cluster参数。主要是假设每个cluster间为non-DE，先对cluster总体的pool-based size factor的标准化；再对每个cluster里的细胞分别进行cell-based size factor标准化。

library(scran)
set.seed(100)
clust.zeisel <- quickCluster(sce.zeisel) 
table(clust.zeisel)
deconv.sf.zeisel <- calculateSumFactors(sce.zeisel, cluster=clust.zeisel)

sce.zeisel <- computeSumFactors(sce.zeisel, cluster=clust.zeisel, min.mean=0.1)
sce.zeisel <- logNormCounts(sce.zeisel)

3.2 根据外参转录本的标准化

（1）前提假设

• 每个细胞初始水平添加的spike-in转录本的含量是一致的；
• 测序技术误差无论对于内源基因还是外参转录本都有相同的模式的表达影响。

（2）适用情形

• 之前提到的细胞文库size factor标准化相当于是把所有细胞的文库大小拉到同一水平进行基因的比较；
• 但是如果观察的目的之一就是细胞内基因表达总量的变化，例如T细胞激活前后等，具体要结合生物背景知识了；
• 相关函数有computeSpikeFactors()，示例在原教程里有，这里就不多展开了。