编辑推荐：
本文提出一种基于机器学习和归因模型的反事实推断工具 CIMLA，用于识别生物条件间基因调控网络（GRN）的差异。通过理论推导、模拟数据和真实数据验证，CIMLA 在处理混杂变量时更稳健、准确，为研究阿尔茨海默病（AD）等疾病的调控机制提供了新方法。
一、研究背景
在生物信息学领域，从组学数据中提取因果分子关系是一个关键挑战。基因调控网络（GRN）作为理解基因表达变化分子机制的重要框架，其推断一直是研究热点。然而，目前从转录组数据推断 GRN 的主要生物信息学方法，在将特征重要性分数解释为转录因子（TF）对基因的因果影响方面存在不确定性，并且生物学家和生物信息学家对 GRN 中调控关系的理解存在概念差距。此外，定量定义 GRN 边缘缺乏精确性，从表达数据单独推断 GRN 时，TF 与基因的调控关系难以明确界定。
二、研究内容与方法
- 定义调控关系并建立与机器学习的联系
- 研究人员使用因果推理语言，仅通过表达干预和测量，正式定义了 TF - 基因调控关系。假设存在m个候选 TF，用Xt表示其表达水平，Yg表示目标基因g的表达。为量化 TF t对目标基因g的因果影响，提出 “局部治疗效应”（local treatment effect，LTE），公式为LTEt,g(x,x^t)=E[Yg∣do(Xt=xt),do(Xj∈M?{t}=xj∈M?{t})]?E[Yg∣do(Xt=x^t),do(Xj∈M?{t}=xj∈M?{t})] ，其中do是 Pearl 的 “do - 算子”，代表干预。但由于 LTE 难以从观测数据估计，研究人员对其进行了重新定义和调整，引入了αt,g(x)作为可从观测数据估计的 TF - 基因调控关系的度量。
- 研究发现，在某些假设下，训练一个能从X预测Y的机器学习模型，并使用 SHapley Additive exPlanations（SHAP）特征归因模型量化协变量Xt对模型输出的贡献，得到的 SHAP 分数等于αt,g(x)。这建立了上述定义的调控关系量化与基于机器学习推断 GRN 边缘策略之间的联系。
- CIMLA 工具概述
- CIMLA 是一种用于识别两个转录组数据集之间 TF - 基因调控关系变化的计算工具。其输入为两个条件下多个样本的表达谱，目标是检测特定 TF 与目标基因之间的直接调控关系在两个群体中是否发生变化，即 “差异调控关系”。
- CIMLA 主要包括三个模块：
- ML 模块：分别对两个群体的转录组数据训练机器学习模型，以预测目标基因的表达。当前实现依赖随机森林（RF）和带随机失活（dropout）的神经网络（NNs）来处理协变量之间的多重共线性。
- 解释模块：将训练好的模型和数据传递到该模块，使用 TreeSHAP 和 DeepSHAP 分别估计 RF 和 NN 模型中每个协变量对预测结果的局部贡献，即 SHAP 值。这一步基于理论洞察，将 SHAP 分数与反映 TF 对基因影响的因果量的可估计代理联系起来。
- 聚合模块：计算两个模型中局部贡献分数的差异，并通过 “均方根” 函数在所有样本上进行聚合，得到 CIMLA 分数。该分数表示协变量与结果之间因果关联在两个群体中变化的程度。
- 模拟转录组数据的基准测试
- 由于缺乏评估差异基因调控网络（dGRN）推断方法的金标准数据集，研究使用 SERGIO 模拟器生成代表两种不同条件的单细胞表达数据集，模拟具有不同共享边比例的 GRN。这些数据集分为 “低混杂” 和 “高混杂” 两种设置，在低混杂设置中，使用可比的主调节因子（MRs）表达谱；在高混杂设置中，设置 MRs 表达谱不同，以模拟未观察到的混杂因素。
- 使用 CIMLA 和其他方法对模拟数据进行 dGRN 边缘推断，与六种基于共表达的方法、GENIE3 - diff、BoostDiff 和 DoubleML - diff 等方法进行比较。评估指标包括受试者工作特征曲线下面积（AUROC）和归一化精度召回曲线下面积（AUPRC）。
- CIMLA 在相关性中优先考虑因果关系
- 研究定义了 “delta - 相关性” 来衡量 TF - 基因对在病例和对照组中的相关性差异。在低混杂设置中，真正的差异对和非差异对的 delta - 相关性分布差异明显；在高混杂设置中，两者分布更相似，且非差异对的 delta - 相关性更高，这表明混杂因素会导致基于相关性的方法出错。
- 比较 CIMLA 和基于相关性的z - score - S 方法对差异 TF - 基因对的评分，发现 CIMLA 分数受 delta - 相关性的影响较小，在高混杂设置中，CIMLA 能为 delta - 相关性较小的差异对分配更高分数，具有更高的真阳性率（TPR）和更低的假阳性率（FPR）。
- 噪声模拟数据的评估
- 考虑到真实单细胞 RNA 测序数据存在技术噪声，研究在具有 dropout 噪声的合成单细胞表达数据集上重复部分评估。在高混杂设置下，对一个测试添加不同水平的 dropout，并使用 MAGIC 对缺失值进行插补。除 RF 外，还测试了全连接 NN 作为 CIMLA 的基础模型，并使用 DeepSHAP 计算 SHAP 分数。
- 结果表明，随着噪声水平增加，dGRN 推断性能下降，但 CIMLA - RF 和 CIMLA - NN 在几乎所有 dropout 水平上均优于基于共表达的方法。在最高 dropout 水平下，CIMLA - NN 表现最佳。此外，发现 CIMLA - RF 和 CIMLA - NN 分数具有一定互补性，结合两者的 CIMLA - MeanRank 方法在不同噪声水平下均表现出色。
- 阿尔茨海默病（AD）的差异 GRN 案例研究
- 使用 CIMLA 分析先前发表的 AD 患者和非 AD 患者前额叶皮层的单细胞核 RNA 测序（snRNA - seq）数据集。对数据进行预处理，包括数据标准化、缺失值插补和基因筛选。使用 CIMLA - RF 和 CIMLA - NN 分别识别两组之间的差异 TF - 基因对。
- 评估 ML 模型的准确性，发现 1803（RF）和 2079（NN）个基因可被可靠建模。通过随机打乱细胞组标签获得背景分布，评估预测的差异调节因子的统计显著性。限制 dGRN 为差异表达基因，并确定可能在 AD 相关失调中起重要作用的枢纽 TF 和高靶向基因。
- CIMLA 揭示 CREB3 和 NEUROD6 为 AD 的潜在关键调节因子
- 将 CIMLA 的结果与已发表的人类大脑 GRN（PsychGRN）整合，以丰富直接调控关系。从 CIMLA - RF 和 CIMLA - NN 的 dGRN 中取前 10% 评分的 TF - 基因边的并集，得到组合 dGRN，再与 PsychGRN 取交集，得到高置信度 dGRN。
- 分析高置信度 dGRN 中 TF 的调控子，发现 CREB3（环磷腺苷 5′ - 单磷酸反应元件结合蛋白 3）具有最大的调控子，其参与了与 AD 相关的多种生物学过程；NEUROD6（神经元分化 6）的调控子在 dGRN 评分中位居第三，其下调是 AD 的生物标志物和认知下降的关键预测指标。此外，ELK1（ETS 转录因子 ELK1）的调控子也在 dGRN 评分中位居前列，其与神经元死亡和 AD 相关。同时，还发现了如 GATA3 等潜在的 AD 调节因子。
三、研究结论
- CIMLA 在现实合成数据集上进行 dGRN 推断时，优于现有基于线性和非线性模型的方法，尤其在处理强混杂因素时表现出色。这表明因果方法在 dGRN 推断中具有重要意义。
- CIMLA 能够识别在两个群体中 TF - 基因相关性差异较小的差异调控关系，相比基于共表达的方法更具优势。
- 通过 CIMLA 对 AD 的研究，发现 CREB3 和 NEUROD6 是 AD 的重要调节因子，并揭示了其他潜在的调节因子，为深入理解 AD 的调控机制提供了新线索。
- 尽管 CIMLA 在 dGRN 推断方面取得了较好的成果，但仍存在一些局限性。例如，CIMLA 使用机器学习模型进行反事实推断时，依赖于模型在未见过的数据分布上的可转移性假设，这在实际应用中可能存在问题。未来可通过使用能同时从两个群体的观测数据中学习的机器学习模型来改进 CIMLA，以更好地分离不同条件下的调节因子，并将该工具扩展到从更多条件进行 dGRN 推断。同时，对 SHAP 值的因果解释虽然取得了初步进展，但仍需进一步研究和完善，以更准确地捕捉真正的因果效应。
