Single-cell best practices

7 归一化

7.1 动机

到目前为止，我们从数据集中删除了低质量细胞、环境RNA污染和双联体，并且数据以“细胞 x 基因”数字矩阵的形式呈现，作为计数矩阵。 在scRNA-seq实验中转录本的捕获，逆转录和测序，其中每一步都会对相同细胞的测量深度增加变异性，这让数据集和计数矩阵间存在很大方差。

标准化的目的就是消除其他技术引起的变异，如测序深度、细胞裂解和逆转录效率的差异，使得我们得到的分析结果不受技术噪音的影响。

Poisson分布用于描述UMI计数数据的平均值，即每个基因在样本中的平均表达水平。Poisson分布是一种描述离散随机变量的分布，它的参数即为平均值。
Gamma分布用于描述UMI计数数据的过度分散，即UMI计数数据的变异性。Gamma分布是一种连续分布，通常用于描述正数数据的变异性。

Poisson分布用于描述UMI计数数据的基本分布特征，而Gamma分布描述了在Poisson分布之上的额外方差。这种组合模型更能够捕捉到UMI计数数据中的变异性和分布特性，从而提高了对数据的建模精度和分析效果。

“归一化”的预处理步骤旨在通过将“UMI的方差”缩放到指定范围，来调整数据集中的原始UMI计数以实现模型建模。

本篇将介绍三种不同的归一化技术：移位对数变换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方法用于计算序列中相邻元素之间的差异或变化。
移位对数变换是一种平滑数据的方法，通常用于减小数据的变化幅度，使其更加稳定。
使用delta方法来计算数据的变化量，然后对这些变化量应用移位对数变换方法进行平滑处理。

log1p()？ (伪计数:数据平滑处理) log1p(x)=log(1+x)
处理小于1的值，避免产生无限大或非数值的结果,保持相对比例稳定。

移位对数是一种快速归一化技术，优于其他揭示数据集潜在结构的方法（特别是在进行主成分分析时），并且有利于方差的稳定性，以进行后续的降维和差异表达基因的识别。

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

##进行总数标准化，然后对标准化后的数据进行对数变换
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 方法。

移位对数变换与scran的区别是：

移位对数变换通过取对数操作(通常使用log1p()函数来处理小于1的值)降低数据的方差来调整数据的分布形态，以提高数据的稳定性和可分析性。

Scran 使用反卷积方法估计每个细胞的大小因子（size factor），然后对原始的UMI计数数据进行归一化，从而消除细胞间的技术差异，旨在消除由于细胞大小差异引起的技术噪音和偏差，对单细胞数据进行批次效应和细胞大小效应的校正，从而提高数据的质量和可比性。

将细胞划分为池，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残差法是一种基于统计学原理的数据标准化方法，通过计算每个基因在每个细胞中的Pearson残差，利用“正则化负二项式回归”的Pearson残差来计算数据中的技术噪声模型，来消除细胞之间的技术噪音和批次效应同时保留了数据集中的细胞异质性。

##计算每个基因在每个细胞中的 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 特征选择

8.1 动机

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

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

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

根据偏差值进行特征选择。

##设置环境
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()