标准归一，一对CP

豆豆写于19.8.19
归一化和标准化是两个常用的数据变换需求，但它们常常被混淆，我也是经常分不清楚，因此有必要记录一下它们的区别

数据归一化（normalize）

为何normalize要叫“归一化”？
这个名称不必纠结，看它使用的方法应该就好理解：因为我们经常关心的就是看差异表达基因，也就是找一个基因在不同样本的差异，而这种差异，往往不能直接进行比较，因为不同样本的测序深度和文库大小可能不能，那么基因表达量差异就存在一部分因素来自这里。为了更真实反映基因差异的生物学意义，就要将数据进行一个转换（例如RPKM、TPM等），可以让同一基因在不同样本中具有可比性。

“归一”归的就是这个文库大小。

另外呢，原始表达量是一个离散程度很高的值，有的高表达基因表达量可能成千上万，可有的却只有几十。我们即使采用了一些归一化的算法消除了一些技术性误差，但真实存在的表达量极差往往会影响之后绘制热图、箱线图的美观性，于是可以另外采用log 进行一个区间缩放（却不会影响真实值），比如原来表达量为1的，用log2(1+1)转换后结果依然为1；原来表达量为1000的，log2(1000+1)转换后约为10。那么原来相差1000倍的变化，现在只差10倍，在不破坏原有数据差异的同时，使数据更加集中。

关于Seurat归一化原理，可以看这一篇：https://www.biorxiv.org/content/biorxiv/early/2019/03/18/576827.full.pdf

当过滤掉一部分细胞以后，就要会数据进行归一化。默认是进行一个全局的LogNormalize操作，具体方法是：

normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result
翻译成公式就是：log1p(value/colSums[cell-idx] *scale_factor) ，其中log1p指的是log(x + 1)，(Tip:只看到log就表示以e为底的log值，用log1p(exp(1)-1)测试一下，看结果是1 )

这里也有解释：https://github.com/satijalab/seurat/issues/833

当然，除了这一种默认的算法，还有：

• CLR: Applies a centered log ratio transformation

• RC: Relative counts. Feature counts for each cell are divided by the total counts for that cell and multiplied by the scale.factor. No log-transformation is applied. For counts per million (CPM) set scale.factor = 1e6

得到的结果存储在：pbmc[["RNA"]]@data

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
# 当然其中都是默认参数，于是直接写这种也是可以的
pbmc <- NormalizeData(pbmc)

这一步要因人而异，可以看：https://bioinformatics.stackexchange.com/questions/5115/seurat-with-normalized-count-matrix 假入有一个TPM的count矩阵，那么就可以不需要使用Seurat::NormalizeData()操作了，因为TPM已经根据测序深度进行了归一化，只需要进行log降一下维度即可。如果后续进行ScaleData操作，程序会检测是否使用了Seurat提供的标准化方法，如果我们使用自己的的归一化数据，那么就可能出现一个warning提醒，不过到时候不想被提醒，可以设置check.for.norm =F

数据标准化（scale）

它的目的是实现数据的线性转换(scaling)，这也是降维处理(如PCA)之前的标准预处理。
Seurat包主要利用了ScaleData函数，回忆一下，之前归一化（normalize）这里做的操作是log处理，它是对所有基因的表达量统一对待的，最后放在了一个起跑线上。但是为了真正找到差异，我们还要基于这个起跑线，考虑不同样本对表达量的影响

它做了：

• 将每个基因在所有细胞的表达量进行调整，使得每个基因在所有细胞的表达量均值为0

• 将每个基因的表达量进行标准化，让每个基因在所有细胞中的表达量**方差为1

# 先使用全部基因
all.genes <- rownames(pbmc)
> length(all.genes)
[1] 13714
pbmc <- ScaleData(pbmc, features = all.genes)
# 结果存储在pbmc[["RNA"]]@scale.data中

ScaleData的操作过程是怎样的?

看一下帮助文档就能了解：
ScaleData这个函数有两个默认参数：do.scale = TRUE和do.center = TRUE，然后需要输入进行scale/center的基因名（默认是取前面FindVariableFeatures的结果）。
scale和center这两个默认参数应该不陌生了，center就是对每个基因在不同细胞的表达量都减去均值；
scale就是对每个进行center后的值再除以标准差（就是进行了一个z-score的操作）

