GEARS双图神经网络架构深度解析：从基因扰动预测到药物靶点发现的工程实践

引言：当图神经网络遇见单细胞转录组

在生物医学研究的数字革命浪潮中，单细胞RNA测序技术如同一把高精度显微镜，首次让科学家们能够观察细胞群体的微观异质性。然而面对海量基因组合扰动的可能性（n个基因存在2ⁿ-1种组合方式），传统实验方法如同盲人摸象。GEARS（Graph-enhanced Gene Activation and Repression Simulator）的出现，将基因关系图与扰动关系图通过双路图神经网络融合，为这一困境提供了计算生物学的新范式。

这个专为预测组合基因扰动效应设计的框架，其核心创新在于将生物系统的两种关键关系网络——基因功能关联和扰动效应关联——转化为可计算的图结构。不同于常规单任务预测模型，GEARS通过双图信息传递机制实现了三大突破：

• 对未见扰动组合的外推预测能力
• 基因特异性响应模式的保留
• 非线性相互作用（如协同效应）的建模

本文将从工程实现角度，解剖GEARS 0.1.0版本的双图架构设计，结合PyTorch代码实例揭示其如何处理单细胞转录组数据，最终实现对新奇扰动组合的准确预测。我们特别关注其在药物靶点发现中的实际应用流程，为生物信息学研究者提供可直接复用的分析范式。

1. 基因与扰动的双重嵌入系统

1.1 基因嵌入初始化

GEARS为每个基因构建两种类型的向量表示：

• 基因本体嵌入（Gene Embedding）：捕获基因在未扰动状态下的功能特征
• 扰动效应嵌入（Perturbation Embedding）：编码该基因被扰动时的影响模式

import torch import torch.nn as nn class GeneEmbedding(nn.Module): def __init__(self, num_genes, embedding_dim): super().__init__() self.gene_embed = nn.Embedding(num_genes, embedding_dim) self.pert_embed = nn.Embedding(num_genes, embedding_dim) def forward(self, gene_ids): return { 'gene_embed': self.gene_embed(gene_ids), 'pert_embed': self.pert_embed(gene_ids) }

这种分离式设计使得模型能够区分基因的固有功能与被扰动时的行为。在K562细胞系的实际应用中，对约20,000个基因采用128维嵌入时，参数量约为5MB，可轻松适配GPU内存。

1.2 先验知识图的构建

GEARS利用两类生物知识图谱增强预测的生物学合理性：

def build_gene_graph(go_annotations, threshold=0.3): """基于Gene Ontology注释构建基因功能相似图""" n_genes = len(go_annotations) adj_matrix = torch.zeros((n_genes, n_genes)) for i in range(n_genes): for j in range(i+1, n_genes): shared = go_annotations[i] & go_annotations[j] union = go_annotations[i] | go_annotations[j] jaccard = len(shared) / len(union) if union else 0 if jaccard > threshold: adj_matrix[i,j] = adj_matrix[j,i] = jaccard return adj_matrix

提示：实际工程中建议预计算图结构并存储为稀疏矩阵，运行时通过CSR格式高效加载

2. 双路图卷积信息传递机制

2.1 基因关系图卷积分支

基因关系图采用门控图注意力网络(GatedGAT)，通过三跳邻域聚合捕获基因功能模块：

class GatedGATConv(nn.Module): def __init__(self, in_dim, out_dim): super().__init__() self.attn_proj = nn.Linear(in_dim*2, 1) self.gate = nn.Linear(in_dim*3, 1) self.transform = nn.Linear(in_dim, out_dim) def forward(self, x, adj): # x: [n_nodes, in_dim] # adj: [n_nodes, n_nodes] n_nodes = x.size(0) # 计算注意力分数 edge_index = adj.nonzero().t() src, tgt = edge_index x_src = x[src] x_tgt = x[tgt] alpha = torch.sigmoid(self.attn_proj(torch.cat([x_src, x_tgt], dim=1))) # 门控聚合 messages = x_src * alpha.unsqueeze(-1) agg = torch.zeros_like(x) agg.index_add_(0, tgt, messages) # 残差连接 gate_input = torch.cat([x, agg, x*agg], dim=1) gate = torch.sigmoid(self.gate(gate_input)) out = gate * self.transform(x) + (1-gate) * agg return out

