基于图神经网络的多组学癌症转移预测模型PATH构建与应用

1. 项目缘起：当癌症转移遇上多组学与图神经网络

在肿瘤研究领域，预测癌症是否会转移，以及会转移到哪个器官，一直是个“老大难”问题。传统的临床病理指标，比如肿瘤大小、分级、淋巴结状态，虽然有用，但预测精度常常撞上天花板，很多患者明明指标看着还行，术后几年却出现了远处转移，让人措手不及。这背后的根本原因在于，癌症转移是一个极其复杂的生物学过程，它不仅仅是癌细胞自身的“叛变”，更是癌细胞与肿瘤微环境（包括免疫细胞、基质细胞、血管等）以及全身系统之间一场精密的、多层次的“合谋”。

这就引出了“多组学数据”的价值。单一维度的数据，比如只看基因突变（基因组），或者只看基因表达量（转录组），就像盲人摸象，只能看到一个侧面。而多组学整合分析，则试图把基因组、转录组、表观基因组（如DNA甲基化）、蛋白质组甚至代谢组的数据拼凑起来，形成一个更立体的“肿瘤全景图”。理论上，这张图里应该藏着决定癌细胞是否具备“远征”能力的密码。

然而，问题也随之而来。这些不同来源的数据，维度各异、尺度不同、噪声巨大，且彼此间存在着非线性的、错综复杂的相互作用关系。用传统的统计学方法或者简单的机器学习模型（如逻辑回归、随机森林）来处理这种高维、异构、关系复杂的数据，往往力不从心。模型要么容易过拟合，要么无法有效捕捉不同组学层面之间深层次的生物学关联。

正是在这个背景下，图神经网络进入了我们的视野。GNN的核心能力就是处理关系数据。我们可以很自然地将一个患者的各种分子数据构建成一张图：每个基因、每个蛋白、每个代谢物是一个“节点”，它们之间已知的或预测的生物学相互作用（如蛋白质互作、代谢通路、共表达关系）就是连接节点的“边”。GNN能够在这种图结构上有效地进行信息传递和聚合，从而学习到节点和整个图的表征。这对于挖掘驱动癌症转移的、隐藏在复杂分子网络中的“关键通路”和“协同模块”来说，简直是量身定做的工具。

因此，这个名为“PATH”的模型项目，其核心动机就是：利用图神经网络作为“解码器”，去破译多组学数据构成的复杂生物网络，从而实现对癌症转移风险更精准、更可解释的预测。它不是为了追求模型复杂度而复杂，而是为了解决传统方法在应对此类数据时固有的建模瓶颈。

2. 数据基石：多组学数据的获取、预处理与图构建

任何模型的上限都取决于数据质量。对于PATH模型，数据工程是整个流程中最耗时、也最需要严谨对待的环节。这一步没做好，后面再精巧的模型也是空中楼阁。

2.1 多组学数据来源与挑战

通常，我们会从公共数据库如TCGA、ICGC、GEO等获取癌症队列的多组学数据。一个理想的数据集应包含同一批患者的：

• 基因组：体细胞突变（SNV, Indel）、拷贝数变异（CNV）。
• 转录组：RNA-Seq得到的基因表达量（TPM或FPKM值）。
• 表观基因组：DNA甲基化芯片或测序数据。
• 临床信息：至关重要的标签，包括是否发生转移、转移部位、无转移生存期等。

面临的第一个挑战是数据缺失与异质性。不是每个患者都有全部组学数据，不同平台、不同批次的数据存在巨大的技术偏差。常见的处理流程是：

1. 批次校正：对于表达和甲基化数据，使用ComBat或Harmony等算法校正批次效应。
2. 缺失值处理：对于少量缺失，可采用KNN插补；若某组学数据大量缺失，有时不得不舍弃该组学或该样本，这是一个需要权衡的决策。
3. 特征筛选：高维是灾难的源头。不能把几万个基因全部扔进模型。我们需要进行严格的筛选：
   • 方差过滤：去除在所有样本中表达几乎没有变化的基因。
   • 与表型相关性过滤：通过Cox回归或差异分析，筛选与转移生存显著相关的基因。
   • 通路/网络先验知识：从KEGG、Reactome等数据库聚焦于已知的癌症转移相关通路（如EMT、血管生成、细胞粘附）中的基因。

