年底复盘，发现很多知识点理解并不深入，想对这些不清晰的知识点重新梳理一下。
一、NormalizeData()归一化
考虑到文库测序深度的影响，我们需要对单细胞counts矩阵数据进行归一化处理。NormalizeData()默认方式是LogNormalize， 其他方法有CLR，RC，详细说明详见链接。
- LogNormalize: 每个细胞的特征计数除以该细胞的总计数并乘以scale.factor。 然后使用log1p对其进行自然对数转换；
- CLR (centered log ratio transformation): 应用中心对数比转换；
- RC: 相对计数。 每个细胞的特征计数除以该细胞的总计数并乘以scale.factor，不应用对数转换(log-transformation)。 对于每百万计数 (CPM) 设置 scale.factor = 1e6
LogNormalize的做法是，每个细胞，每个基因的count除以该细胞的总counts，乘以scale.factor（默认是 10,000，就好像所有单细胞总共有10k UMI），并对获得的值进行 log1p转换，进行自然对数转换。归一化的数据存储在“RNA” assay的 seurat_obj[['RNA']]@data中。
示例：
data("pbmc_small")
pbmc_small
pmbc_small <- NormalizeData(object = pbmc_small)
1.1 NormalizeData()源码解析
本想学习下seurat的NormalizeData()，发现并没有那么轻松，追根溯源，NormalizeData()通过Rcpp方式调用C++函数LogNorm，加快运行速度。
LogNorm <- function(data, scale_factor, display_progress = TRUE) {
.Call('_Seurat_LogNorm', PACKAGE = 'Seurat', data, scale_factor, display_progress)
}
// [[Rcpp::export(rng = false)]]
Eigen::SparseMatrix<double> LogNorm(Eigen::SparseMatrix<double> data, int scale_factor, bool display_progress = true){
Progress p(data.outerSize(), display_progress);
Eigen::VectorXd colSums = data.transpose() * Eigen::VectorXd::Ones(data.rows());
for (int k=0; k < data.outerSize(); ++k){
p.increment();
for (Eigen::SparseMatrix<double>::InnerIterator it(data, k); it; ++it){
it.valueRef() = log1p(double(it.value()) / colSums[k] * scale_factor);
}
}
return data;
}
看这段C++代码，理解还是不够深入。通过一个小案例，我们看下LogNormalize()具体的运算。
LogNormalize(data, scale.factor = 10000, verbose = TRUE, ...)
示例：
step1：构造一个细胞基因表达矩阵5*5，列为细胞，行为基因；如cell1的gene1的count为2；cell1的总counts数目为6(2+3+1);
step2: 我们计算cell1的gene1归一化之后的数值：log1p(2* 1e4/6)=8.112028;
step3: 对所有细胞的每个基因归一化，log1p(value/colSums[cell-idx] *scale_factor)，
R语言执行：log1p(sweep(mat, 2, Matrix::colSums(mat), FUN = "/") * 1e4)
mat <- matrix(data = rbinom(n = 25, size = 5, prob = 0.2), nrow = 5)
rownames(mat) <- paste0("gene", c(1:5))
colnames(mat) <- paste0("cell", c(1:5))
mat
# cell1 cell2 cell3 cell4 cell5
# gene1 2 2 0 3 0
# gene2 0 1 0 1 1
# gene3 0 1 0 0 1
# gene4 3 0 0 3 0
# gene5 1 1 1 3 0
# 运用LogNormalize函数
mat_norm <- LogNormalize(data = mat)
mat_norm
# 5 x 5 sparse Matrix of class "dgCMatrix"
# cell1 cell2 cell3 cell4 cell5
# gene1 8.112028 8.294300 . 8.006701 .
# gene2 . 7.601402 . 6.908755 8.517393
# gene3 . 7.601402 . . 8.517393
# gene4 8.517393 . . 8.006701 .
# gene5 7.419181 7.601402 9.21044 8.006701 .
# 手动执行
mat_norm.r <- log1p(sweep(mat, 2, Matrix::colSums(mat), FUN = "/") * 1e4)
mat_norm.r
# cell1 cell2 cell3 cell4 cell5
# gene1 8.112028 8.294300 0.00000 8.006701 0.000000
# gene2 0.000000 7.601402 0.00000 6.908755 8.517393
# gene3 0.000000 7.601402 0.00000 0.000000 8.517393
# gene4 8.517393 0.000000 0.00000 8.006701 0.000000
# gene5 7.419181 7.601402 9.21044 8.006701 0.000000
step4: 我们测试下log1p处理之后的数值进行expm1还原；
cell1的gene1归一化之后进行expm1还原的数值：expm1(8.112028)=3333.333
cell1的gene1的count除以细胞总counts，乘以scale.factor =(2* 1e4/6)=3333.333，和上面的expm1还原数值是一样的；
# 进行expm1转换
expm1(mat_norm.r)
# cell1 cell2 cell3 cell4 cell5
# gene1 3333.333 4000 0 3000 0
# gene2 0.000 2000 0 1000 5000
# gene3 0.000 2000 0 0 5000
# gene4 5000.000 0 0 3000 0
# gene5 1666.667 2000 10000 3000 0
1.2 NormalizeData()相关疑问
问题1：初学单细胞分析的时候，比较疑惑下面两种处理方式，两种计算的归一化数值是否一致?
方式1：
pbmc_list <- list()
pbmc_list[["pbmc10k_3p"]] <- srat_3p
pbmc_list[["pbmc10k_5p"]] <- srat_5p
for (i in 1:length(pbmc_list)) {
pbmc_list[[i]] <- NormalizeData(pbmc_list[[i]], verbose = F)
}
方式2：
pbmc_seurat <- merge(srat_3p,srat_5p)
pbmc_seurat <- NormalizeData(pbmc_seurat, verbose = F)
答案是：两种归一化后的数值是一样的，因为归一化LogNormalize的变量是单个细胞某基因的count除以该细胞的总counts数目，乘以恒量scale.factor，再取 log1p。无论两个样本是否merge，都不影响该细胞的总counts数目，只是对单个细胞的基因count去偏态。
问题2：缩放因子scale.factor为什么默认10^4？
在网上找到关于此问题的解释：
缩放因子为 10^4 的原因很简单：
至少在使用 10X-Genomics 技术时，每个细胞的总RNA计数大约为1000到 50000。如果将一个基因的count（通常为 在 0 到 200 的范围内；尽管大多数count低于 10）通过这样的计数，除以细胞总counts，然后取 log1p，得到的数值很小，小数点后有许多零，在绘图上难以理解。那么使用 10^4 的比例因子，基因的count除以细胞总counts，乘以比例因子10^4 ，可以避免这种极小数字的出现，似乎是合适的。 您可以获得 0.X 到 10 范围内的对数归一化计数，或者我们人类很容易理解的东西。
问题3：怎么理解log1p？
在seurat的源码中，会出现大量关于log1p和expm1的计算。
log1p = log(number+1)，等价的写为ln(1+x)，当 number 的值接近零也能计算出准确结果；
expm1 = exp(number) - 1，等价的写为 ex +1，当 number 的值接近零也能计算出准确结果；
由上面的log1p曲线可以看出，曲线上升而后平缓。log1p( ) 的使用就像是一个数据压缩到了一个区间，与数据的标准类似。其逆运算就是expm1的函数;
由于使用的log1p()对数据进行了压缩，最后需要将预测出的平滑数据进行一个还原，而还原过程就是log1p的逆运算expm1。
如果数据高度偏态，则使用对数变换，有两种处理方式：
1）对数变换 即将原始数据X的对数值作为新的分布数据：
x1 = log(x)
当原始数据中有小值及零时：
x1 = log1p(x)
应用场景
如果数据非正态，有左偏情况，可以使用log1p进行平滑
1.去偏度：
在进行数据预处理的时候，我们对偏度比较大的数据进行log1p(x)处理，即取对数，使得数据在一定程度上符合正态分布（Normal distribution)的特征。正态分布分布也称高斯分布（Gaussian distribution）。此步处理可能会使我们后续的分类结果得到一个好的结果。
不同细胞的起始底物浓度不同，导致cDNA捕获或PCR扩增效率差异。原始表达量是一个离散程度很高的值，有的高表达基因表达量可能成千上万，可有的却只有几十。归一化的目的就是去除细胞间与真实表达量无关的技术因素，方便后续比较。
2.考虑极小值：
log1p函数有它存在的意义，即保证了x数据的有效性，当x很小时，由于太小超过数值有效性，用log(x+1)，计算得到结果为0，换作log1p则计算得到一个很小却不为0的结果。同时，使用log1p能避免复值（复值指一个自变量对应多个因变量）的问题。
问题4：什么时候使用expm1？
在上面归一化的处理中，我们应用了log1p做数值变化，以使基因的原始count去偏态，符合正态分布。那么，也引出来一个问题，我们什么时候使用expm1？
场景1：我在《seurat-AverageExpression()源码解析》记录了一条，在AverageExpression()源码中，对归一化的基因表达矩阵seurat_obj[['RNA']]@data进行expm1处理，这个怎么理解呢？我在网上也看到网友对此问题的疑惑，为什么“take un-log data when calculate average expression”。
场景2：在seurat教程中，计算线粒体比例也用到了expm1
mito.genes <- grep("^MT-", rownames(pbmc@data), value = T)
percent.mito <- colSums(expm1(pbmc@data[mito.genes, ]))/colSums(expm1(pbmc@data))
#AddMetaData adds columns to object@data.info, and is a great place to stash QC stats
pbmc <- AddMetaData(pbmc, percent.mito, "percent.mito")
VlnPlot(pbmc, c("nGene", "nUMI", "percent.mito"), nCol = 3)
在seurat的代码中会出现对已经log1p归一化的数值进行expm1，比如线粒体的占比和AverageExpression的处理，背后的逻辑怎么理解？
在网上找到一些零星的信息：
The expm1 does un-log the data, but the normalization persists (this would be lost in the counts slot)
expm1() transformed in order to recover normalized values not in log scale
这块没有想的很明白，后面再补充
参考链接：
https://blog.csdn.net/weixin_48135624/article/details/114714449
https://www.researchgate.net/post/For-Seurat-in-the-log-normalize-step-of-sc-RNA-seq-data-what-does-the-scaling-value-imply
