虚拟扰动分析实战：从Geneformer到SHAP的单细胞基因筛选

1. 先搞清楚为什么我们需要“虚拟扰动”这种实验大约两年前我第一次在一个单细胞转录组项目里遇到“必须做基因扰动实验但实验室条件根本不允许”的困境。导师给的课题方向是研究某种罕见疾病相关基因的网络调控机制传统的做法是先敲除基因再看细胞转录组的变化。但问题是这批样本本身已经非常稀缺细胞培养条件极其苛刻更不用说做 CRISPR 筛选了——成本和时间都是我们根本承担不起的。就在那个阶段我开始接触一种从命名到思路都颠覆传统实验逻辑的实践方式Geneformer 虚拟扰动分析。简而言之它不是物理意义上“敲掉”某个基因而是通过一个大规模预训练的 AI 模型在计算环境中模拟基因被扰动之后的转录组状态。你不需要碰移液枪不需要养细胞不需要等几周甚至几个月的实验周期你只需要把已经有的单细胞转录组数据输入模型让模型在计算空间里对该基因做一次“虚拟敲除”然后观察下游基因表达的变化。这个方向之所以值得讨论并不是因为它已经能完全替代湿实验。从我的实际体感看它更像是把一个原本“先做实验再看结果”的过程变成了“先在计算空间里做几十种假设验证再挑出最值得做的几只去实验室落地”。本质上它是把实验假设的筛选过程大幅前置了。当然虚拟扰动也有它明确的能力边界。我先说结论**这套方法真正解决的不是把基因敲除实验省掉而是把“找该敲哪个基因”的成本降下来。**接下来我会从原理、完整实操流程、SHAP 分析结合以及常见的坑这四个维度展开把一套可用到单细胞项目里的工作框架整理出来。2. 从单细胞数据到虚拟基因敲除这些概念必须串起来理解2.1 Geneformer 的基本定位一种单细胞转录组的 Transformer很多人第一次看到 Geneformer 这个名字容易把它理解成另一个普通的深度学习分类模型。实际上它更像是把自然语言处理里被验证过的预训练范式迁移到了基因表达数据上。在自然语言处理里Transformer 的输入通常是一段句子句子里的“词”是离散的文本符号。Geneformer 借用了这个结构但把“词”换成了基因。一个细胞中的基因表达数据通过排序和编码变成一串携带表达量信息的“基因序列”。模型通过海量单细胞转录组数据的预训练学会了基因之间的共表达关系、调控关系和上下级联关系。正因为有这个预训练步骤Geneformer 在面临一个并不包含在训练标签里的新任务时仍然可以利用已有知识做预测。这个过程和“看图识物”有一点类似。预训练阶段模型见过大量的“细胞状态”它对哪些基因在什么状态下会同涨同跌、哪些基因是上游调控节点、哪些是下游效应基因已经形成了一个内部表示。到了虚拟扰动任务里把一个目标基因标记为被扰动状态模型就基于它对细胞网络的整体理解去重构在这个扰动假设下其他基因的表达趋势。2.2 单次跑通不等于理解原理先理清要经历的四个步骤很多教程喜欢直接把代码放出来跑通一个示例就觉得事情结束了。但如果你要真的把这个方法用在自己的数据上没有理解它的流程骨架后续报错了都不知道该修哪一层。从实践看虚拟扰动分析的主干步骤可以拆成四个阶段数据准备把原始单细胞数据从表达矩阵形式处理成 Geneformer 能够接受的输入格式。这个阶段涉及基因符号标准化、表达量排序、Token 化编码。模型加载与预训练权重准备确认当前环境里的模型版本加载对应的预训练权重避免后续因为权重新旧不一致导致指标不可比。扰动模拟确定要虚拟敲除的目标基因构造扰动状态下的输入序列。这一步并不是简单删掉该基因而是要在一个更合理的转录组语境下模拟它不表达之后的影响。输出解读对比扰动前后模型输出的表达状态计算受影响基因列表通常配合 SHAP 分析归因找到哪些基因在网络中贡献最高。在这四个步骤中最容易出问题的其实不是模型代码本身而是第一步的数据准备工作。原因很简单单细胞数据的基因命名在不同版本、不同数据库之间常有差异如果一开始用了一组过时的基因名后续所有分析都会安静地跑到错误的结果上而不会报错提示你。2.3 “虚拟敲低”和“真正敲除”到底是不是一回事这是很多新入坑的朋友最容易被绕晕的概念问题。基因敲除的本意是让目标基因完全不表达但虚拟扰动在实现上更接近“部分扰动”或“网络状态重构”。因为模型学到的是数据分布和共表达规律它不是生物学意义上的 Cas9它并不知道把一个基因完全移除后细胞内是否会发生补偿机制、反馈调节或者细胞死亡。因此虚拟扰动更适合用来回答的问题是“如果这个基因的功能被干扰最可能在转录组层面看到哪些下游变化”而不适合回答“这个基因突变后细胞在表型层面会不会死亡。”这也决定了适用边界。在初期做候选基因筛选时虚拟扰动可以帮你把一百个候选基因收敛到十个高价值目标。但是到了验证阶段仍然需要回到湿实验里用真实的敲除或敲低去确认。它不是替代而是排序。一个很重要的实操提醒在对目标基因做扰动之前先观察一下这个基因在你自己的数据集里是否真的有足够的表达量。如果它在大多数细胞里本身就是低表达或零表达虚拟扰动结果可能会非常难以解释因为模型能捕捉到的信号本身就极弱。3. 一套可落地的虚拟扰动分析流程从环境准备到结果解读3.1 环境与前置依赖先确认版本再开始复制代码不管你使用的是本地 GPU 工作站还是云服务器环境准备始终是一个绕不开的环节。Geneformer 基于 PyTorch 和 Hugging Face Transformers 生态构建所以依赖版本是否匹配直接决定了你后面会不会遇到一堆莫明其妙的报错。我的习惯是先建立独立的虚拟环境避免和项目里的其他深度学习依赖互相污染。基础依赖通常包括 Python 3.9 或 3.10、PyTorch、Transformers、Datasets、Scikit-learn、Pandas 和 SHAP 相关的依赖包。下面这段是常见的基础安装示例结构实际版本号建议结合你的 GPU 驱动和 CUDA 版本来确定conda create -n geneformer python3.10 conda activate geneformer pip install torch pip install transformers datasets scikit-learn pandas pip install shap这里多说一句不要盲目追求最新版本。Geneformer 的预训练权重可能是在固定版本组合下导出的如果你把 Transformers 升级到非常新的主版本有可能出现权重名称不匹配或者 API 不兼容的问题。建议先查清当前加载权重要求的版本范围再创建环境。3.2 单细胞数据处理这一层决定了虚拟扰动的质量下限在实际项目里我见过不少用户直接拿一个 Seurat 对象或者 Scanpy 的 AnnData 对象就准备往模型里怼。这种做法在少数演示场景里也许能跑动但一旦数据里混入了不规范的基因名、包含了大段低质量细胞或者存在明显的批次效应模型输出的扰动结果就会非常不稳定。更稳妥的做法是先把原始表达矩阵标准化到 Geneformer 期望的输入格式。标准流程通常包括以下几个环节基因符号标准化确保所有基因名都使用同一套标准比如 HGNC Symbol避免出现“同一基因不同写法”被模型当成两个基因。过滤低质量细胞按基因数、UMI 数和线粒体基因比例过滤掉质量较差的细胞这一步在虚拟扰动前非常重要因为低质量细胞会引入大量噪声表达。表达量排序Geneformer 的输入核心是把每个细胞里的基因按表达量降序排列每个基因对应一个排序后的 Token 序列。这里的 Token 化过程会和普通 NLP Token 化有一些差异核心在于保留排序信息而非原始数据逐位拼接。用代码表示大致是这个思路def sort_genes_by_expression(cell_expression): sorted_genes [ gene for gene, _ in sorted( cell_expression.items(), keylambda item: item[1], reverseTrue ) ] return sorted_genes当然实际处理时还会涉及 Padding、截断长度等细节但核心逻辑就是“排序后编码”。把这一步做扎实模型后续看到的每个细胞都是一个结构化的、有序的基因序列。3.3 虚拟扰动与随机对照不能只做一次要形成扰动-对照组虚拟扰动实验必须有一个对照组否则你很难判断模型输出的变化到底来自目标基因扰动还是来自模型自身的不确定性。比较常见的做法是构造两组输入对照组对同一个细胞的基因序列做多次随机扰动但不是目标基因本身而是随机选取等量的无关基因进行同样的 Token 替换。目标扰动组对目标基因所在的 Token 位置做扰动替换保持其他条件不变。随后比较两组输出转录组状态的差异用差异来推断目标基因特异的扰动效应。这里的核心思想是控制变量只有目标基因不同的两组输入其输出差异才可以归因到目标基因。实际执行时随机扰动次数不宜太少。从我的经验看每个目标基因至少重复十次以上再对结果取平均能够得到相对稳定的扰动响应轮廓。如果资源允许重复二十次会更稳。3.4 结果解读不要只看差异基因列表还要看功能方向完成虚拟扰动之后最常见的操作就是对扰动前后的差异表达基因做富集分析。这个环节虽然看起来简单但有一个很容易掉进去的坑富集分析使用的背景基因集必须和你的数据来源一致。如果你拿小鼠的数据做扰动但使用了人的 GO 注释集那结果只能用来参考不能直接作为生物学结论。更关键的是不要只看单一时间点的扰动效应。虚拟扰动机制上类似一个静态快照它给出的是模型在给定状态下的推断而不是细胞培养若干小时后的动态响应。所以解释结果时把重点放在“哪些通路发生了方向性变化”上而不是纠结于“这个基因在第 3 天表达量会变成多少倍”。我自己通常会把扰动结果整理成一张表至少包含以下字段字段名含义目标基因被虚拟敲除的基因符号下游基因扰动后表达变化最显著的基因变化方向上调或下调差异分数扰动前后模型输出差异度量所属通路基于功能注释或先验知识补充置信度多次随机扰动后的稳定性这张表的价值在于它能把模型的潜空间输出变成可供生物学判断的候选清单。后续做湿实验设计时我一直是从这张表里挑基因验证而不是重新面对成百上千个差异基因。4. SHAP 在虚拟扰动里的角色从“扰动结果”到“归因解释”4.1 为什么单看扰动结果还不够需要 SHAP 来解释如果说虚拟扰动回答的是“如果该基因被破坏细胞状态会怎样变化”那 SHAP 回答的则是“在这个预测过程中到底是哪些基因在发挥主要推动力”。两者互不替代但放在一起看能帮我们构建一个更完整的叙事。很多做单细胞分析的人对 SHAP 的第一反应是“这不就是机器学习里的可解释性工具吗和虚拟扰动有什么关系”确实有关系。在机器学习模型里SHAP 基于博弈论中的 Shapley 值可以把模型对某个样本的预测结果分解到每一个输入特征上给出该特征对预测结果的贡献方向和大小。放到 Geneformer 虚拟扰动场景里输入特征是基因 Token 序列中的不同位置。当我们用模型对比扰动前后状态时可以在输出层对受影响的下游基因计算 SHAP 值从而知道哪些位置的变化在模型看来影响最大。这里的位置可以对应到具体基因。简单地说SHAP 形成了一个“可解释层”它不改变扰动结果但能让结果从“一堆基因变化”变成“一个归因清单”。4.2 SHAP 分析和扰动分析的正确协作顺序我第一次做这个流程时犯过一个错误先跑扰动然后再把 SHAP 套在扰动结果上。实际上更好的顺序应该是分成两轮。第一轮先做虚拟扰动筛选得到受影响的下游基因列表。这时候用的模型可以是 Geneformer也可以是在你的数据集上做过微调的版本。第二轮基于这个候选基因列表使用 SHAP 对模型做归因分析。目标不是重新跑一遍所有基因而是要看目标基因在下游基因预测中的贡献分布。具体操作上选定下游基因集合后对模型输入做 SHAP 分解得到每个输入 Token 的贡献值。贡献值大且方向一致的基因大概率是网络中的关键节点。4.3 SHAP 图怎么画才能既美观又有信息量搜索热词里反复出现“shap 图怎么画”说明大家在可视化阶段普遍卡住了。SHAP 图的意义不只是好看而是让审稿人或者同事一眼看出哪些基因贡献最大、方向如何。常用的几种图各有各的优点我按优先级介绍。第一推荐的是Shapley 值条形图。它把每个基因的 Shapley 值取平均绝对值后按降序排列能快速展示出影响最大的基因列表。在候选基因很多时这张图的信息密度最高。import shap explainer shap.TreeExplainer(model) shap_values explainer.shap_values(X_sample) shap.summary_plot(shap_values, X_sample, feature_namesgene_names)第二推荐的是蜜蜂图也就是散点图。它可以同时展示每个样本的 SHAP 值分布和基因表达量高低对贡献的影响。横向坐标是 SHAP 值纵向是基因颜色深浅表示原始表达量。这张图适合展示“该基因是主要在低表达细胞中发挥贡献还是在高表达细胞中发挥贡献”。第三推荐的是瀑布图。瀑布图适合单个样本的归因展示。如果有一个特别典型的目标细胞状态可以用它展示各个基因从基线到最终预测结果的贡献路径。这个图在论文的案例展示部分非常好用。画图的时候有两点经验第一先控制基因数量只展示贡献值绝对值排前十五到三十的基因不然图会糊成一片第二不要只画一张总图最好按细胞类型或样本分组各画一张因为不同细胞类型对同一目标基因扰动的响应机制可能差异明显。shap.summary_plot( shap_values, X_sample, feature_namesgene_names, max_display20, showFalse )4.4 SHAP 结合虚拟扰动的常见误区很重要的一点SHAP 反映的是模型内部的预测归因不是真实的因果效应。一个基因在 SHAP 里具有高贡献只代表模型在推断该状态时高度依赖这个特征并不等于它在细胞内就是真正的上游调控因子。这句话我想用最大的字体刻在虚拟扰动分析的操作文档里。尤其是当你对模型做过微调时SHAP 的值会受训练数据分布和标签定义影响。如果训练数据里某种细胞类型比例过高SHAP 结果就可能偏向于该类别的基因网络。因此在写论文或者做技术汇报时一定要把“模型依据”和“生物学证据”分开表述。模型依据用 SHAP 归因生物学证据用实验验证或已有文献支撑。5. 从单次 Demo 走向项目级应用必须补齐的四块拼图5.1 日志与可追踪性每次扰动结果都要能复现如果是自己学习跑一两次 demo代码不写日志问题不大。但一旦进入真实项目哪怕只是三个细胞亚群的扰动比较都要有完整的日志记录。我建议至少给每次扰动分析记录以下信息Geneformer 模型版本和预训练权重来源数据预处理的时间、基因映射标准版本目标基因名称及 ID 映射结果随机扰动次数和随机种子输出文件路径和运行时间戳SHAP 计算时的样本选择逻辑这些信息能保证你在一个月后重新打开结果时还能准确还原当时的分析条件。很多人高估了自己的记忆能力低估了结果复现的难度这是项目推进中的隐形风险。5.2 批量策略逐基因扰动比一次扰动大批基因更有用真实场景下我们往往有多个候选基因需要评估。但不需要把所有基因塞进一次扰动里。更合理的方式是设计一个循环对每个目标基因单独做扰动、记录结果最后整合成一个大表。每个目标基因单独做有几个明显好处第一出了问题能快速定位到具体基因不会因为某个异常基因污染整批结果第二能分别计算每个基因的扰动稳定性为后续置信度排序提供数据第三更方便和 SHAP 归因配合因为每个基因都可以单独生成一组图。下面是一个循环处理的示例结构target_genes [GENE1, GENE2, GENE3] results [] for gene in target_genes: perturbed_states run_virtual_perturbation(cell_data, target_genegene) summary summarize_perturbation(perturbed_states) results.append(summary)5.3 权限与资源管理GPU 不是无限使用的Geneformer 虽然通过预训练降低了大量下游任务的计算成本但虚拟扰动需要在每个样本上做多次模型推理如果数据量是数万个细胞计算开销会迅速上升。实际操作时我一般先取一个小型子集做扰动验证确认流程无误后再扩展到全量数据。同时要关注 GPU 显存占用。如果你的模型加载后显存过高可以尝试把批量大小拉低或者按细胞亚群分批处理不要在同一时间内把所有细胞塞进模型。如果只是做 SHAP 归因有些环节甚至可以考虑在 CPU 上运行因为速度虽然慢一点但可以避免 GPU 资源被长时间占用导致其他任务排队。5.4 结果验证闭环虚拟结果最终还是要回到真实数据里我见过一些项目虚拟扰动结果特别漂亮但一旦和真实实验数据对比对上的比例却不理想。原因通常是数据偏见虚拟扰动模型依赖原始训练数据分布如果训练数据里没有某些罕见细胞类型或疾病状态它对这些状态的推演就是在用自己的幻想来填空。所以落地时建议至少做一个验证闭环第一层虚拟扰动是否能在已有真实扰动数据集里复现已知的调控关系。第二层在与目标疾病或状态相关的患者数据里虚拟扰动信号是否和临床趋势一致。第三层挑选排名靠前的下游基因回到已有文献里确认是否有独立证据支持。这一层验证做完后虚拟扰动的结论才算真正从“模型输出”变成了“可用假设”。6. 一个更清醒的判断虚拟扰动是假设引擎不是实验替代品虚拟扰动分析最吸引人的地方是它提供了一种低成本、高通量、可重复的“先验假设生成”能力。它可以把过去需要几个月才能完成的候选基因筛选压缩到几天甚至几小时。但它也有非常明确的边界。它不是晶体结构预测也不是蛋白-蛋白相互作用模拟它本质上是一个基于大规模单细胞表达数据的表示学习模型。它在基因调控网络推断、细胞状态转换假设、候选靶点排序这些任务上表现出了独特的价值但它不能回答蛋白翻译后修饰、代谢通路动态变化、细胞间微环境相互作用等超出转录组信息范围的问题。我经常把虚拟扰动理解成导航地图里的智能推荐路径。地图会根据历史数据推荐几条最可能到达的路线但它不会替你把车开完。最后时刻的你仍然需要自己上路、遇到真实路口、处理真实路况。所以在项目规划里虚拟扰动应该占据的位置是早期筛选阶段。它帮你快速收敛候选基因建立优先级清单减少无效实验。到了验证阶段它让位于更经典的扰动表达实验。这个顺序才是最经济、最科学的用法。如果你正在接触单细胞数据、基因调控网络或候选靶点筛选我建议你不要把精力全部放在调整模型参数上。先找一小套已经有部分真实扰动结果的数据作为对照把“虚拟 vs 真实”的一致性做出来你才算真正掌握这个方法。因为只有当你亲眼看到虚拟扰动结果在多大程度上贴近真实实验时你才能对这套模型产生合理的使用直觉而不是盲目信任或者随意否定。下一步优先做的事情不是下载更多数据也不是换一个更大规模的模型而是先把你手里最熟悉的一个基因完整地走一遍虚拟扰动到 SHAP 归因的流程再把输出结果和已知文献对照一遍。这个方法本身不难难的是在每一步都保持对结论边界的分寸感。