2.3 Feature selection¶

Motivation¶

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

特征选择通常描述只选择一部分相关特征的过程，这些特征可以是最具信息性的、最具变异性的或最偏离的。

通常，scRNA-seq实验和结果数据集集中在一个特定组织上，因此只有一小部分基因具有信息性和生物变异性。传统方法和流程要么计算所有基因的变异系数（高变异基因），要么计算平均表达水平（高表达基因），以获得500-2000个选择基因，并使用这些特征进行下游分析步骤。然而，这些方法对之前使用的标准化技术非常敏感。如前所述，以前的预处理工作流包括使用CPM进行标准化和随后的对数变换。但由于对数变换对精确的零值不可行，分析人员通常在对数据进行对数变换之前向所有标准化计数添加一个小的“伪计数”，例如1（log1p）。然而，选择伪计数是任意的，可能会引入偏差到变换后的数据中。这种任意性也会影响特征选择，因为观察到的变异性取决于选择的伪计数。接近零的小伪计数值会增加零计数基因的方差【2019】。

Germain等人提议使用偏差进行特征选择，这适用于原始计数【2020】。偏差可以通过闭式形式计算，并量化基因是否在细胞间显示出恒定的表达谱，因为这些基因不具信息性。具有恒定表达的基因由多项式零模型描述，它们通过二项式偏差进行近似。跨细胞高度信息性的基因将具有高偏差值，这表明零模型的拟合度较差（即，它们在细胞间不显示恒定表达）。根据偏差值，该方法对所有基因进行排名，仅获取高度偏离的基因。

如前所述，偏差可以通过闭式形式计算，并在R包scry中提供。

我们首先设置我们的环境。

In [1]:

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

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

In [2]:

%%R
library(scry)

%%R library(scry)

接下来，我们加载已标准化的数据集。由于偏差计算适用于原始计数，因此不需要用标准化层之一替换adata.X，可以直接使用在标准化笔记本中存储的对象。

In [3]:

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

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

与之前类似，我们将AnnData对象保存到R环境中。

In [4]:

ro.globalenv["adata"] = adata

ro.globalenv["adata"] = adata

我们现在可以直接在未标准化的计数矩阵上使用偏差进行特征选择，并将二项式偏差值导出为一个向量。

In [5]:

%%R
sce = devianceFeatureSelection(adata, assay="X")

%%R sce = devianceFeatureSelection(adata, assay="X")

In [6]:

binomial_deviance = ro.r("rowData(sce)$binomial_deviance").T

binomial_deviance = ro.r("rowData(sce)$binomial_deviance").T

下一步，我们对向量进行排序，选择前4000个高度偏离的基因，并将它们保存为.var中的一个额外列，命名为'highly_deviant'。此外，我们保存计算出的二项式偏差，以便以后可以选择不同数量的高度变异基因。

In [7]:

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

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函数计算每个基因在所有细胞中的均值和离散度。

In [8]:

sc.pp.highly_variable_genes(adata, layer="scran_normalization")

sc.pp.highly_variable_genes(adata, layer="scran_normalization")

我们通过绘制基因的离散度与均值的关系图，并用'高度偏离'进行着色来检查结果。

In [9]:

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()

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()

我们观察到，具有高平均表达的基因被选择为高度偏离的基因。这与【2019】的经验观察一致。

In [10]:

adata.write("s4d8_feature_selection.h5ad")

adata.write("s4d8_feature_selection.h5ad")