Scanpy中HVG函数返回基因数多于预期的问题分析

问题背景

在单细胞RNA测序数据分析中，筛选高变异基因(HVG)是一个关键预处理步骤。Scanpy作为Python生态中广泛使用的单细胞分析工具，其highly_variable_genes函数用于执行这一任务。然而，在某些情况下，该函数会返回比用户指定数量更多的基因，这可能导致下游分析出现意外结果。

问题现象

当用户指定n_top_genes=10000时，函数实际返回了13355个基因，远超过预期数量。这种情况在特定数据集上可复现，特别是当数据中存在大量零表达或恒定表达的基因时。

技术原因分析

问题根源在于highly_variable_genes函数的实现逻辑：

1. 归一化离散度计算：函数首先计算基因表达值的归一化离散度(dispersion_norm)
2. 截断值确定：通过_nth_highest函数找到第n_top_genes高的离散度值作为截断值
3. 基因筛选：选择离散度大于等于截断值的基因

关键问题出现在以下两个环节：

1. NaN值处理不一致：_nth_highest函数首先移除NaN值再计算，而后续比较时NaN被转换为0
2. 离散度相同处理：当多个基因具有相同的离散度值时，函数会保留所有这些基因

解决方案

针对NaN值处理问题，开发团队提出了以下修复方案：

# 原代码
return np.nan_to_num(dispersion_norm) >= disp_cut_off

# 修复代码
return np.nan_to_num(dispersion_norm, nan=-np.inf) >= disp_cut_off

这一修改确保NaN值不会被误认为满足条件，因为负无穷大永远不会大于等于任何实际计算得到的截断值。

测试验证

为了验证修复效果，开发团队设计了以下测试用例：

def test_no_filter_genes(flavor):
    """测试即使在数据包含0表达基因时，n_top_genes参数也能被正确遵守"""
    adata = sc.datasets.pbmc3k()
    means, _ = _get_mean_var(adata.X)
    assert (means == 0).any()
    sc.pp.normalize_total(adata, target_sum=10000)
    sc.pp.log1p(adata)
    sc.pp.highly_variable_genes(adata, flavor=flavor, n_top_genes=10000)
    assert adata.var["highly_variable"].sum() == 10000

该测试验证了在包含零表达基因的数据集上，修复后的函数能准确返回用户指定数量的高变异基因。

相关讨论

在问题分析过程中，还发现了另一个相关现象：当多个基因具有相同的离散度值时，函数会返回多于指定数量的基因。这与Seurat等工具的行为不同，后者会应用额外的阈值限制。Scanpy团队认为这种行为是合理的，因为相同离散度的基因应当被视为同等重要。

总结

Scanpy的highly_variable_genes函数在处理包含NaN值或相同离散度基因时存在返回基因数多于预期的问题。通过将NaN值转换为负无穷大而非零，可以解决主要问题。对于相同离散度基因的处理，则被视为合理行为而非缺陷。这一修复确保了函数行为的可预测性，同时保持了生物学合理性。