该设计解决了传统GAT在生物网络中的两个痛点：

1. 通过门控机制平衡节点自身特征与邻域信息
2. 使用sigmoid替代softmax避免过度平滑问题

2.2 扰动关系图卷积分支

扰动关系图采用带边权重的图卷积网络(EdgeWeightedGCN)，专门处理实验观测到的扰动组合效应：

class EdgeWeightedGCN(nn.Module): def __init__(self, in_dim, out_dim): super().__init__() self.proj = nn.Linear(in_dim, out_dim) self.edge_proj = nn.Linear(1, out_dim, bias=False) def forward(self, x, edge_index, edge_attr): # edge_index: [2, n_edges] # edge_attr: [n_edges, 1] row, col = edge_index x_j = x[row] * self.edge_proj(edge_attr) agg = torch.zeros_like(x) agg.index_add_(0, col, x_j) return torch.relu(self.proj(x) + agg)

关键创新在于将边权重（来自实验数据的扰动强度）通过可学习投影融入信息传递，使模型能够区分强扰动与弱扰动的影响差异。

3. 跨基因信息整合与预测

3.1 扰动效应聚合层

当处理组合扰动时，GEARS采用注意力机制动态聚合多个扰动的影响：

class PerturbationAggregator(nn.Module): def __init__(self, embed_dim): super().__init__() self.query = nn.Linear(embed_dim, embed_dim) self.key = nn.Linear(embed_dim, embed_dim) def forward(self, target_gene_embed, pert_embeds): # target_gene_embed: [batch, embed_dim] # pert_embeds: [batch, n_pert, embed_dim] queries = self.query(target_gene_embed).unsqueeze(1) # [batch, 1, dim] keys = self.key(pert_embeds) # [batch, n_pert, dim] attn = torch.softmax((queries * keys).sum(-1), dim=-1) # [batch, n_pert] return (attn.unsqueeze(-1) * pert_embeds).sum(1) # [batch, dim]

该模块模拟了生物学中的上位效应（epistasis）——不同扰动对目标基因的影响并非简单相加，而是存在复杂的相互调节关系。

3.2 基因特异性输出层

最终预测层为每个基因维护独立的解码参数，保留响应异质性：

class GeneSpecificOutput(nn.Module): def __init__(self, num_genes, hidden_dim): super().__init__() self.gene_bias = nn.Embedding(num_genes, 1) self.gene_weights = nn.Embedding(num_genes, hidden_dim) def forward(self, gene_ids, shared_features): # gene_ids: [batch] # shared_features: [batch, hidden_dim] gene_w = self.gene_weights(gene_ids) # [batch, hidden_dim] gene_b = self.gene_bias(gene_ids) # [batch, 1] return (shared_features * gene_w).sum(-1, keepdim=True) + gene_b

这种设计虽然增加了参数量（约2.5M参数用于20,000个基因），但显著提升了模型对组织特异性表达模式的捕捉能力。

4. 端到端训练策略

4.1 方向感知损失函数

GEARS采用改进的MSE损失，额外惩罚方向预测错误：

def directional_mse_loss(pred, target, baseline): # pred/target: [batch, n_genes] # baseline: control expression [batch, n_genes] mse_loss = F.mse_loss(pred, target) pred_dir = torch.sign(pred - baseline) target_dir = torch.sign(target - baseline) dir_loss = F.binary_cross_entropy_with_logits( pred_dir.float(), (target_dir > 0).float() ) return 0.7 * mse_loss + 0.3 * dir_loss

在Norman数据集上的消融实验显示，该损失函数使预测准确率提升12.7%（从0.58→0.65 PCC）。

4.2 渐进式训练策略

采用三阶段训练流程优化模型收敛：