经过这些步骤，我们可能将每个组学的维度从数万降至数百到一两千个关键特征。

2.2 从表格数据到生物分子网络图

这是PATH模型区别于传统方法的核心一步。我们需要将每个患者的分子数据，从一张特征表格，转化为一张图。

节点定义：最直接的方式是将筛选后的每个基因（或蛋白、代谢物）作为一个节点。一个更精细的做法是，考虑不同组学数据的融合。例如，一个基因节点可以同时拥有多个特征属性：它的表达量（转录组）、它的突变状态（基因组，可编码为0/1或突变类型）、它的启动子区域甲基化水平（表观组）。这样，每个节点就是一个多维向量。

边（关系）定义：这是赋予图以生物学意义的关键。边的构建依赖于先验知识库，主要有以下几种来源：

1. 蛋白质-蛋白质相互作用网络：从STRING、BioGRID等数据库获取。这是最常用、也是最可靠的边来源，代表了蛋白质之间直接的物理或功能关联。
2. 基因共表达网络：可以利用当前数据集或大型参考数据集计算基因间的表达相关性，保留显著的正/负相关关系作为边。但要注意，这可能导致数据泄露，需谨慎。
3. 代谢通路关系：来自KEGG、Reactome。将同一通路内的基因连接起来，或者将具有酶-底物关系的基因连接起来。
4. 调控关系：如转录因子与靶基因的调控关系（来自TRRUST等数据库）。

通常，我们会整合以上多种来源，构建一个相对全面且稳健的生物分子相互作用网络，作为所有患者共享的背景图拓扑结构。

图构建实操：最终，对于第 (i) 个患者，我们构建一张图 (G_i = (V, E, X_i))。

• (V)：节点集合（例如，1500个关键基因），所有患者共享。
• (E)：边集合，由先验知识库定义，所有患者共享相同的拓扑结构。
• (X_i)：节点特征矩阵，维度为 ([节点数, 特征维度])。每一行对应一个节点的多组学特征向量，这个向量是患者特异性的。例如，基因A在患者1中可能是[高表达， 无突变， 低甲基化]，在患者2中则是[低表达， 错义突变， 高甲基化]。

注意：共享拓扑结构是GNN处理这类问题的常见假设，它意味着生物网络的基本架构是保守的，但每个患者在这个网络上的“活动状态”（节点特征）不同，正是这种不同的活动模式决定了其表型（是否转移）。

3. 模型核心：PATH的图神经网络架构设计

PATH模型不是一个单一的GNN，而是一个为多组学癌症数据定制的架构流水线。其核心思想是：先分别从不同组学视角学习特征，再在生物网络图上进行融合与深层推理。

3.1 多组学特征编码器

由于基因组（突变）、转录组（连续值）、甲基化组（比例值）数据分布不同，直接拼接并不合理。PATH模型的第一层通常是针对不同组学数据的独立编码器。

• 对于连续值特征（如表达量、甲基化β值）：可以使用简单的全连接层（Linear）或一维卷积层（Conv1D）进行非线性变换和降维。

  # 伪代码示例：表达量编码器 class ExpressionEncoder(nn.Module): def __init__(self, input_dim, hidden_dim): super().__init__() self.fc1 = nn.Linear(input_dim, hidden_dim) self.bn1 = nn.BatchNorm1d(hidden_dim) self.act = nn.ReLU() self.dropout = nn.Dropout(0.3) def forward(self, x): # x: [batch_size, num_genes] x = self.fc1(x) x = self.bn1(x) x = self.act(x) x = self.dropout(x) return x # 输出统一维度的特征

• 对于离散值特征（如突变类型）：可以先用嵌入层（Embedding Layer）将其转化为稠密向量，再经过全连接层处理。

每个组学编码器将原始高维特征映射到一个统一且较低维度的语义空间（例如，128维）。然后，我们将同一个基因对应的不同组学编码特征拼接起来，形成该基因最终的节点初始特征向量 (h_v^{(0)})。

3.2 图卷积层与信息传递

这是GNN发挥魔力的地方。我们采用图卷积网络（GCN）或图注意力网络（GAT）层作为核心模块。

