Scanpy中HVG函数返回基因数超限问题的分析与解决

问题背景

在单细胞RNA测序数据分析中，筛选高变基因(HVG)是一个关键预处理步骤。Scanpy作为Python生态中广泛使用的单细胞分析工具，其highly_variable_genes函数用于识别数据集中变异程度最高的基因。然而，用户在使用过程中发现，当指定n_top_genes=10000时，实际返回的高变基因数量可能达到13355个，远超预期。

问题根源分析

经过深入代码审查，发现问题出在离散度计算和阈值确定的逻辑上。具体来说：

1. 在计算归一化离散度(dispersion_norm)时，部分基因可能产生NaN值
2. 当前实现中，首先过滤掉NaN值计算离散度阈值(disp_cut_off)
3. 随后将剩余基因的NaN值转换为0进行比较
4. 当disp_cut_off为负值时，所有被转换为0的NaN值都会通过阈值检查

这种不一致的处理方式导致最终筛选出的高变基因数量超出预期。

技术细节

核心问题出现在以下代码逻辑中：

disp_cut_off = _nth_highest(dispersion_norm, n_top_genes)  # 计算时不考虑NaN
return np.nan_to_num(dispersion_norm) >= disp_cut_off  # 比较时将NaN转为0

当数据集中存在大量低表达或零表达基因时，它们的离散度计算可能产生NaN值。这些基因在阈值计算阶段被排除，但在最终比较阶段又被包含进来，导致结果数量膨胀。

解决方案

经过讨论，社区确定了两种可能的修复方案：

1. 将NaN值转换为负无穷大：np.nan_to_num(dispersion_norm, nan=-np.inf)
2. 在计算阈值前就将NaN值转换为0

最终采用了第一种方案，因为它更符合数学逻辑：NaN值代表的"不可计算"状态应该被视为最不可能成为高变基因的情况。

测试验证

为了验证修复效果，开发者设计了专门的测试用例：

def test_hvg_gene_count():
    adata = sc.datasets.pbmc3k()
    sc.pp.normalize_total(adata, target_sum=10000)
    sc.pp.log1p(adata)
    sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=10000)
    assert adata.var["highly_variable"].sum() == 10000

该测试模拟了实际使用场景，确保修复后函数能准确返回指定数量的高变基因。

相关讨论

在问题排查过程中，开发者还发现了另一个相关现象：当多个基因具有相同的离散度值时，函数会返回多于指定数量的基因。这与统计学中处理并列排名的惯例一致，被认为是合理行为而非缺陷。

总结

本次问题修复确保了Scanpy中高变基因筛选功能的准确性，特别是在处理包含大量低表达基因的数据集时。这一改进使得分析结果更加可靠，为下游分析提供了更好的基础。对于单细胞数据分析人员来说，理解这一修复有助于更准确地解释分析结果，避免因工具行为异常导致的错误结论。