Single-cell best practices

7 Normalization

7.1 动机

到目前为止，我们从数据集中删除了低质量细胞、环境RNA污染和双联体，并且数据以“细胞 x 基因”数字矩阵的形式呈现，作为计数矩阵。

计数矩阵中的每个计数代表成功捕获、逆转录和测序一个细胞 mRNA 分子。

技术变异：在 scRNA-seq 实验中，数据会受到多种技术变异的影响。例如，不同细胞中的 RNA 捕获效率、逆转录效率和测序深度可能不同，这些技术变异会导致原始计数数据的偏差。

归一化的目的：为了比较不同细胞之间的基因表达水平，必须对这些技术变异进行校正。归一化步骤就是为了消除上述技术误差，使不同细胞的总 RNA 计数达到相同的水平，从而使细胞间的基因表达水平具有可比性。

处理这种异方差的一种方法是对采样分布进行显式建模：负二项分布（Gamma-Poisson 分布），其中参数推断实现繁琐，计算成本高。没有充分利用数据中经常出现零和其他小数这一事实。这些局限性阻碍了该模型的采用。（Ahlmann-Eltze, C. & Huber, W. glmGamPoi: fitting gamma-Poisson generalized linear models on single cell count data. Bioinformatics https://doi.org/10.1093/bioinformatics/btaa1009 (2020).）

Var[Y] ，指的是某一基因在不同细胞中的表达水平的变异程度。

Var[Y]=μ，方差等于均值，这意味着，对于单个基因而言，在不同细胞中的表达水平变异性完全由其平均表达水平决定，没有额外的变异来源。

实际应用中的偏差：单细胞RNA-seq数据分析中，往往观察到数据的方差比均值大得多，种现象称为过度离散，这时，Poisson分布模型就不再适用了，需要采用更复杂的模型，如负二项分布（Gamma-Poisson 分布），来捕捉这种额外的变异性。

Var[Y]=μ+αμ 2 反映了负二项分布中方差随均值变化的情况，μ 项表示 Poisson 分布中的基础变异性，αμ 2 表示额外的变异性，反映了数据中存在的过度离散。

另外一种方法：使用方差稳定变换作为预处理步骤，然后使用许多现有的统计方法。

本篇将介绍三种不同的归一化技术：移位对数变换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)

counts

7.2.1 移位对数变换shifted logarithm transformation

移位对数变换 基于delta方法。

delta方法：
1、非线性函数的选择，如对数函数log(1+x)，来变换原始计数数据。
2、近似方差计算，Var(g(X))≈(g ′ (μ)) 2 Var(X)， Var(g(X)) 是应用非线性函数变换后的方差，μ 是 X 的均值，g ′ (μ) 是在均值 μ 处的导数，Var(X) 是原始计数数据的方差。

Delta方法应用了非线性函数g(X) ，旨在使原始计数X中的差异更加相似。

log1p()？ (伪计数:数据平滑处理) log1p(x)=log(1+x)
1、处理零值，避免了无法计算的问题。
2、缩小动态范围，减少了高表达基因的极端值的影响，使得数据更加平滑。
3、均衡平滑数据，使数据接近正态分布。

(Ahlmann-Eltze C, Huber W. Comparison of transformations for single-cell RNA-seq data. Nat Methods. 2023;20(5):665-672. doi:10.1038/s41592-023-01814-1)

伪计数y 0 = 1 / (4 α ) ，过度离散α可以用于确定合适的L值

每个细胞的大小因子表示该细胞中相对偏差的估计，因此，将其计数除以其大小因子应消除该偏差。

大小因子的计算：分子是细胞C的UMI总数，g索引基因，L是所有细胞umi总数的平均值。 有时，使用固定值代替L。例如，Seurat使用L=10,000，尽管L的选择看似任意，但它非常重要。

将原始计数除以大小因子（size factor）是一种标准化步骤，目的是为了校正不同细胞之间测序深度的差异，使得基因表达水平在不同细胞之间具有可比性。

