特征选择

10. 特征选择#

10.1. 动机#

我们现在得到了一份归一化后的数据表示，它在保留生物学异质性的同时，降低了基因表达中的技术性采样效应。scRNA-seq 数据集通常包含多达 30,000 个基因（取决于物种），而到目前为止，我们只去除了在少于 20 个细胞中检测到的基因。然而，剩余基因中有许多并不具信息量，且大多为零计数。因此，标准的预处理流程会包含特征选择这一步，旨在剔除那些不具信息量、可能无法代表样本间有意义生物学变异的基因。

Feature selection — 图10.1 特征选择通常指只选取一部分相关特征（基因）的过程，这些特征可以是最具信息量的、变异最大的，或偏差最高的。#

通常，scRNA -seq 实验聚焦于某一特定组织，因此只有一小部分基因具有信息量，且在生物学上存在变异。传统方法和流程要么计算所有基因的变异系数（高变基因），要么计算所有基因的平均表达水平（高表达基因），从中选出 500–2000 个基因，并将这些特征用于下游分析步骤。然而，这些方法对此前所用的归一化技术非常敏感。如前所述，以往的预处理流程包括用 CPM 归一化、随后做对数变换。但由于对数变换无法作用于精确的零值，分析人员往往会先给所有归一化计数都加上一个小的 伪计数（例如 1，即 log1p），再对数据取对数。然而，伪计数的选择是任意的，可能给变换后的数据引入偏倚。这种任意性进而也会影响特征选择，因为观测到的变异性取决于所选的伪计数。一个接近零的小伪计数值，会增大零计数基因的方差 [Townes et al., 2019]。

Germain 等人转而提出使用偏差（deviance）来进行特征选择，它作用于原始计数 [Germain et al., 2020]。偏差可以用闭式（closed form）计算，用来衡量某个基因在各细胞间是否呈现恒定的表达谱——这类基因不具信息量。表达恒定的基因可由一个多项式零模型描述，并用二项偏差来近似。在各细胞间信息量很高的基因会有较高的偏差值，表明零模型对其拟合不佳（即它们在各细胞间并非恒定表达）。该方法据此按偏差值对所有基因排序，只选取高偏差的基因。

如前所述，偏差可以用闭式计算，并由 R 包 scry 提供。

我们从建立环境开始。

import logging

import lamindb as ln
import matplotlib.pyplot as plt
import numpy as np
import rpy2.rinterface_lib.callbacks as rcb
import rpy2.robjects as ro
import rpy2.robjects.packages as rpackages
import scanpy as sc
import seaborn as sns
from rpy2.robjects import default_converter, numpy2ri, pandas2ri, r
from rpy2.robjects.conversion import localconverter

# Suppress verbose logging from Scanpy
sc.settings.verbosity = 0

# Set figure parameters for clean, minimal plots
sc.settings.set_figure_params(dpi=80, facecolor="white", frameon=False)

assert ln.setup.settings.instance.slug == "theislab/sc-best-practices"

ln.track()

rcb.logger.setLevel(logging.ERROR)


%load_ext rpy2.ipython

%%R
library(scry)
library(SingleCellExperiment)

接着，我们加载在前一个章节中已经归一化好的数据集。偏差作用于原始计数，因此无需把 adata.X 替换为某个归一化层，我们可以直接使用归一化笔记本中所存储的对象。

af = ln.Artifact.connect("theislab/sc-best-practices").get(
    key="preprocessing_visualization/s4d8_normalization.h5ad", is_latest=True
)
adata = af.load()
adata

AnnData object with n_obs × n_vars = 14814 × 20223
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'outlier', 'mt_outlier', 'soupx_groups', 'scDblFinder_score', 'scDblFinder_class', 'size_factors'
    var: 'gene_ids', 'feature_types', 'genome', 'interval', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
    layers: 'analytic_pearson_residuals', 'counts', 'log1p_norm', 'scran_normalization', 'soupX_counts'

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

X_sparse = adata.X.T.tocoo()

Matrix = rpackages.importr("Matrix")

with localconverter(ro.default_converter + pandas2ri.converter + numpy2ri.converter):
    ro.globalenv["obs"] = adata.obs
    ro.globalenv["var"] = adata.var

i, j = X_sparse.row, X_sparse.col
x = X_sparse.data

ro.globalenv["i"] = ro.IntVector((i + 1).tolist())  # R is 1-indexed
ro.globalenv["j"] = ro.IntVector((j + 1).tolist())
ro.globalenv["x"] = ro.FloatVector(x.tolist())

r("X <- sparseMatrix(i = i, j = j, x = x, dims = c({}, {}))".format(*X_sparse.shape))

<rpy2.robjects.methods.RS4 object at 0x164985550> [25]
R classes: ('dgCMatrix',)

现在我们可以直接在未归一化的计数矩阵上调用基于偏差的特征选择，并将二项偏差值导出为一个向量。

%%R
sce <- SingleCellExperiment(
  assays = list(X = X),
  colData = obs,
  rowData = var
)

sce <- devianceFeatureSelection(sce, assay = "X")

with localconverter(default_converter + pandas2ri.converter + numpy2ri.converter):
    binomial_deviance = ro.r("rowData(sce)$binomial_deviance")

接下来，我们对该向量进行排序，选出偏差最高的前 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()

../_images/12a462e82a4bbd79a6adff85dde93b6ac9c3e48dd769265800bd4ca0f573947f.png

我们观察到，平均表达较高的基因被选为高偏差基因。这与以下研究的经验观察一致：[Townes et al., 2019]。

af = ln.Artifact.from_anndata(
    adata,
    key="preprocessing_visualization/s4d8_feature_selection.h5ad",
    description="anndata after feature selection",
).save()
af

10.2. 参考文献#

[fsGSR20]

Pierre-Luc Germain, Anthony Sonrel, and Mark D. Robinson. pipeComp, a general framework for the evaluation of computational pipelines, reveals performant single cell rna-seq preprocessing tools. Genome Biology, 21(1):227, September 2020. URL: https://doi.org/10.1186/s13059-020-02136-7, doi:10.1186/s13059-020-02136-7.

[fsTHAI19] (1,2)

F. William Townes, Stephanie C. Hicks, Martin J. Aryee, and Rafael A. Irizarry. Feature selection and dimension reduction for single-cell rna-seq based on a multinomial model. Genome Biology, 20(1):295, Dec 2019. URL: https://doi.org/10.1186/s13059-019-1861-6, doi:10.1186/s13059-019-1861-6.

10.3. 贡献者#

我们衷心感谢以下人员的贡献：

10.3.1. 作者#

Anna Schaar
Seo H. Kim

10.3.2. 审阅者#

Lukas Heumos