• GCN层：它对每个节点邻居的特征进行平均聚合，然后通过一个可学习的权重矩阵进行变换。 [ h_v^{(l+1)} = \sigma \left( W^{(l)} \cdot \text{AGGREGATE} \left( { h_u^{(l)}, \forall u \in \mathcal{N}(v) } \right) \right) ] 其中，(\mathcal{N}(v)) 是节点 (v) 的邻居集合，(\sigma) 是非线性激活函数。GCN的优点是简单高效，但它平等对待所有邻居。

• GAT层：它引入了注意力机制，让节点在聚合邻居信息时，可以“有选择地”关注更重要的邻居。 [ \alpha_{vu} = \frac{\exp(\text{LeakyReLU}(a^T[Wh_v || Wh_u]))}{\sum_{k \in \mathcal{N}(v)} \exp(\text{LeakyReLU}(a^T[Wh_v || Wh_k]))} ] [ h_v^{(l+1)} = \sigma \left( \sum_{u \in \mathcal{N}(v)} \alpha_{vu} W h_u^{(l)} \right) ] GAT理论上更强大，能学习生物网络中不同相互作用的重要性权重，但参数量稍大，对数据量要求更高。

在PATH中，通常会堆叠2-3层这样的图卷积层。第一层学习局部邻域特征（如一阶邻居的影响），第二层则可以聚合到更广范围的网络信息（二阶邻居）。经过几层传播后，每个节点的特征 (h_v^{(L)}) 都包含了其自身及其所在局部网络的多组学信息。

3.3 图级读出与预测头

我们的目标是预测患者级别的结局（转移/不转移），因此需要将整张图的信息汇总成一个全局向量，即“图嵌入”。

读出函数至关重要。常见的方法有：

1. 全局平均/最大池化：对所有节点的最终特征取平均或最大值。简单，但可能丢失结构信息。
2. 全局注意力池化：引入一个可学习的查询向量，计算每个节点对图级任务的贡献度（注意力权重），然后加权求和。

   # 伪代码示例：注意力池化 class GlobalAttentionPooling(nn.Module): def __init__(self, in_dim): super().__init__() self.attn_proj = nn.Linear(in_dim, 1) # 计算每个节点的重要性得分 def forward(self, node_features): # node_features: [num_nodes, in_dim] attn_scores = torch.softmax(self.attn_proj(node_features), dim=0) # [num_nodes, 1] graph_embedding = torch.sum(attn_scores * node_features, dim=0) # [in_dim] return graph_embedding

   这种方法允许模型自动识别哪些基因或子网络对预测转移最为关键，增强了模型的可解释性。

得到图嵌入向量后，后面接上传统的全连接层分类器（或Cox比例风险模型用于生存预测），输出最终的预测概率。

一个完整的PATH模型前向传播流程可以概括为：原始多组学表格数据->按组学类型分别编码->拼接成节点特征->输入到共享拓扑的生物网络图->经过2-3层GCN/GAT进行信息传播->通过注意力池化得到图嵌入->全连接层分类/回归。

4. 实战演练：从零搭建PATH模型的代码框架与调优

理论说得再多，不如一行代码。这里我们勾勒一个使用PyTorch Geometric（一个非常流行的图神经网络库）搭建PATH模型简化版的核心框架。

4.1 环境搭建与数据准备

首先，确保安装必要的库。

pip install torch torchvision torchaudio pip install torch-geometric

数据准备部分最为繁琐。假设我们已经有了：

• patient_features.npy: 形状为[num_patients, num_genes, num_omics]的特征矩阵。
• adjacency_matrix.npy: 形状为[num_genes, num_genes]的邻接矩阵（0/1或权重），由先验知识生成。
• labels.npy: 形状为[num_patients]的标签（0/1）。

我们需要将其转换为PyG需要的Data对象列表。

import torch import numpy as np from torch_geometric.data import Data def create_pyg_data(features, adj_matrix, label): """ 为单个患者创建PyG Data对象。 features: [num_genes, num_omics] adj_matrix: [num_genes, num_genes] (scipy sparse matrix format is better) label: scalar """ # 将邻接矩阵转换为边索引（COO格式） edge_index = torch.tensor(np.array(adj_matrix.nonzero()), dtype=torch.long) # 节点特征 x = torch.tensor(features, dtype=torch.float) # 图标签 y = torch.tensor([label], dtype=torch.float) # 回归任务用float，分类用long return Data(x=x, edge_index=edge_index, y=y) # 假设我们已经加载了所有数据 all_patient_data = [] for i in range(num_patients): feat = patient_features[i] # [num_genes, num_omics] data = create_pyg_data(feat, adjacency_matrix, labels[i]) all_patient_data.append(data)

