2.2 Normalization¶
Motivation¶
到目前为止，我们从数据集中移除了低质量细胞、环境RNA污染和双重细胞，并将数据存储为形状为细胞 x 基因的数值矩阵的计数矩阵。这些计数表示在scRNA-seq实验中对分子进行捕获、反转录和测序的过程。每个步骤都会为相同细胞的测量计数深度增加一定的变异性，因此计数数据中细胞间的基因表达差异可能只是由于采样效应造成的。这意味着数据集以及计数矩阵仍然包含广泛变化的方差项。分析数据集通常是具有挑战性的，因为许多统计方法假设数据具有均匀的方差结构。
{admonition}
对于UMI数据，理论上和实证上建立的模型是Gamma-Poisson分布，这意味着均值-方差关系是二次的，$Var[Y] = \mu + \alpha \mu^2$，其中均值为$\mu$，过度离散参数为$\alpha$。当$\alpha=0$时，这是泊松分布，$\alpha$描述了在泊松分布之上的附加方差。
“标准化”的预处理步骤旨在通过将可观测方差调整到指定范围来校正数据集中由于采样效应导致的原始计数变动。实践中使用了几种不同的标准化技术，它们在复杂性上有所不同。它们的设计主要是为了使后续分析任务及其基础统计方法适用。
Ahlmann-Eltze和Huber最近发表的一项基准研究【2023】比较了22种不同的单细胞数据转换方法。该基准研究根据与真实情况的细胞图重叠情况比较了不同标准化技术的性能。我们要强调的是，尚缺乏一个完整的基准来比较标准化对各种不同后续分析任务的影响。我们建议分析师根据后续分析任务谨慎选择标准化方法。
本章将介绍三种不同的标准化技术，即移位对数转换、scran标准化和Pearson残差的解析近似。移位对数转换有助于稳定后续降维和差异基因表达鉴定的方差。scran广泛用于批次校正任务，而解析的Pearson残差非常适合选择生物学上可变的基因和识别稀有细胞类型。
首先，我们将导入所有必需的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)
移位对数变换¶
我们将介绍的第一种标准化技术是基于delta方法的移位对数变换【1938】。delta方法对原始计数$Y$应用非线性函数$f(Y)$，旨在使数据集中的方差更加相似。
移位对数变换通过以下方式解决这个问题：
$$f(y) = \log(\frac{y}{s}+y_0)$$
其中$y$是原始计数，$s$是所谓的大小因子，$y_0$描述了一个伪计数。大小因子是为每个细胞确定的，以解释采样效应和不同细胞大小的变异。细胞$c$的大小因子可以通过以下公式计算：
$$s_c = \frac{\sum_g y_{gc}}{L}$$
其中$g$索引不同的基因，$L$描述目标总和。有不同的方法可以从数据中确定大小因子。在本节中，我们将利用scanpy的默认设置，其中$L$是数据集中原始计数深度的中位数。许多分析模板使用固定值$L$，例如$L=10^5$或$L=10^6$，这些值通常被称为每百万计数（CPM）。对于初学者来说，这些值可能看起来是任意的，但它可能导致比单细胞数据集中通常看到的更大的过度离散。
{admonition}
过度离散描述了数据集中存在比预期更大的变异性。
移位对数变换是一种快速的标准化技术，在揭示数据集的潜在结构（特别是在随后进行主成分分析时）方面表现优于其他方法，并且有助于稳定后续降维和差异表达基因鉴定的方差。我们现在将检查如何将这种标准化方法应用于我们的数据集。通过运行pp.normalized_total并将target_sum设置为None，可以方便地调用移位对数变换。我们将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通过计算 $f(y) = \log(\frac{y}{s}+y_0)$ 来遵循与移位对数变换相同的原理，其中$y$是原始计数，$s$是大小因子，$y_0$描述了一个伪计数。唯一的不同是，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环境中。以下是将Python中的数据导入到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()
分析Pearson残差¶
我们在本章介绍的第三种标准化技术是Pearson残差的解析近似。该标准化技术的动机源于观察到scRNA-seq数据中的细胞间变异可能受到生物异质性与技术效应的混淆。该方法利用“正则化负二项回归”的Pearson残差来计算数据中技术噪声的模型。它在广义线性模型中明确地将计数深度作为协变量添加进去。{cite}norm:germain_pipecomp_2020在独立比较不同标准化技术的研究中表明，该方法在去除采样效应影响的同时保留了数据集中的细胞异质性。值得注意的是，分析Pearson残差不需要下游的启发式步骤，如伪计数添加或对数变换。
该方法的输出是可以为正或负的标准化值。对于某个细胞和基因，负残差表明观察到的计数比基因的平均表达和细胞测序深度所期望的要少。正残差则表明观察到的计数比期望的要多。分析Pearson残差在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")