6. 特征选择Feature selection
6.1 前言
我们现在有了标准化的数据表示，它仍然保留了生物异质性，但减少了基因表达中的技术采样效应。单细胞RNA-seq数据集通常包含多达30000个基因，到目前为止，我们仅删除了至少20个细胞中未检测到的基因。然而，许多剩余的基因没有提供信息，并且大多包含零计数。因此，标准预处理流程涉及特征选择步骤，旨在排除可能不代表样本间有意义的生物变异的无信息基因。
通常，scRNA-seq实验和生成的数据集集中于一个特定组织，因此，只有一小部分基因具有信息性和生物学可变性。传统的方法和流程要么计算所有基因的变异系数（高度可变的基因）或平均表达水平（高度表达的基因），以获得500-2000个选定的基因，并将这些特征用于下游分析步骤。然而，这些方法对之前使用的标准化技术高度敏感。如前所述，以前的预处理工作流程包括使用CPM进行规范化以及随后的log转换。但由于log转换不可能精确为零，因此分析人员通常会在对数据进行对数转换之前向所有归一化计数添加一个小的伪计数，例如 1 (log1p)。然而，伪计数的选择是任意的，并且可能会给转换后的数据带来偏差。这种任意性也会影响特征选择，因为观察到的变化取决于所选的伪计数。接近于零的小伪计数值会增加计数为零的基因的方差。
相反，也有研究建议使用偏差deviance进行特征选择，该特征选择适用于原始计数。偏差可以以封闭形式计算，并量化基因是否在细胞中表现出恒定的表达谱。具有恒定表达的基因由多项零模型描述，它们通过二项式偏差进行近似。细胞间信息丰富的基因将具有较高的偏差值，这表明零模型拟合不佳（即它们不显示细胞间的恒定表达）。根据偏差值，该方法然后对所有基因进行排序并仅获得高度偏差的基因。
如前所述，偏差可以以封闭形式计算，并在 R 包 scry 中提供。
我们首先设置我们的环境。
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)
接下来，我们加载已经标准化的数据集。偏差deviance适用于原始计数，因此无需将adata.X替换为标准化层之一，但我们可以直接使用该对象。
adata = sc.read(
filename="s4d8_normalization.h5ad",
backup_url="https://figshare.com/ndownloader/files/40015741",
)
和前面类似，将AnnData对象存储在R环境中。
ro.globalenv["adata"] = adata
我们现在可以直接在非归一化计数矩阵上调用带有偏差的特征选择，并将生物偏差值导出为向量。
%%R
sce = devianceFeatureSelection(adata, assay="X")
binomial_deviance = ro.r("rowData(sce)$binomial_deviance").T
下一步，我们现在对向量进行排序，选择前4000个高度异常基因，并将它们保存为.var中的附加列“highly_deviant”。我们还保存计算出的二项式偏差，以防我们之后想要子选择不同数量的高度可变的基因。
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")
我们通过绘制基因的分散度与均值的关系图来检查我们的结果，并通过“highly_deviant”绘制颜色。
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()
我们观察到，具有高平均表达的基因被选择为高度异常的。
adata.write("s4d8_feature_selection.h5ad")