4.2 模型定义

下面是一个结合了多组学编码和GAT的简化PATH模型定义。

import torch.nn as nn import torch.nn.functional as F from torch_geometric.nn import GATConv, global_mean_pool, global_max_pool class OmicsEncoder(nn.Module): """编码单个组学数据""" def __init__(self, input_dim, hidden_dim): super().__init__() self.fc1 = nn.Linear(input_dim, hidden_dim) self.bn1 = nn.BatchNorm1d(hidden_dim) self.dropout = nn.Dropout(0.2) def forward(self, x): # x: [batch_size * num_nodes, input_dim] 或 [num_nodes, input_dim] x = self.fc1(x) x = self.bn1(x) x = F.relu(x) x = self.dropout(x) return x class PATHModel(nn.Module): def __init__(self, omics_dims, # 字典，如 {'expr': 500, 'mut': 50, 'meth': 300} hidden_dim=128, gat_heads=4, # GAT注意力头数 num_gat_layers=2, dropout=0.3): super().__init__() # 1. 多组学编码器 self.omics_encoders = nn.ModuleDict() for omic_name, dim in omics_dims.items(): self.omics_encoders[omic_name] = OmicsEncoder(dim, hidden_dim) # 节点特征融合后的维度 node_in_dim = hidden_dim * len(omics_dims) # 2. 图卷积层（使用GAT） self.gat_convs = nn.ModuleList() self.gat_convs.append(GATConv(node_in_dim, hidden_dim, heads=gat_heads, dropout=dropout)) for _ in range(num_gat_layers - 1): self.gat_convs.append(GATConv(hidden_dim * gat_heads, hidden_dim, heads=gat_heads, dropout=dropout)) # 3. 图读出层（结合平均和最大池化） self.pooling_hidden_dim = hidden_dim * gat_heads # 4. 预测头 self.fc1 = nn.Linear(self.pooling_hidden_dim * 2, 64) # *2因为concat了mean和max self.fc2 = nn.Linear(64, 1) # 二分类输出一个logit self.dropout = nn.Dropout(dropout) def forward(self, data): x, edge_index, batch = data.x, data.edge_index, data.batch # 拆分多组学特征 (假设x是按[expr, mut, meth]顺序拼接的) omics_features = {} start_idx = 0 # 这里需要根据实际特征顺序和维度来拆分，此处为示例 # 实际中可能需要更精细的拆分逻辑 omics_features['expr'] = x[:, :500] omics_features['mut'] = x[:, 500:550] omics_features['meth'] = x[:, 550:] # 分别编码并拼接 encoded_features = [] for omic_name in self.omics_encoders.keys(): enc = self.omics_encoders[omic_name](omics_features[omic_name]) encoded_features.append(enc) x = torch.cat(encoded_features, dim=-1) # [num_nodes_total, node_in_dim] # GAT层前向传播 for i, conv in enumerate(self.gat_convs): x = conv(x, edge_index) if i != len(self.gat_convs) - 1: x = F.elu(x) # GAT常用ELU x = self.dropout(x) # 图级读出 x_mean = global_mean_pool(x, batch) # [batch_size, pooling_hidden_dim] x_max = global_max_pool(x, batch) # [batch_size, pooling_hidden_dim] x_pooled = torch.cat([x_mean, x_max], dim=1) # [batch_size, pooling_hidden_dim*2] # 预测 x = self.fc1(x_pooled) x = F.relu(x) x = self.dropout(x) out = self.fc2(x) # [batch_size, 1] return torch.sigmoid(out) # 输出概率

4.3 训练、验证与关键调优技巧

定义了模型和数据后，训练循环和传统深度学习类似，但有一些需要特别注意的地方。

1. 数据划分的陷阱：千万不能随机打散节点或边来划分训练集和测试集！因为所有患者共享同一张图拓扑。正确的做法是在患者（图）级别进行划分。例如，70%患者的图数据用于训练，30%用于测试。这确保了模型是在未见过的患者上进行评估，模拟真实应用场景。

