参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

归一化的引入

在质量控制中，已经从数据集删除了低质量细胞。然而由于测序技术的限制，我们在样本中获得RNA时，经过分子捕获，逆转录和测序，这些步骤会影响同一种细胞的细胞间测序深度的变异性，因此，数据中的细胞间差异包含了这部分误差，等价于counts矩阵包含了变化很大的方差项。

归一化旨在通过将UMI counts的方差缩放到指定范围，以调整原始矩阵的UMI counts。目前有两种归一化方法：

• 移位对数：在大部分数据中表现良好，有利于稳定方差，进而利于降维和差异基因识别；
• 皮尔森残差的近似解析：保留生物学差异，有利于鉴定稀有细胞类型。

首先，我们加载数据：

import omicverse as ov
import scanpy as sc
import matplotlib.pyplot as plt

ov.utils.ov_plot_set()

adata = sc.read("./data/s4d8_quality_control.h5ad")
print(adata)

然后，可视化total_counts，这是描述一个细胞中发现的分子数量（UMI），通常也可以被认为是这个细胞的文库大小：

import seaborn as sns
plt.figure(figsize=(8, 6))
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False)
plt.show()

这可视化了原始计数UMI的分布，可以用于和之后归一化的分布对比。

移位对数

这里介绍基于delta方法的移位对数，delta方法应用$f(Y)$，使得原始计数$Y$中的差异被缩小：$$f(y)=log(\frac{y}{s}+y_0)$$其中，$s$是缩放因子，$y_0$是伪计数。每个细胞都有对应的缩放因子，细胞$c$的缩放因子记为：$$s_c=\frac{{\sum }_g{y}_{gc}}{L}$$其中，$g$代表不同的基因，$L$代表基因的计数总和。

使用移位对数归一化：

scales_counts = sc.pp.normalize_total(adata, target_sum=None, inplace=False)
print(scales_counts)
# log1p transform
adata.layers["log1p_norm"] = sc.pp.log1p(scales_counts["X"], copy=True)

可视化对比归一化前后：

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False, ax=axes[0])
axes[0].set_title("Total counts")
p2 = sns.histplot(adata.layers["log1p_norm"].sum(1), bins=100, kde=False, ax=axes[1])
axes[1].set_title("Shifted logarithm")
plt.savefig("./result/2-3.png")

我们发现UMI的最大值在1000左右，经过移位对数化后，UMI的分布近似正态分布。

皮尔森近似残差

scRNA-seq包含生物异质性和批次效应，移位对数更倾向于消除批次差距，皮尔森近似残差可以保留移位对数去除的信息。实验中发现，皮尔森近似残差计算非常慢。对于14814×20171的adata，移位对数花费5秒，皮尔森近似残差花费3分33秒。

归一化与可视化为：

from scipy.sparse import csr_matrix
analytic_pearson = sc.experimental.pp.normalize_pearson_residuals(adata, inplace=False)
adata.layers["analytic_pearson_residuals"] = csr_matrix(analytic_pearson["X"])

fig, axes = plt.subplots(1, 2, figsize=(8, 4))
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False, ax=axes[0])
axes[0].set_title("Total counts")
p2 = sns.histplot(adata.layers["analytic_pearson_residuals"].sum(1), bins=100, kde=False, ax=axes[1])
axes[1].set_title("Analytic Pearson residuals")
plt.savefig("./result/2-4.png")

注意，如果我们设置inplace=True时，我们归一化的计数矩阵会取代原anndata文件中的计数矩阵，即更改adata.X的内容。

相比移位对数，皮尔森近似残差归一化后的数据分布与原始数据更相似，所以保留了更多信息。

两个归一化方法的总结

移位对数和皮尔逊近似残差是两种用于归一化数据的方法，它们各自具有不同的特点：

1. 移位对数（Log-transformation）：

   • 特点：将原始数据的计数值进行对数转换，通常是加上一个小的常数（如1），以避免计数值为零时出现无穷大的情况。
   • 优点：可以有效地减小数据的偏斜，使其更符合正态分布假设。对于计数数据，对数转换也可以减小计数之间的差异，有助于更好地展现数据的模式和关系。
   • 缺点：对于一些数据分布，特别是存在大量低计数值的情况下，对数转换可能会引入噪音，使数据更难解释。此外，对数转换可能会导致丢失原始数据的一些信息。

2. 皮尔逊近似残差（Analytic Pearson residuals）：

   • 特点：利用正则化负二项回归得到的皮尔逊残差，通过计算数据中的技术噪声模型来归一化数据。
   • 优点：能够更准确地处理数据中的技术效应和生物异质性，避免了一些常见归一化方法可能引入的偏差。不需要额外的启发式步骤（如伪计数添加或对数转换）。
   • 缺点：相对于简单的对数转换方法，计算复杂度较高。

总的来说，移位对数适用于简单的数据集，对数转换可使数据更易于处理和分析；而皮尔逊近似残差则更适用于复杂的数据集，尤其是对于单细胞RNA测序数据很需要生物异质性的情况。