对细胞文库差异进行normalization

豆豆写于2019.9.23  花花已经好几天没写了，她说麦旋风可以召唤她！

• 如何对低质量细胞进行质控检验

• 关于细胞周期推断的知识更新

• scRNA表达矩阵的基因层面检查

再次借助https://github.com/satijalab/seurat/issues/950的一个回答：

• Normalization “normalizes” within the cell for the difference in sequenicng depth / mRNA thruput

• Scaling “normalizes” across the sample for differences in range of variation of expression of genes

关于它们的中文翻译实在是错综复杂，记不清楚就简单记，反正没有正确答案，自己知道这一步做了什么最重要。normalization一般是对文库处理，目的消除一些技术差异；scale一般对基因表达量处理（典型的z-score：表达量减均值再除以标准差），目的是后续分析不受极值影响（比如做热图我们不关心哪个基因表达量最多，而关心哪个基因在样本间最有差异，虽然这个基因不一定表达量很高。比如A基因在样本1中表达量为900，在样本B中表达量为901；而B基因在样本1中表达量为10，在样本B中表达量为50。我们更关心B基因。但直接作图，往往会因为A基因的绝对表达量很高，而将B基因在图中掩盖掉。scale这个操作就是帮我们”拨开云雾见青天“）

1 去卷积方法（deconvolution）根据内源基因计算size factor

表达量的差异除了生物学因素外，还会受到细胞捕获效率、测序深度的影响(Stegle, Teichmann, and Marioni 2015)。在进行下游定量分析之前，需要利用归一化去掉这些“无关因素”产生的细胞差异。

假设大部分基因在细胞间并非差异表达，统计思想是：any systematic difference in count size across the non-DE majority of genes between two cells is assumed to represent bias and is removed by scaling.

更具体说就是计算“size factors” （represent the extent to which counts should be scaled in each library）

size factor的计算方法有多种，例如常规转录组中DESeq2的estimateSizeFactorsFromMatrix 、edgeR的calcNormFactors 。但是这里单细胞数据又和常规转录组不同，它存在了太多的低丰度甚至是0表达（从上面的过滤结果也能感知），一个细胞一个细胞单独进行处理应该是误差较大。

于是开发了一种基于去卷积的新方法：先将许多细胞的counts混合起来，来增加总体的counts值，提高了后续计算size factor的准确度 (Lun, Bach, and Marioni 2016)。计算得到的混合size factor（Pool-based size factors），再“去卷积”，将pool size factor还原给每个细胞的factor。

说起来复杂，操作简单：

sce <- computeSumFactors(sce)
summary(sizeFactors(sce))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3464  0.7117  0.9098  1.0000  1.1396  3.5695

如何知道计算的size factor准确性好不好呢？

做一个散点图，看看size factor与lib size的关系

plot(sce$total_counts/1e6, sizeFactors(sce), log="xy",
    xlab="Library size (millions)", ylab="Size factor",
    col=c("red", "black")[sce$Oncogene], pch=16)
# 右下角添加图例，对数据中两种细胞类型进行了区分
legend("bottomright", col=c("red", "black"), pch=16, cex=1.2,
    legend=levels(sce$Oncogene))

图中可以看到，size factor与文库大小相关性很好，说明大部分文库大小差异都是由于技术因素（捕获效率、测序深度）造成的。

【问：为何size factor不考虑生物因素？
答：因为这个名词创造之前，就已经做了无差异基因的假设条件，它只体现一些技术误差导致的差异，所以只有存在技术误差时，lib size才和size factor的相关趋势接近直线】

另外，从这个图中，看到红点（原癌基因诱导细胞）比黑点（对照组）的size factor要小，这从生物学角度也能解释通：黑点红点总体lib size一样，红点诱导后基因表达量升高（即生物学因素导致的文库变化升高），那么红点非生物学因素导致的文库变化就要降低

Bioconductor作者的观点

• 虽然去卷积法对高频次的0表达量处理效果不错，但如果太多的0表达量最后还是会失败。如果计算的size factor为负数，这很明显没有意义。因此computeSumFactors函数在计算size factor前，会自动去先按平均表达量去除低丰度转录本，一般会以min.mean=1作为过滤条件

• 在标准化之前，细胞应该经常查看一下质控结果，保证去掉那些只有很少基因表达的细胞。否则，`computeSumFactors`就会由于低质量细胞产生负数。因为size factor的计算涉及到了混合细胞的操作，最后操作的是pool library，因此每个细胞都有太多的0，最后校正也很难去掉它们

• computeSumFactors()函数有一个参数`sizes`，它指定每个混池使用多少细胞计算size factor。每个混池使用的细胞数量越多，越准确，但也会消耗计算时间和内存。一般来讲，一个混池至少要用20个以上细胞才能保证有足够的非0表达量去计算。

• 对于高度异质性的数据集来说，因为要保证”没有差异基因“这个原假设。所以可以先利用quickCluster()做一个粗略的聚类，然后将结果通过clusters参数传递给computeSumFactors()再计算

2 根据spike-in转录本计算size factor

根据内源基因count结果计算的size factor是不适用于spike-in的，原因如下：

一般单细胞实验都不是定量实验（也就是说，我们事先不能规定好所有的细胞文库大小都是一样的），我们常常得到的数据是：包含更多转录本的细胞具有更大的文库，因此计算的size factor也更大，目的就是消除这个文库差异的影响。
这一点也很好理解，但是基于这个事实，我们人工在不同细胞中加入的spike-in转录本浓度却是一样的，在小文库细胞中某个ERCC表达量是100，换到另一个大文库细胞中，这个ERCC的表达量也只是近似，不会呈现数量级的差别。再用原来的方法对一个稳定的表达量进行校正，结果会”矫枉过正“。

另外，即使使用了定量的实验，每个细胞的cDNA文库大小都是不变的，那么其中如果内源基因表达量升高，那么势必会导致外源spike-in表达量下降，它们的变化趋势不同，因此利用内源基因计算的size factor也不能直接利用到spike-in身上。

虽然原来方法不适用，但内源基因们都校正过了，spike-in也不能不校正。于是有一个单独为它开发的方法：假设所有的spike-in转录本表达无差异（当然这个假设是符合客观事实的）。For each cell, the spike-in-specific size factor is defined as the total count across all transcripts in the spike-in set

利用一个函数：

sce <- computeSpikeFactors(sce, type="ERCC", general.use=FALSE)

利用参数general.use=FALSE 将结果单独存放在sce对象中，保证了计算的size factor只对ERCC使用，而对内源基因生效

3 将size factor应用到normalization这一步

函数很简单：

sce <- normalize(sce)

它做的事情就是：每个细胞中基因的表达量除以这个细胞计算的size factor，其中内源基因除以computeSumFactors()得到的factor，spike-in除以computeSpikeFactors()得到的factor。然后做个log2( +1)处理，降低数据的维度又不破坏数据原始结构（也就是说，原来大的数log完了还是大，但是在坐标轴上的数字就能减小，更方便展示数据全貌）

最后计算结果保存在：logcounts中，可以用logcounts(sce)调用表达矩阵