2. 损失函数与类别不平衡：癌症转移样本通常远少于未转移样本。使用普通的二元交叉熵损失会导致模型偏向于预测多数类。必须使用带权重的BCE损失或Focal Loss。

pos_weight = torch.tensor([(num_negatives / num_positives)]) # 根据你的数据计算 criterion = nn.BCEWithLogitsLoss(pos_weight=pos_weight) # 或者使用Focal Loss

3. 正则化是关键：图神经网络，特别是GAT，在小规模生物数据上很容易过拟合。除了常用的Dropout和权重衰减（L2正则化），图结构上的正则化特别有效：

• Edge Dropout：在每一层图卷积前，随机丢弃一部分边（例如，丢弃率0.2）。这相当于对生物网络进行随机扰动，可以增强模型的鲁棒性，防止其过度依赖少数几条关键的边。

  # 在forward中，可以简单实现 def forward(self, data, training=True): edge_index = data.edge_index if training: # 随机丢弃一部分边 mask = torch.rand(edge_index.size(1)) > self.edge_dropout_rate edge_index = edge_index[:, mask] # ... 后续操作

• 特征Dropout：对节点输入特征也进行随机丢弃。

4. 评估指标：不要只看准确率。对于不平衡的医疗数据，应重点关注：

• AUROC：受试者工作特征曲线下面积，综合衡量分类器性能。
• AUPRC：精确率-召回率曲线下面积，在不平衡数据上比AUROC更敏感。
• F1-Score：精确率和召回率的调和平均。 在验证集上根据AUROC或AUPRC早停是常见的策略。

5. 超参数调优：对PATH模型性能影响最大的超参数通常包括：

• GNN层数：2-3层通常足够，更深可能导致过度平滑（所有节点特征趋于相同）。
• 隐藏层维度：64, 128, 256。
• 注意力头数：4或8。
• 学习率与优化器：AdamW优化器，学习率1e-3到1e-4，配合余弦退火调度器。
• Dropout率：0.3到0.6。
• 边丢弃率：0.1到0.3。

建议使用贝叶斯优化或超参数搜索库（如Optuna）进行系统调优。

5. 结果解读与模型可解释性：不只是黑箱预测

对于一个旨在辅助临床决策的模型来说，仅仅给出“高转移风险”的预测是不够的。医生和生物学家更关心的是：为什么这个患者风险高？是哪些基因和通路在起作用？因此，模型的可解释性与预测性能同等重要。

5.1 节点与边的重要性分析

得益于GNN的架构，我们可以采用一些方法来追溯模型的决策依据。

• 节点重要性（基因重要性）：在使用了注意力池化的模型中，最终的图嵌入是各个节点特征的加权和。这个注意力权重(\alpha_v) 直接反映了节点 (v) 对最终预测的贡献度。我们可以对权重进行排序，找出对预测转移风险最重要的前N个基因。这些基因很可能是驱动转移的关键因子。

• 边重要性（相互作用重要性）：在GAT中，每一层都计算了节点对之间的注意力系数 (\alpha_{vu})。我们可以分析这些系数，找出那些被模型赋予高权重的边（相互作用）。例如，如果模型高度关注基因A和基因B之间的边，这可能暗示着这两个基因的协同失调在转移过程中扮演重要角色。我们可以将这些重要的边映射回KEGG等通路数据库，看它们是否富集在特定的生物学通路中。

5.2 基于梯度的可解释性方法

对于不使用注意力的GCN，或者为了进行补充验证，我们可以使用基于梯度的方法。

• Saliency Maps：计算模型输出（预测概率）相对于输入节点特征的梯度。梯度绝对值大的特征，意味着该特征的微小变化会对预测结果产生较大影响，从而表明该特征重要。

  model.eval() data.requires_grad_(True) # 需要计算输入的梯度 output = model(data) output.backward() # 反向传播 node_saliency = data.x.grad.abs().sum(dim=1) # 对每个节点的多组学特征梯度取绝对值和

  这可以告诉我们，是哪些基因的哪些组学特征（如某个基因的高表达或特定突变）对预测贡献最大。