移位对数是一种快速归一化技术，优于其他揭示数据集潜在结构的方法（特别是在进行主成分分析时），并且有利于方差的稳定性，以进行后续的降维和差异表达基因的识别。
本节分析中使用了scanpy 默认值作为数据集中的原始计数深度的中位数，即L=10的五次方(由于数据的稀疏性,总计数的中位数作为L可能会受到数据中极端值的影响)

将移位对数变换方法应用于我们的数据集。

##进行总数标准化，将每个细胞的总数标准化到一个相对的尺度上，然后对标准化后的数据进行对数变换
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()

归一化前后对比

7.2.2 scran归一化

高频率的空胞事件会干扰标准化使大小因子无意义，我们可以对每一个细胞去除有问题的基因，但是不同细胞中基因计数为零的数量不同，可能会引入偏差。

scran归一化法也是基于 delta 方法。对多个细胞的合并计数进行归一化，其中通过对细胞求和来减少有问题的零的发生率。然后对合并的大小因子进行去卷积以推断单个细胞的大小因子。

(Lun AT, Bach K, Marioni JC. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Genome Biol. 2016;17:75. Published 2016 Apr 27. doi:10.1186/s13059-016-0947-7)

Scran 方法是基于反卷积的方法：先合并细胞避免drop-out问题，使用基因的线性回归来估计基于池的大小因子，然后将基于池的大小因子解卷积到基于细胞的大小因子，以进行特定细胞的标准化

移位对数变换与scran的区别是：移位对数变换主要针对每个基因的表达量进行调整，而Scran归一化则是针对每个细胞的总计数进行校正。

将细胞划分为池，Scran使用基因的线性回归来估计基于池的大小因子.

from scipy.sparse import csr_matrix, issparse 
%%R
library(scran)
library(BiocParallel) 

##对数据进行总数标准化,对数变换,主成分分析,计算邻居, Leiden 聚类
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")

##转秩后，检查转换后的数据是否为稀疏矩阵，将转换后的数据矩阵传递到 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"]

del adata_pp

##计算单细胞数据中每个细胞的大小因子，并将结果传递回 Python 中的 size_factors 变量中
%%R -o size_factors
size_factors = sizeFactors(
    computeSumFactors(
        SingleCellExperiment(
            list(counts=data_mat)), 
            clusters = input_groups,
            min.mean = 0.1,
            BPPARAM = MulticoreParam()
    )
)

##进行 scran 方法的标准化处理：每个细胞的RNA计数除以其对应的大小因子，后对数变换转为稀疏矩阵。
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))

##绘制两个直方图，比较了原始数据的总计数与经过 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()

7.3 Pearson 残差

正则化负二项式回归(伽马泊松)”的Pearson残差，其中细胞测序深度在广义线性模型中用作协变量，成功地消除了下游分析中技术特征的影响，同时保留了生物学异质性。

(Ahlmann-Eltze C, Huber W. Comparison of transformations for single-cell RNA-seq data. Nat Methods. 2023;20(5):665-672. doi:10.1038/s41592-023-01814-1)

(Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 2019;20(1):296. Published 2019 Dec 23. doi:10.1186/s13059-019-1874-1)

rgc 是某基因在细胞c中的 Pearson 残差
ygc 是某基因在细胞c中的观测值umi数
μgc 是某基因在细胞c中的期望值umi数，βg,intercept 和 βg,slope 是基因 g 的截距和斜率参数
分母是是参数为μgc和αg的 gamma-Poisson 随机变量的标准差

将计数深度添加为广义线性模型中的协变量

我们使用UMI计数的回归模型来校正细胞之间的测序深度差异并标准化数据,单独对每个基因进行建模会导致过度拟合，特别是对于仅在一小部分细胞中检测到的低丰度基因，并且以高方差进行建模。我们认为这是对真实方差的高估，因为这是由样本中的细胞类型异质性驱动的，而不是由于自变量log10m的细胞间变异性。为了避免这种过拟合，我们通过在基因之间共享信息来正则化所有模型参数，包括NB离散参数θ。 分母是是参数为μgc和αg的 gamma-Poisson 随机变量的标准差