再看一个来自BioStar的讨论：

(https://www.biostars.org/p/349824/)

问题1：为什么在`NormalizeData() `之后，还需要进行`ScaleData`操作？

前面NormalizeData进行的归一化主要考虑了测序深度/文库大小的影响（reads*10000 divide by total reads, and then log），但实际上归一化的结果还是存在基因表达量分布区间不同的问题；而ScaleData做的就是在归一化的基础上再添zero-centres （mean/sd =》 z-score），将各个表达量放在了同一个范围中，真正实现了”表达量的同一起跑线“。

另外，新版的ScaleData函数还支持设定回归参数，如（nUMI、percent.mito），来校正一些技术噪音

问题2： `FindVariableGenes() or RunPCA() or FindCluster()`这些参数是基于归一化`Normalized_Data`还是标准化`Scaled_Data` ？

所有操作都是基于标准化的数据Scaled_Data ，因为这个数据是针对基因间比较的。例如有两组基因表达量如下:

g1 10 20 30 40 50
g2 20 40 60 80 100

虽然它们看起来是g2的表达量是g1的两倍，但真要降维聚类时，就要看scale的结果：

> scale(c(10,20,30,40,50))
           [,1]
[1,] -1.2649111
[2,] -0.6324555
[3,]  0.0000000
[4,]  0.6324555
[5,]  1.2649111
attr(,"scaled:center")
[1] 30
attr(,"scaled:scale")
[1] 15.81139
> scale(c(20,40,60,80,100))
           [,1]
[1,] -1.2649111
[2,] -0.6324555
[3,]  0.0000000
[4,]  0.6324555
[5,]  1.2649111
attr(,"scaled:center")
[1] 60
attr(,"scaled:scale")
[1] 31.62278

二者标准化后的分布形态是一样的（只是中位数不同），而我们后来聚类也是看数据的分布，分布相似聚在一起。所以归一化差两倍的两个基因，根据标准化的结果依然聚在了一起

问题3： `ScaleData() `在scRNA分析中一定要用吗？

实际上标准化这个过程不是单细胞分析固有的，在很多机器学习、降维算法中比较不同feature的距离时经常使用。当不进行标准化时，有时feature之间巨大的差异（就像👆问题2中差两倍的基因表达量，实际上可能差的倍数会更多）会影响分析结果，让我们误以为它们的数据分布相差很远。

很多时候，我们会混淆normalization和scale：那么看一句英文解释：

Normalization “normalizes” within the cell for the difference in sequenicng depth / mRNA throughput. 主要着眼于样本的文库大小差异
Scaling “normalizes” across the sample for differences in range of variation of expression of genes . 主要着眼于基因的表达分布差异

另外，看scale()函数的帮助文档就会发现，它实际上是对列进行操作：

那么具体操作中是要先对表达矩阵转置，然后scale，最后转回来

t(scale(t(dat[genes,])))

问题4：`RunCCA() `函数需要归一化数据还是标准化数据？

通过看这个函数的帮助文档：

RunCCA(object, object2, group1, group2, group.by, num.cc = 20, genes.use,
  scale.data = TRUE, rescale.groups = FALSE, ...)

显示scale.data = TRUE，那么就需要标准化后的数据

总结一下

• normalize是归一化，耳熟能详的TPM、FPRM就是做这个事；还有常见的log2、log10处理；

• scale是标准化（standardization），主要采用z-score，就是表达量减均值再除以标准差

• z-score其实是没有具体含义的，它的作用就是反映表达量偏离均值多少个sd值，数值越大表示越离散，对于一个符合正态分布的数据来说，有95%的数据会落在mean -2sd ~ mean+2sd中 ，因此一般做热图时将z-score的最大值设为2，最小值设为-2，也是能代表大多数据的，并且可以不受特别离群的点影响，保证主体的差异

• scale包含中心化，中心化就是表达量减均值

• 降维、聚类之前要进行标准化（z-score）

• 生信中归一化和标准化也没有明确的区别，可以认为归一化是特殊的标准化，所以经常听到的各种标准化算法是做了归一化这件事