• GNNExplainer：这是一个专门为GNN设计的模型无关解释工具。它可以为单个预测样本生成一个解释子图，即一个小的、连通的节点和边集合，这个子图是驱动模型做出该预测的最重要因素。

  from torch_geometric.nn import GNNExplainer explainer = GNNExplainer(model, epochs=200) node_feat_mask, edge_mask = explainer.explain_graph(data.x, data.edge_index)

  edge_mask值高的边构成了解释子图。可视化这个子图，并与已知的癌症转移通路进行比对，能提供非常直观的生物学见解。

5.3 生物学验证与发现

将模型识别出的重要基因和相互作用列表，进行下游的生物学富集分析（如GO、KEGG富集分析），是验证模型是否学到真实生物学信号的关键步骤。如果富集到的通路恰好是已知的转移相关通路（如TGF-β信号通路、Wnt信号通路、细胞粘附分子通路等），那么模型的预测就具备了生物学合理性，增加了可信度。

更进一步，这些分析可能产生新的科学假设。例如，模型可能高亮了一个尚未被充分研究与转移相关的基因或基因互作对，这可以成为后续湿实验验证的起点。这才是AI for Science的真正价值所在——不仅预测，更帮助发现。

6. 局限、挑战与未来方向

尽管PATH模型框架展现了巨大潜力，但在实际研究和应用中，我们必须清醒地认识到其面临的挑战和当前局限。

数据层面的挑战：

1. 多组学数据对齐与质量：不同组学数据来自不同平台，批次效应显著。即便经过校正，残留的噪声仍可能淹没微弱的真实信号。如何更有效地进行跨组学数据整合与清洗，是一个持续的研究课题。
2. 样本量有限：具有完整多组学数据和长期随访信息的癌症队列样本量通常只有数百到数千。在如此高的特征维度（数千个基因节点）下，深度学习模型极易过拟合。除了强力的正则化，利用迁移学习（在大型单组学数据上预训练编码器）或采用更简单的GNN架构是务实的选择。
3. 静态图的假设：PATH模型假设生物网络拓扑是静态且共享的。但实际上，不同亚型、不同阶段的肿瘤，其内部的分子互作网络可能是动态变化的。如何构建患者特异性的、或条件依赖性的图结构，是一个前沿方向。

模型层面的挑战：

1. 可解释性与性能的权衡：更复杂的注意力机制可能提升性能，但有时会降低可解释的清晰度。如何设计既强大又易于理解的GNN架构是关键。
2. 异质图的处理：当前模型将不同组学特征拼接为节点属性，这是一种同质图处理方式。更自然的建模方式是构建异质图，其中节点类型可以是基因、甲基化位点、蛋白等，边类型可以是调控、互作、共表达等。异质图神经网络能更细致地建模不同实体间的关系，但模型复杂度和数据需求也更高。
3. 整合空间与时间信息：目前的模型缺乏空间信息（如肿瘤内不同区域的空间转录组）和时间维度（治疗前后或转移前后的纵向样本）。未来的模型需要向时空图神经网络发展，以捕捉癌症演进的全景。

应用层面的考量：

1. 临床转化路径长：从算法性能到临床可用工具，需要经过严格的独立队列验证、前瞻性临床试验，并满足医疗器械监管要求。模型的稳定性、鲁棒性和公平性（在不同人种、性别间的表现）需要被充分评估。
2. 计算成本：训练深度GNN需要GPU资源，且推理速度需满足临床实时或准实时需求。模型压缩和轻量化是部署前必须考虑的问题。

在我个人的多次尝试中，一个深刻的体会是：不要一味追求最复杂的模型。对于很多中型规模的队列，一个2层GCN配合精心设计的特征筛选和强正则化，其表现往往比一个深层的、参数众多的GAT更稳定、更可解释。生物数据的信噪比低，有时“简单而稳健”比“复杂而脆弱”更有价值。另一个实用的技巧是，在训练初期，可以冻结GNN层的参数，只训练编码器和分类器，待损失平稳后再解冻全部参数进行微调，这有助于稳定训练过程。

PATH模型代表了一种强大的范式转变——从将生物数据视为孤立的特征点，到将其视为相互关联的网络系统进行分析。虽然前路仍有诸多挑战，但它无疑为我们更深入地理解癌症转移这一复杂过程，并最终实现更精准的预后预测和靶向治疗，打开了一扇充满希望的新窗口。