##计算每个基因在每个细胞中的 Pearson 残差，并对这些残差进行标准化处理,转换为稀疏矩阵保存
analytic_pearson = sc.experimental.pp.normalize_pearson_residuals(adata, inplace=False)
adata.layers["analytic_pearson_residuals"] = csr_matrix(analytic_pearson["X"])

##绘制两个直方图，比较原始数据的总计数与Pearson 残差方法处理后的数据的分布情况。
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()
adata.write("s4d8_normalization.h5ad")

Analytic Pearson residuals 此方法的输出是归一化值，可以是正值，也可以是负值。它们代表了每个基因在每个细胞中的 Pearson 残差，即原始数据与期望值之间的偏差(这里的期望值是根据整个数据集中的基因表达量来计算的。)。正值表示相应基因在相应细胞中的表达量高于期望值，而负值表示相应基因在相应细胞中的表达量低于期望值。

8 Feature selection

8.1 动机

我们现在有了标准化的数据表示，它仍然保留了生物异质性，但减少了基因表达中的技术采样效应。单细胞 RNA-seq 数据集通常包含多达 30,000 个基因，到目前为止，我们仅删除了至少 20 个细胞中未检测到的基因。然而，许多剩余的基因没有提供信息，并且大多包含零计数。因此，标准预处理流程涉及特征选择步骤，旨在排除可能不代表样本间有意义的生物变异的无信息基因。

特征选择通常描述了仅选择相关特征子集的过程，这些特征可以是信息量最大、变化最大或偏差最大的特征。

传统方法和管道要么计算变异系数（高变基因）要么计算所有基因的平均表达水平（高表达基因），以获得 500-2000 个选定的基因，并将这些特征用于其下游分析步骤。但是，这些方法对之前使用的归一化技术高度敏感。

log1p：接近于零的小伪计数值会增加计数为零的基因的方差。由于计数为零的基因通常是在大多数细胞中表达较低的基因，因此对这些基因添加小的伪计数值可能会导致其在对数转换后的数据中得到过度放大，进而增加其方差。

Germain 等人建议使用偏差进行特征选择，这适用于原始计数。
根据偏差值进行特征选择。

首先使用多项式零模型来描述具有恒定表达的基因，并使用二项式偏差来近似这些基因。然后，根据基因的偏差值对所有基因进行排名，仅选择具有高度偏差的基因作为最终的特征集合。
这种基于偏差的特征选择方法可以帮助识别在样本中具有显著变化的基因，从而过滤掉那些在不同样本间表达变化较小的基因。

##设置环境
import scanpy as sc
import anndata2ri
import logging
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

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",
    frameon=False,
)

rcb.logger.setLevel(logging.ERROR)
ro.pandas2ri.activate()
anndata2ri.activate()

%load_ext rpy2.ipython

%%R
library(scry)


##读入数据
adata = sc.read(
    filename="s4d8_normalization.h5ad",
    backup_url="https://figshare.com/ndownloader/files/40015741",
)

##将 AnnData 对象保存在 R 环境中
ro.globalenv["adata"] = adata 


##在R中进行特征选择
%%R
sce = devianceFeatureSelection(adata, assay="X")

##提取sce的binomial_deviance，转置，保存高度偏差基因和偏差值。
binomial_deviance = ro.r("rowData(sce)$binomial_deviance").T

##排序，选择前4000个高变基因，保存。
idx = binomial_deviance.argsort()[-4000:]
mask = np.zeros(adata.var_names.shape, dtype=bool)
mask[idx] = True

adata.var["highly_deviant"] = mask
adata.var["binomial_deviance"] = binomial_deviance

##可视化结果，使用 scanpy 函数来计算所有细胞中每个基因的平均值和离散度。  
sc.pp.highly_variable_genes(adata, layer="scran_normalization")

ax = sns.scatterplot(
    data=adata.var, x="means", y="dispersions", hue="highly_deviant", s=5
)
ax.set_xlim(None, 1.5)
ax.set_ylim(None, 3)
plt.show()

我们观察到，具有高平均表达的基因被选为高度偏差的基因。