📚 课程目录第4章特征选择基础

特征选择基础：方差过滤、基于均值的过滤、稀疏性处理

特征选择，说白了就是替数据“瘦身”。高维转录组数据里，常常有大量基因表达量要么恒高、要么恒低——这些特征对模型毫无贡献。我入行那会儿接手过一个单细胞项目，13000个基因扔进聚类模型，跑了整宿没跑完，结果一看，超过一半是“噪声基因”。从那天起，特征选择就成了我每一套 pipeline 里雷打不动的第一步。

为什么要给 AnnData 做特征过滤？

AnnData 对象里存着基因×细胞矩阵，维度动辄数万×数千。直接拿来做特征学习，一是计算开销大，二是过拟合风险高。你可以这么理解：如果某个基因在所有细胞里表达量几乎一样，那它提供的信息量接近于零。这样的特征留着，只会稀释真正的信号。

我个人的经验是，先做粗筛，再精调。粗筛阶段，重点关注三类指标：方差、均值、稀疏度。下面一个个说。

方差过滤：剔除“平缓”基因

方差衡量的是基因表达值在细胞间的波动幅度。方差太小，说明这个基因要么一直高表达，要么一直低表达——对区分细胞亚群没用。

在 AnnData 里，我习惯用 sc.pp.filter_genes_dispersion 来按离散度筛选。不过方差过滤更直接：直接计算每个基因的方差，然后保留方差最大的前 N 个。

import scanpy as sc
import numpy as np

# 假设 adata 已经过基本预处理
sc.pp.highly_variable_genes(adata, n_top_genes=2000, subset=False)
# 这会计算离散度（方差/均值归一化后的指标）
# 实际方差可以用 adata.var['variances'] 拿到
adata.var['variances'] = np.var(adata.X.toarray(), axis=0)

为什么我不用纯方差？你看，均值高的基因天然方差大，但未必有价值。Scanpy 的离散度指标其实是对数归一化后的方差与均值关系拟合残差——效果更稳健。

我的习惯：保留前2000～5000个高变异基因。对于10x数据，这个数量通常能覆盖大部分生物学信号。如果后续做聚类，建议保持2000～3000个；如果做轨迹推断，可以放宽到5000个。

基于均值的过滤：去掉“沉默”基因

均值过滤解决另一个问题：表达量普遍极低的基因。很多基因在大多数细胞里计数为零，偶尔一两个细胞有微弱表达。这种稀疏特征会导致模型不稳定。

我早期做过一个错误的尝试：用全局均值>0.1的阈值过滤，结果把一个稀有的B细胞亚群的特异标志物全都切掉了。后来学乖了——阈值要根据数据分布动态设定。

# 基于平均表达量过滤，保留均值 > 0.01 的基因
min_mean = 0.01
sc.pp.filter_genes(adata, min_cells=20)  # 先按细胞数过滤
# 或者手动
mean_expr = adata.X.mean(axis=0)
keep = np.array(mean_expr > min_mean).flatten()
adata = adata[:, keep].copy()
print(f"保留基因数: {adata.n_vars}")

你可以用直方图观察均值分布。数据稀疏时，大部分基因的均值集中在0附近。一个技巧：先用 sc.pp.filter_genes(adata, min_cells=3) 去掉在少于3个细胞中有表达的基因，再根据剩余的分布设阈值。

避坑指南：我曾经直接用全局均值>0.1过滤，结果丢失了所有低表达的转录因子基因。后来改为按细胞类型分别过滤，或者用“至少在X%的细胞中表达”这个条件，效果好了很多。

稀疏性处理：应对零膨胀数据

单细胞数据是个典型的稀疏矩阵——60～90% 的元素是零。这些零不全是“未表达”，也可能是技术性缺失（dropout）。特征选择时如果不处理稀疏性，零基因会大量涌入模型，干扰训练。

常见的稀疏性处理策略有三种：

策略	方法	适用场景
硬阈值法	要求基因在至少k个细胞中表达	快速初筛，通用
矩阵补全	用imputation估算缺失值	数据质量差、期望保留更多特征
二值化+计数过滤	表达>0视为1，再统计细胞数	处理dropout时保留定性信息

我个人最常用的是第一种。AnnData 里直接用：

sc.pp.filter_genes(adata, min_cells=10)

min_cells=10 表示至少在10个细胞中有非零表达。这个值取决于你的细胞总数：对于1000个细胞，可以设20～50；对于10000个细胞，100～200比较合理。我一般设为细胞总数的0.5%～2%。

注意：min_cells 设太高会丢失罕见群体的标志基因。我有一个真实案例——某个细胞类型只占3%，但它的 Marker 基因在6%的细胞中表达，我一刀切设10%就给砍没了。后来我改为两步：先按全局过滤，再单独保留那些在很小但生物学确定的群里有表达的基因。

三步走的特征选择流程

把方差过滤、均值过滤、稀疏性处理串起来，就成了一个稳定的初筛 pipeline。下面这个流程图展示了我的标准做法：

流程顺序有讲究。先做稀疏性过滤，避免零值干扰均值和方差的计算。第二步均值过滤，剔除大部分“沉默基因”。最后方差过滤，精挑最有信息量的特征。三步走完，维度通常能从 2 万降到 2000～5000，同时保留 90% 以上的聚类结构。

实战：基于 sc.pp.highly_variable_genes 的一步法

如果你不想手工分步，Scanpy 提供了 sc.pp.highly_variable_genes，它内部整合了均值、离散度、以及一个可选的 min_mean / max_mean / min_disp 阈值。我最常用的调用是：

sc.pp.highly_variable_genes(
    adata,
    n_top_genes=3000,
    flavor='seurat_v3',    # 使用 Seurat 的方差稳定变换方法
    subset=True,
    min_mean=0.0125,
    max_mean=3,
    min_disp=0.5
)
print(adata.shape)
# 输出类似：(1000 cells, 3000 genes)

这个接口把三步压缩成了一行，但内部逻辑和手动流程完全一致。我建议初学者手动过一遍，再去用高级接口。否则哪天换了个数据集，发现 n_top_genes=3000 选出了一大堆低表达的基因——你都不知道该调哪个参数。

我个人的偏好：对于10x染色质可及性数据（如scATAC-seq），稀疏性远高于RNA-seq。我的min_cells通常设为细胞总数的5%。对于RNA-seq，0.5%~1%就够了。你可以先跑一次全pipeline，检查过滤后的基因数，假如剩余基因不到1000，说明过滤太狠；超过8000，说明还有大量噪声。

特征选择不是一步到位的艺术。好在我把常见的坑都踩了一遍，上面这些步骤基本能应对 80% 的数据场景。记住一个原则：宁少勿滥，但少到丢失生物学信号就得不偿失了。你可以在不同阈值下跑一轮聚类，观察聚类稳定性和 marker 基因的保留情况——这是我validation阶段的标配。