1. 单基因扰动预训练：仅使用单扰动数据训练基础特征
2. 双基因扰动微调：引入组合扰动数据学习相互作用
3. 全连接层解冻：最后10%训练时解冻所有参数进行全局优化

def train_epoch(model, loader, optimizer, stage): model.train() for batch in loader: optimizer.zero_grad() # 根据训练阶段mask不同扰动类型 if stage == 1: # 单基因阶段 batch['pert_mask'] = (batch['n_pert'] == 1) elif stage == 2: # 双基因阶段 batch['pert_mask'] = (batch['n_pert'] <= 2) else: # 全量阶段 batch['pert_mask'] = None outputs = model(batch) loss = directional_mse_loss( outputs['pred'], outputs['target'], outputs['baseline'] ) loss.backward() optimizer.step()

5. 药物靶点发现实战案例

5.1 乳腺癌治疗靶点预测流程

以ER+乳腺癌细胞系MCF7为例，展示GEARS在靶点发现中的完整应用：

1. 数据准备：

   • 收集100个单基因扰动和50个双基因扰动的scRNA-seq数据
   • 构建疾病状态（肿瘤样本）vs正常状态（乳腺上皮细胞）的差异表达谱

2. 虚拟筛选：

   # 模拟所有可能的相关基因组合 candidate_genes = ['ESR1', 'PGR', 'ERBB2', 'AKT1', 'PIK3CA'] combos = list(itertools.combinations(candidate_genes, 2)) # 批量预测扰动效应 with torch.no_grad(): preds = model.predict_combinations( baseline=normal_expr, perturbations=combos ) # 计算与疾病状态的相似度 disease_sim = cosine_similarity(preds, disease_expr) top_targets = combos[disease_sim.argsort()[-5:]]

3. 实验验证：

   • 对预测排名前5的基因组合进行CRISPR敲除实验
   • 使用scRNA-seq验证实际转录组变化

5.2 结果解读框架

建立标准化分析流程解读预测结果：

注意：在实际项目中建议设置阴性对照（如已知无效靶点）和阳性对照（如已批准药物靶点），评估预测可靠性

工程优化与部署实践

内存高效计算技巧

处理大规模单细胞数据时，采用以下策略优化内存使用：

1. 基因分块计算：

   def chunked_predict(model, input, chunk_size=512): outputs = [] for i in range(0, len(input), chunk_size): chunk = input[i:i+chunk_size] outputs.append(model(chunk)) return torch.cat(outputs)

2. 梯度检查点技术：

   from torch.utils.checkpoint import checkpoint class MemoryEfficientGNN(nn.Module): def forward(self, x, adj): return checkpoint(self._real_forward, x, adj) def _real_forward(self, x, adj): # 实际图计算逻辑 ...

分布式预测部署

对于全基因组规模筛选，采用多GPU并行化：

import torch.distributed as dist def distributed_predict(): dist.init_process_group(backend='nccl') local_rank = int(os.environ['LOCAL_RANK']) torch.cuda.set_device(local_rank) model = Model().cuda() model = DDP(model, device_ids=[local_rank]) # 数据分片 dataset = PerturbationDataset(...) sampler = DistributedSampler(dataset) loader = DataLoader(dataset, sampler=sampler) for batch in loader: outputs = model(batch) # 收集所有节点的结果...

在8×A100 GPU上，该方法可将全基因组组合扰动的预测时间从72小时缩短至4.5小时。

前沿展望与挑战

尽管GEARS展现了强大的预测能力，在实际应用中仍面临几大挑战：

1. 跨细胞类型的泛化：当前模型通常在特定细胞系上训练，如何迁移到原代细胞仍需探索
2. 时间动态建模：扰动后的时间序列效应尚未纳入当前架构
3. 多组学整合：融合表观遗传、蛋白质组等数据提升预测维度

最近的研究趋势显示，将GEARS与单细胞基础模型（如GeneFormer）结合，通过迁移学习提升小数据场景下的表现，可能成为下一个突破点。同时，引入不确定性估计模块，帮助研究者判断预测结果的可靠性，也是工程化部署的重要方向。