5. 归一化Normalization
5.1 前言
到目前为止，我们从数据集中删除了低质量细胞、环境RNA污染和双联体，并且数据以“细胞 x 基因”数字矩阵的形式呈现，作为计数矩阵。这些计数代表scRNA-seq实验中捕获、逆转录和测序的分子。每个步骤都会给相同细胞的测量计数深度增加一定程度的变异性，这意味着数据集以及计数矩阵仍然包含变化很大的方差项。分析数据集通常具有挑战性，因为许多统计方法假设数据具有统一的方差结构。
“归一化”的预处理步骤旨在通过将可观察方差缩放到指定范围来调整数据集中的原始计数以实现可变采样效果。实践中使用了多种复杂性各异的归一化技术。它们的设计方式大多是为了后续分析任务及其适用于基础统计方法。
本篇将介绍三种不同的归一化技术：移位对数变换shifted logarithm transformation、scran归一化和皮尔逊残差的解析逼近analytic approximation of Pearson residuals。移位对数变换有利于稳定方差，以利于后续降维和差异表达基因的识别。Scran经过广泛测试并用于批量校正任务，皮尔逊残差的解析逼近非常适合选择生物可变基因和鉴定稀有细胞类型。
我们首先导入所有必需的Python包并加载我们过滤低质量细胞、去除环境 RNA 和评分双联体的数据集。
import scanpy as sc
import numpy as np
import seaborn as sns
from matplotlib import pyplot as plt
import anndata2ri
import logging
from scipy.sparse import issparse
import rpy2.rinterface_lib.callbacks as rcb
import rpy2.robjects as ro
sc.settings.verbosity = 0
sc.settings.set_figure_params(
dpi=80,
facecolor="white",
# color_map="YlGnBu",
frameon=False,
)
rcb.logger.setLevel(logging.ERROR)
ro.pandas2ri.activate()
anndata2ri.activate()
%load_ext rpy2.ipython
adata = sc.read(
filename="s4d8_quality_control.h5ad",
backup_url="https://figshare.com/ndownloader/files/40014331",
)
我们现在可以检查我们在质量控制过程中已经计算出的原始计数的分布。在标准单细胞分析流程中可以忽略此步骤，但有助于理解不同的标准化概念。
p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False)
5.2 移位对数变换
我们介绍的第一个归一化技术是基于delta方法的移位对数。Delta方法将非线性函数 f(Y) 应用于原始计数 Y，旨在使数据集中的方差更加相似。
移位对数变换是一种快速归一化技术，优于其他揭示数据集潜在结构的方法（特别是在进行主成分分析时），并且有利于稳定方差，以进行后续的降维和差异表达基因的识别。我们现在将检查如何将此标准化方法应用于我们的数据集。通过运行pp.normalized_total设置为target_sum=None，可以使用scanpy方便地调用移位对数。我们将inplace参数设置为False，因为我们想在本教程中探索三种不同的标准化技术。第二步,使用缩放后的计数，获得了第一个归一化计数矩阵。
scales_counts = sc.pp.normalize_total(adata, target_sum=None, inplace=False)
# log1p transform
adata.layers["log1p_norm"] = sc.pp.log1p(scales_counts["X"], copy=True)
现在，我们可以检查应用移位对数变换后计数分布如何变化，并将其与原始（但已过滤）数据集中的总计数进行比较。
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
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.show()
第二种归一化方法也是基于delta方法，是scran归一化方法。它与移位对数变换唯一的区别是Scran 利用反卷积方法根据细胞池基因的线性回归来估计大小因子。这种方法旨在更好地解释数据集中所有细胞计数深度的差异。
细胞被分成多个池，Scran使用基因的线性回归来估计基于池的大小因子。Scran针对批量校正任务进行了广泛的测试，并且可以使用相应的R包轻松调用。
from scipy.sparse import csr_matrix, issparse
%%R
library(scran)
library(BiocParallel)
scran需要粗聚类输入来提高尺寸因子估计性能。在本教程中，我们使用简单的预处理方法并以低分辨率对数据进行聚类，以获得尺寸因子估计的输入。基本预处理包括假设所有大小因子相等（文库大小归一化为每百万计数 - CPM）以及对计数数据进行对数转换。
# Preliminary clustering for differentiated normalisation
adata_pp = adata.copy()
sc.pp.normalize_total(adata_pp)
sc.pp.log1p(adata_pp)
sc.pp.pca(adata_pp, n_comps=15)
sc.pp.neighbors(adata_pp)
sc.tl.leiden(adata_pp, key_added="groups")
现在，我们将data_mat和计算组添加到R环境中。
data_mat = adata_pp.X.T
# convert to CSC if possible. See https://github.com/MarioniLab/scran/issues/70
if issparse(data_mat):
if data_mat.nnz > 2**31 - 1:
data_mat = data_mat.tocoo()
else:
data_mat = data_mat.tocsc()
ro.globalenv["data_mat"] = data_mat
ro.globalenv["input_groups"] = adata_pp.obs["groups"]
我们现在还可以删除anndata对象的副本，因为我们获得了运行scran所需的所有对象。
del adata_pp
现在，我们根据之前计算的细胞组来计算尺寸因子。
%%R -o size_factors
size_factors = sizeFactors(
computeSumFactors(
SingleCellExperiment(
list(counts=data_mat)),
clusters = input_groups,
min.mean = 0.1,
BPPARAM = MulticoreParam()
)
)
我们将size_factors保存在.obs中，现在能够标准化数据并随后应用log1p转换。
adata.obs["size_factors"] = size_factors
scran = adata.X / adata.obs["size_factors"].values[:, None]
adata.layers["scran_normalization"] = csr_matrix(sc.pp.log1p(scran))
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
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["scran_normalization"].sum(1), bins=100, kde=False, ax=axes[1]
)
axes[1].set_title("log1p with Scran estimated size factors")
plt.show()
5.3 皮尔逊残差的解析逼近
我们在本篇中介绍的第三种归一化技术是皮尔逊残差的解析逼近。这种标准化技术的原因是观察到scRNA-seq数据中的细胞间变异可能会因生物异质性与技术效应而混淆。该方法利用“正则化负二项式回归”中的皮尔逊残差来计算数据中的技术噪声模型。它明确地将计数深度添加为广义线性模型中的协变量。在对不同标准化技术的独立比较中表明，该方法消除了采样效应的影响，同时保留了数据集中的细胞异质性。解析Pearon残差在scanpy中实现，可以直接在原始计数矩阵上计算。
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=(10, 5))
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.show()
我们对数据集应用了不同的标准化技术，并将它们作为单独的层保存到我们的anndata对象中。根据下游分析任务，可以使用不同的归一化层并评估结果。
adata.write("s4d8_normalization.h5ad")
