1. 引言
历经数个世纪的生物学探索,人类对于生命基本单位——细胞在人体内的运作机制,仍面临着根本性的认知挑战。尽管现代生物技术已实现对基因组、转录组、蛋白质组及代谢组等分子构成及其空间分布的静态图谱解析,但是这些细胞的复杂性难以解析。这些复杂性主要来源于分子之间通过高度动态的相互作用实现细胞功能分化的内在机制[1]。作为解析遗传相互作用网络及其对细胞功能影响的核心手段,基因扰动技术的相关研究发现,正逐步成为现代生物技术创新与疾病治疗理论的重要科学基础[2]。
近年来,单细胞RNA测序技术(scRNA-seq)的突破性发展[3],并与高通量遗传微扰筛选[4]和CRISPR介导的基因编辑技术[5]实现联合应用,推动该领域实现技术革新并形成标准化研究范式[6]。这种技术融合使研究者能够在单细胞分辨率下系统性解析基因扰动对细胞状态的下游调控效应[7]。然而,人类细胞类型的多样性及其潜在扰动组合的计算复杂度随指数级增长,使得通过实验手段穷尽解析所有扰动条件下的细胞响应,在时间与人力层面均不可行。在此背景下,多基因扰动转录效应预测算法的研究范式经历了根本性转变:从早期的传统统计模型迭代发展为融合生物学先验知识的深度学习框架,并在可解释性和计算效率等关键维度取得突破性进展。尤为重要的是,单细胞多组学数据为深度学习技术提供了理想载体,使其能够在细胞分割[8] [9]、亚群聚类[10]及药物效应预测[11]等核心任务中实现高效建模。
在早期生成模型和线性传播阶段,研究主要依赖生成模型捕获基因表达分布特征。2019年Lotfollahi等人提出了scGen,首次将变分自编码器(VAE)引入单细胞扰动预测,通过潜在空间插值估计基因表达变化,在跨细胞类型预测中达到0.78的Pearson相关系数[12]。但该方法仅能处理单基因扰动,且无法解释调控机制。与此同时,Aibar等人提出了基于基因调控网络(GRN)的方法通过线性传播扰动效应,利用基因共表达矩阵构建调控关系[13],但在多基因组合预测中误差增加42%,且线性假设难以捕捉非线性遗传的相互作用。
计算复杂度随基因数量呈指数增长,进一步推动了知识增强型模型的兴起。2020年后,有研究者开始整合生物通路知识提升模型泛化能力。Roohani等人提出的GEARS模型利用基因共表达图编码转录协同性构建双路知识图谱、通过GO通路相似性图指导扰动传播[14]、通过图神经网络(GNN)聚合邻域信息,在131组双基因扰动实验中,MSE降低37%,对未见过的基因组合的预测精度提升53%。同期,Kamimoto等人提出了CellOracle模型,该模型将GRN与因果推理结合,通过扰动响应梯度反推调控路径,但仅限于单基因分析[15]。该模型的优点在于引入了Reactome等通路数据库构建注意力掩码、开发自聚焦损失函数、增强差异表达基因权重、方向感知损失纠正表达趋势偏差。然而,GNN的消息传递机制在处理层级通路结构时效率低下,且无法建模长程基因互作用。
这些方法中的大多数,包括scGEN,PerturbNet [15] [16]和CPA [17],主要预测化学处理的转录效应,而不是遗传干扰。它们通常将遗传扰动作为不同剂量的类似物,这不能充分捕获不同基因和特定基因之间的复杂相互作用,导致在该特定领域的性能有限。
直到Transformer架构的引入显著提升了跨基因关系的建模能力,其自注意力机制能够捕捉长程依赖关系。2024年,Bai等人提出AttentionPert,将注意力机制和图神经网络引入扰动分析,突破了传统线性模型或单一尺度方法的局限[18]。其提出的多尺度建模框架不仅适用于遗传扰动研究,还可迁移至药物组合预测、环境刺激响应等场景,为系统生物学提供了通用方法论。同年,Tan等人提出双流网络BioDSNN [19],采用变分自编码器构建数据驱动流,从而重构基因表达全局特征,并且通过注意力机制处理Reactome通路构建由生物先验知识构成的数据流。BioDSNN模型双流设计避免了单一数据或知识的局限性,同时利用数据统计规律与生物学逻辑,提升模型鲁棒性。相较于上一阶段,该阶段在基因预测的各个指标方面都有了长足的进步,不仅使用了混合架构实现数据驱动与知识驱动学习,还引入扩散模型提升生成质量,并在此基础上进一步开发不确定性量化指标(如GEARS的贝叶斯置信度)提高模型鲁棒性。但Transformer的二次计算复杂度以及长序列信息处理的缺陷限制了其在全基因组规模的应用[20]。
Mamba状态空间模型的提出为长序列处理带来革新。其选择性记忆机制通过输入依赖的参数化策略,在保持线性复杂度的同时实现94%的序列信息捕获率,能够挖掘更加复杂的数据特征[21]。在本文中,我们提出使用一种新型基于Mamba的双流框架模型,用于预测多基因扰动响应,以展现其在多基因扰动响应预测中具有独特优势。
2. 模型算法
本文构建的双流Mamba神经网络可简略分为三层,如图1所示,分别是输入层、编码层、输出层。下面我们将对该模型的三层进行详细介绍。
2.1. 输入层
输入层主要包含三部分数据,分别是基因表达数据、基因关系图谱数据、扰动关系图谱数据。这三部分数据构成了整个模型的输入。
基因表达量数据以及基因扰动数据来自Norman、RPE1、K562数据集。本文为了预测由扰动引起的转录结果,选择了这三个公开可用的扰动集。这三个数据集由单基因或多个基因组成。对于单基因,本文使用来自两种不同的遗传干扰筛选的数据,分别由1543 (RPE1细胞)和1092 (K562细胞)干扰组成,每个干扰测量超过170,000个细胞[22]。此外,Norman数据集[23]包含131个双基因扰动,那么我们可以评估不同扰动水平的性能。数据集总结见表1。
Figure 1. Diagram of the dual-stream mamba neural network model
图1. 基于双流Mamba神经网络模型流程图
Table 1. Dataset information
表1. 数据集信息
数据集 |
细胞数 |
基因数 |
扰动数 |
K562 |
162,751 |
5000 |
1092 |
RPE1 |
162,733 |
5000 |
1543 |
Norman |
91,205 |
5045 |
131 |
基因关系图谱数据来自Gears中对GO数据的处理结果。因为GO数据被用于测量基因之间的功能相似性,所以本文用来衡量基因之间的功能关系,以及探索更多的扰动信息。例如,
是基因u的Go途径集,一对基因u和v之间的Jaccard指数计算为:
(1)
Jaccard指数用于测量它们之间共享路径的比例。对每个基因,本文选择具有最高Jaccard指数的基因来构建扰动相似性图。然后,所有初始化的可能的基因扰动嵌入将被馈送到图卷积网络中以进行扩充。这样,在每次扰动嵌入时,都能综合相邻扰动的信息。
2.2. 编码层
编码层由两条数据流组成,分别是:基因表达量数据通过图卷积层及变分自编码器的数据流,以及由基因关系图谱数据和扰动关系图谱数据通过图卷积层以及Mamba编码层的数据流。接下来,本文将详细阐述每条数据流的构成。
两条数据流均选择了图卷积网络作为数据嵌入层,在捕获基因表达量数据特征时,使用了变分编码器。变分自编码器作为生成型网络具有多方面优势:能学习数据潜在分布,生成多样且连续变化的新样本;可将高维数据压缩到低维潜在空间,提取有语义的特征;结构灵活,易与其他技术结合;作为概率模型,既能解释数据生成过程,还能对生成结果进行不确定性估计。因此,本文采用变分编码器来捕获基因表达数据的不确定性和可变性。
变分自编码器能够测量潜在变量的学习分布(给定原始基因表达数据)与先验分布(通常为高斯先验分布)之间的偏差。具体来说,变分自编码器最大化证据的下界,称为证据下界(ELBO),因为边际似然是难以处理的,输入数据的边际似然的变分下界是变分自编码器的目标函数。通过对不同数据点的边际似然性求和获得边际似然性,如下所示:
(2)
本文采用了Mamba编码器来捕获基因关系图谱数据与扰动图谱数据的关联,Mamba作为新一种新兴的架构,更加聚焦于数据内部的相关性,使模型更加聚焦于重点信息。同时,Mamba与注意力机制相比,Mamba拥有更大的吞吐量以及更快的处理速度,且具有更好的通用性与适应性。标准状态空间模型(SSM)通常假设参数时不变,即认为状态转移矩阵等在序列处理时不改变。Mamba编码器通过引入时变参数,允许参数随着迭代步长进行变化。这种设计在保持线性复杂计算度的同时,为序列建模带来了更强的灵活性和表达能力,是Mamba编码器区别于其他编码器的关键创新点之一。时变参数的生成公式(3)所示。
(3)
其中
为激活函数。
然后,Mamba模型通过零阶保持法(Zero-Order Hold, ZOH)将连续状态空间模型(SSM)离散化。零阶保持法离散化支撑输入感知的参数生成,实现选择性机制,同时实现线性复杂度及全并行计算,成功突破长序列瓶颈。具体如公式(4)所示。
(4)
之后进行并行状态更新(如公式(5)所示),并使用并行扫描算法(Parallel Scan)避免递归计算。并行扫描算法具有颠覆性价值,在训练速度提升的同时,大幅度节约内存,是唯一支持动态参数化状态转移的并行方案,为Mamba的选择性机制提供计算基础。
(5)
(6)
最后通过门控机制以及残差融合,得到最后的输出:
(7)
本研究为确定Mamba模型的核心参数——状态维度(d_state)的最优取值,选取4、16、128三组数值,在REP1和K652数据集上开展实验。实验结果如图2所示:当状态维度值为16时,模型表现达到最优。
2.3. 输出层
在双流数据经过编码之后,再通过全连接层以求和的方式构建输出,即本模型预测的基因扰动数据。
2.4. 参数设置
在模型训练过程中,学习率调度采用基于gamma系数的衰减策略(gamma = 0.5)。当满足预设衰减条件(如验证集性能停滞)时,学习率将乘以gamma因子进行调整:初始学习率为1e−3,首次衰减后降至5e−4,二次衰减后进一步降至2.5e−4。这种渐进式衰减机制使模型在接近最优解时能够进行更精细的参数调整,有效提升收敛稳定性和泛化能力。
主要训练参数配置如下:学习率1e−3,训练轮次20,损失函数采用KL散度与均方误差的组合形式,以精确量化预测结果与真实标签的差异并引导优化方向。针对训练初期梯度波动剧烈的问题,设置较高的动量参数beta1 = 0.9,既抑制参数更新震荡风险,又保障后期精准收敛。为防止过拟合,引入L2正则化并设置权重衰减系数5e−4,通过约束模型复杂度,避免对训练噪声的过度学习。如图3所示,训练
Figure 2. Comparison chart of results from different state dimension settings. (A) K562 dataset; (B) RPE1 dataset
图2. 不同状态维度设置结果比较图。(A) K562数据集;(B) RPE1数据集
Figure 3. Training loss curves for different datasets. (A) Noman dataset; (B) K562 dataset; (C) RPE1 dataset
图3. 不同数据集的训练loss曲线图。(A) Norman数据集;(B) K562数据集;(C) RPE1数据集
过程呈现稳定且最终收敛。
3. 分析结果
如图4所示,为了评估扰动前后单细胞基因表达的差异,本文计算了对照组和实验组之间每个基因的MSE。由于绝大多数基因在未扰动和扰动状态之间没有显示出实质性的变化,因此本文只考虑前20个
Figure 4. Distribution of Mean Squared Error (MSE) in gene perturbation response prediction. (A) Noman dataset; (B) K562 dataset; (C) RPE1 dataset
图4. 基因扰动响应预测均方误差(MSE)分布图。(A) Norman数据集;(B) K562数据集;(C) RPE1数据集
差异表达的基因。图中“0/2 unseen”、“1/2 unseen”和“2/2 unseen”分别对应训练集中观测到两个、一个或零个扰动目标的实验条件。该分组设置精准模拟了基因扰动预测任务中数据稀疏性的挑战——训练时可见的扰动目标数量越少,模型需要填补的信息缺口越大,相应的预测难度呈指数级上升。从图4的实验结果来看,本文提出的模型在全组别数据上均展现出显著优势,预测性能显著超越所有基线模型。即使在“0/2 unseen”这一极端困难场景下,该模型的预测误差仍低于传统方法,充分验证了其在复杂基因扰动预测任务中的卓越泛化能力。
模型的优异表现源于两大核心创新设计。首先,双流网络架构打破了单一数据输入的局限性,通过构建并行信息通道,实现了基因表达数据与领域知识的深度融合:一条通路专注于从海量转录组数据中提取统计特征,另一条通路则将基因调控网络、功能注释等结构化知识编码为可学习的嵌入向量,二者通过注意力机制实现动态交互,确保模型能够兼顾数据驱动的模式识别与生物学逻辑约束。其次,新型Mamba编码器凭借其独特的线性复杂度架构与动态注意力机制,有效突破了传统Transformer在长序列建模中的效率瓶颈。该编码器不仅能捕捉基因表达数据中的长程依赖关系,还能根据不同扰动场景自适应地调整特征权重,在复杂基因调控网络中精准识别关键调控节点,从而显著提升模型的预测精度与鲁棒性。
本文采用一个不确定性分数来衡量模型预测对新扰动的置信度。采用高斯似然函数来对在扰动下的基因表达值进行建模。同时,添加一个额外的基因特定层来预测对数方差项(公式(8)),并通过修改的贝叶斯神经网络损失(公式(9))来修改模型的不确定性分数[24]。
(8)
(9)
Figure 5. Distribution of model uncertainty scores
图5. 模型不确定性评分分布图
如图5所示为本文提出的模型不确定性评分分布。分析该图可知,模型展现出较低的预测不确定性水平,且预测的基因表达差异与实际值之间存在显著相关性。这一结果表明,模型不仅能对基因扰动后的表达变化进行可靠预测,还能有效捕捉基因互作的复杂动态特征,为多基因扰动效应的精准解析提供了具有高可信度的分析工具。
在疾病研究领域,基因表达模式的动态变化蕴含着丰富信息,不仅能够为疾病诊断与预后评估提供极具价值的潜在生物标志物,还能通过筛选表达水平显著改变的基因,挖掘出治疗干预与药物研发的关键靶点。为深入探究单个基因在不同扰动场景下的响应机制,本研究针对三种典型基因组合的扰动结果展开预测分析,这些组合分别对应训练集可见性的不同层次:双可见基因组合(训练过程中两个基因均存在观测数据)、单可见基因组合(仅一个基因存在训练数据)以及零可见基因组合(两个基因均不存在训练数据)。以具体实验为例,ETS2与CEBPE基因组合属于双可见类型,训练期间可获取二者完整的扰动数据;FOSB与CEBPB组合则为单可见类型,其中CEBPB在训练阶段未经历实验扰动;而FOXL2与MEIS1组合处于零可见状态,训练过程中完全缺乏这两个基因的相关数据[25] [26]。
如图6所示,本研究提出的双流模型在三种复杂扰动场景下均展现出卓越性能,实现了对基因
注:红点表示本文模型的预测值,绿色虚线表示在未受干扰的控制条件下的平均基因表达。箱式图表示真实的基因表达情况。
Figure 6. Distribution of gene expression levels under dual-gene perturbation. (A) ETS2 and CEBPE; (B) POSB and CEBPB; (C) FOXL2 and MEIS1
图6. 双基因干扰下基因表达水平分布图。(A) ETS2和CEBPE;(B) POSB和CEBPB;(C) FOXL2和MEIS1
扰动后表达量的精准且稳定预测。该模型不仅在多基因交互作用的解析中表现优异,其预测结果与真实数据之间的高度相关性,更充分印证了模型在捕捉基因扰动复杂机制时的稳健性与可靠性。
4. 结论
本研究创新性地构建了基于Mamba状态空间模型的双流深度学习框架,实现了多基因扰动响应的精准预测。该架构通过数据驱动流解析基因表达谱的统计规律,同时借助知识驱动流整合Reactome通路数据库及基因本体(GO)功能注释等生物先验知识,并采用动态跨模态注意力机制实现异构信息的自适应融合,有效克服了单一数据源或知识驱动的建模局限性。核心的Mamba模块通过其线性计算复杂度与选择性记忆机制,在保证运算效率的同时精确建模了基因间长程互作关系,其输入依赖的参数化策略显著增强了复杂生物学特征的提取能力。
尽管双流Mamba框架在多基因扰动预测中展现出优越性能,其在处理超高维组合扰动时仍受限于生物系统的固有复杂性。未来的改进应聚焦于三个战略方向:(1) 系统整合染色质可及性、蛋白质互作网络等多组学特征,构建全景式基因调控图谱;(2) 开发关键生物通路的动态优先级学习机制;(3) 建立不确定性量化框架以增强临床转化场景下的模型可解释性。这些技术突破将推动下一代算法在精准医学中解决更复杂的基因调控难题。
NOTES
*通讯作者。