去年冬天接手一个注释收尾的活儿:一批从差异表达里筛出来的候选蛋白,序列比对跑完,接近四成在公共库里只能匹配到"功能未知"或者"某类结构域超家族成员"。这种结果交上去基本等于没做,导师看一眼就会问"那它到底干什么"。后来我把思路从"看它像谁"换成"看它跟谁一起干活"——把PPI 网络(蛋白-蛋白相互作用网络)拉出来构造特征,再挂一个机器学习模型去预测它可能落在哪些功能条目上。这套蛋白质功能预测的玩法不算新,但要真做出能上报表、经得起追问的结果,中间要处理的问题比想象的多得多。这篇就把我从数据清洗、特征构造、模型选型到评估踩坑的完整链路写下来,适合已经会用 Python 跑基础建模、想把这套思路落到自己数据上的人看,也适合刚进这个方向、被各种数据库和指标绕晕的同学当作路线图。
1. 从"序列搜不到同源"到"从网络里找线索"
1.1 序列相似性搜索的天花板在哪
做功能注释的第一反应永远是 BLAST 或者更灵敏的比对工具,这套逻辑的底层假设是"序列相似则功能相似"。它在很多场景下确实好用,尤其是酶家族和保守结构域,但天花板也很明显。第一类是远缘同源,序列分歧到 20% 以下同一性时,比对已经抓不住信号,但结构折叠和功能位点可能还在。第二类是多结构域蛋白,一段序列里拼了三四个结构域,整体比对的最佳命中往往只对应其中一个结构域的功能,剩下的部分被浪费掉。第三类是非酶类调控蛋白,它们的功能不在于催化什么反应,而在于"和谁结合、在什么复合体里出现",这类功能信息在序列里几乎不编码。
我印象最深的是一个只有两百多个氨基酸的候选蛋白,比对结果清一色是"未表征"。但把它丢进互作网络里,第一层邻居里有四个都是某个转录调控复合体的已知成员,第二层邻居大量富集在某个信号通路上。这种"位置信息"是序列比对完全给不出来的,也是我认为这套方法真正有价值的地方。
提示:不要一上来就否定序列比对,正确做法是把它当作其中一个特征来源,而不是唯一判据。我的习惯是先跑序列相似性,把"高置信度同源"的部分直接继承注释,剩下的"灰色地带"才交给网络模型,这样能把算力和人工核对成本压下来一大截。
1.2 "关联有罪"假设能走多远
PPI 网络做功能预测的核心假设非常朴素:互作的蛋白倾向于参与同一或相关的生物学过程。行业内一般叫它 guilt by association,直译过来就是"关联有罪"。这个假设在复合体成员、同一通路的上下游激酶-底物这类关系上非常扎实,因为物理互作本身就意味着空间和时间上的共定位。
但它有明确的失效边界,这点很多教程不讲。第一种失效是瞬时互作与调控型互作,比如一个泛素连接酶和它的底物,两者功能方向可能完全相反,一个是降解执行者,一个是被降解的靶子。第二种是跨区室的偶发接触,高通量实验里假阳性互作往往就长这样。第三种是枢纽蛋白(hub),一个连接度上千的蛋白几乎和所有通路都沾边,它的邻居注释分布极其分散,用它来传播信息反而会污染周围的低度节点。
所以我在实际项目里做了两件事:一是只用高置信度子集构图,把综合数据库里分数低于某个阈值的边砍掉;二是在特征里显式加入邻居的度分布统计,让模型自己学会判断"这个邻居的注释可信度有多高"。这两步做下来,模型在低度节点上的预测质量提升有限,但在中枢节点上的误报明显减少。
1.3 两条技术路线的分野
把网络信息用进预测,业内其实有两条差异很大的路线,混着用很容易出问题。路线 A 是传播式:不做训练,直接把已知注释当作种子,在网络上游走或扩散,得到每个节点对每个功能条目的得分。带重启的随机游走、k 步邻居加权投票都属于这一类,优点是无需标签、可解释、跑得快,缺点是无法融合序列和表达这类非网络特征。路线 B 是监督学习式:把每个"蛋白-功能条目"对当作一个样本,构造特征向量,训练分类器判断这个条目是否属于这个蛋白。它能随意融合多源特征,性能上限更高,但对负样本构造和划分方式极度敏感。
我自己最终采用的是 B 为主、A 为辅的混合方案:先用传播算法生成一批"传播得分"特征,再连同拓扑特征、序列特征一起喂给梯度提升树。这么设计的原因是传播得分本质上把网络的多跳结构压缩成了低维数值,比让树模型自己去从原始邻接矩阵里找规律要高效得多,而且它天然携带了"已知注释的分布先验",在小样本类别上特别管用。
2. 把一张互作图变成可训练矩阵:数据准备的四道坎
2.1 数据源选型:覆盖度、置信度与物种一致性
公开的互作数据来源大致分三类:综合型数据库(聚合多个来源并给出统一打分)、策展型数据库(人工从文献里抠出来的实验验证互作)、高通量实验沉淀库(双杂交、亲和纯化质谱的原始结果)。综合型的覆盖最广,人源能到一万九千个蛋白左右,高置信度子集的边数在几十万量级;策展型的质量最高但规模小一个数量级;原始高通量库噪声大,非专业人士不建议直接拿来构图。
我一般的取舍是:主图用综合库的高置信度子集,再用策展库做交叉验证。具体做法是先用综合库构图跑一遍模型,然后用策展库里那些"只在策展库出现"的边做外部测试——如果模型在这批边对应的功能预测上也表现不错,说明它不是只拟合了某一个数据源的偏好。选阈值这件事没有标准答案,我的经验是从"能保留住你关心的那批已知复合体"这个角度反推,比如某个你确定存在的复合体,它的成员之间至少要能连成一片,如果阈值高到把它打散了,那就是砍过头了。
| 数据类型 | 典型规模(人源) | 主要优点 | 主要风险 |
|---|---|---|---|
| 综合型互作库高置信子集 | 约 1.5 万–1.9 万节点,数十万条边 | 覆盖广,格式统一,自带打分 | 打分标准不透明,不同版本差异大 |
| 人工策展实验互作 | 数千到一万余节点 | 可信度高,多带实验方法标注 | 规模小,偏向研究热点 |
| 高通量原始结果 | 规模最大 | 可能包含新发现 | 假阳性率高,需要专门清洗 |
| 共表达/共现推断 | 覆盖几乎所有基因 | 能补上孤立节点 | 只反映相关性,不能当物理互作 |
注意:不同数据库对"互作"的定义完全不同,物理结合、遗传互作、通路共现、文本共现经常被混在一个字段里。构图之前务必确认你用的是哪一种,否则"关联有罪"的前提就直接塌了。
2.2 标识符统一:一个被严重低估的坑
这一步看着琐碎,但它是整个流程里最容易埋雷的地方。同一个蛋白可能以基因符号、旧版基因符号、UniProt 登录号、Ensembl 基因 ID、RefSeq 编号五种形式出现在不同表里,而基因符号是会变的——有些符号因为历史原因被重新分配给了别的基因。我吃过一次亏:两个不同蛋白因为共用一个过时符号,在合并时被当成了同一个节点,导致整张图的局部结构错乱,模型在那一片区域的预测全部偏移。
我的固定做法是统一到 UniProt 登录号,并且在项目开始时把映射表连同版本号一起冻结下来,存一份到项目目录里。所有后续处理都读这份冻结的映射,不再联网更新。原因很简单:你在三月和六月各跑一次流程,如果映射表变了,两次结果就不可比,而做科研最怕的就是"同样的代码跑出不同的数"。对于映射不上的节点,直接丢弃而不是猜,宁可少几千个节点,也不要引入身份错误。
2.3 清洗细节:自环、重复边与最大连通片
拿到边表之后的标准清洗步骤我列一下,每步都说明为什么。去自环:自己和自己互作的记录基本是数据处理错误或者自调控的特殊情形,留在图里会干扰聚类系数这类局部指标。合并重复边:同一对蛋白被多个实验或数据库重复记录,这里有个选择——是留一条还是累加权重。我倾向于保留权重(重复次数或打分均值),因为重复证据本身是可信度的信号,但要注意后续做嵌入时要对权重做归一化,否则高度重复的边会主导随机游走。
对称化处理:多数互作数据库的方向性不可靠,统一按无向图处理更稳妥,除非你专门做的是激酶-底物这类有向关系。取最大连通片:人源高置信度网络的最大连通片通常能覆盖 95% 以上的节点,剩下那些零散的小组件和孤立点,说实话对训练没有贡献,但不要删掉它们——它们恰恰是你实际要预测的新蛋白的典型形态。正确做法是训练时排除、推理时保留,并且单独报告模型在低度节点上的性能。
import pandas as pd import networkx as nx # edges: columns = [src_uniprot, dst_uniprot, score, n_evidence] edges = pd.read_csv("interactions_filtered.tsv", sep="\t") # 1) 去掉自环 edges = edges[edges["src_uniprot"] != edges["dst_uniprot"]] # 2) 同一对蛋白合并,方向无关(先排序再聚合) edges["a"] = edges[["src_uniprot", "dst_uniprot"]].min(axis=1) edges["b"] = edges[["src_uniprot", "dst_uniprot"]].max(axis=1) edges = (edges.groupby(["a", "b"]) .agg(score=("score", "mean"), n_ev=("n_evidence", "sum")) .reset_index()) # 3) 建无向图,边权保留 G = nx.from_pandas_edgelist( edges, "a", "b", edge_attr=["score", "n_ev"], create_using=nx.Graph() ) # 4) 标记最大连通片,但不要物理删除孤立点 components = sorted(nx.connected_components(G), key=len, reverse=True) lcc_nodes = components[0] print(f"节点 {G.number_of_nodes()}, 边 {G.number_of_edges()}, " f"最大连通片占比 {len(lcc_nodes)/G.number_of_nodes():.3f}")2.4 标签空间:GO 三分支与稀疏性现实
功能预测的标签用的是基因本体(GO),分三个分支:生物学过程(BP)、细胞组分(CC)、分子功能(MF)。三个分支的预测难度差别很大,CC 最容易,因为定位信息在互作网络里体现得最直接——一个蛋白和核孔复合体成员绑在一起,那它大概率在核膜附近。BP 最难,因为同一个蛋白可能参与十几个过程,且过程之间有层级和时序关系。
这里有个非常现实的问题:GO 条目总数在四万以上,但绝大多数条目的已知注释蛋白数不到十个。这类条目根本没法学,正样本太少,模型只能记住训练集里那几个蛋白。我的做法是按正样本数量卡阈值,每个分支单独定,正样本少于 20–30 个的条目直接排除在训练目标之外。这么砍完,实际可训练的条目大概只剩几百到一两千个,但模型在这些条目上的表现才是有意义的。
另外要注意 GO 的层级结构:一个蛋白被注释到某个子节点,按 true path rule 它必然也该被注释到所有祖先节点。有些数据源直接给了带祖先的完整注释(也叫"传播后注释"),有些只给最具体的叶子。用哪种必须明确,因为两种标签空间下的性能数字完全不可比。我习惯用传播后的,然后在预测阶段再做一次向上传播的一致性后处理,这两步配合能把整体的 Fmax 抬上去几个百分点。
3. 特征工程:手工拓扑与网络嵌入怎么配
3.1 手工拓扑特征清单及其物理含义
拓扑特征的好处是快、可解释、对数据规模不敏感。我常用的那批,每一个都能对应到某种生物学直觉,这也是我在写报告时能解释清楚的地方。
| 特征 | 计算方式 | 生物学含义 | 注意事项 |
|---|---|---|---|
| 度 / 加权度 | 直接统计邻居数与边权之和 | 参与互作的广度 | 与研究热度强相关,慎单独使用 |
| 聚类系数 | 邻居之间实际连边比例 | 是否处在紧密复合体中 | 高度节点普遍偏低 |
| k-core | 迭代剥离度数小于 k 的节点得到的层级 | 网络核心程度 | 对阈值敏感 |
| PageRank | 幂迭代(阻尼 0.85,迭代 100) | 全局重要性 | 与度高度共线 |
| 介数中心性 | 采样近似(k=500 源点) | 是否处在通路交叉口 | 精确计算 O(nm),必须近似 |
| 邻居平均度 | 一阶邻居度的均值 | 邻居是枢纽还是边缘 | 对区分枢纽蛋白有用 |
| 二步邻居数 | 邻接矩阵平方的行和 | 局部模块规模 | 稀疏矩阵下用 A@A 求行和 |
| 三角计数 | 局部三角形数量 | 模块紧密度 | 与聚类系数共线 |
算这些特征时有个性能问题要注意:两万节点的图,介数中心性的精确算法在单机上跑要几个小时。用k=500的采样近似能把时间压到分钟级,误差对树模型来说可以接受。PageRank 和 k-core 都是秒级,可以放心算。一个实操建议:把所有这些特征算完存成一个 parquet,后面所有模型实验都读这份缓存,别每次重算,否则调参调到一半时间全花在算图上。
3.2 传播得分:把多跳结构压成低维向量
这是我个人认为性价比最高的一类特征。核心思路是:把已知注释的蛋白当作种子,让"注释信息"沿着网络扩散,扩散完之后每个节点在每个 GO 条目上都有一个得分,这个得分本身就是极好的特征。实现上我用带重启的随机游走,公式是r = α P r + (1-α) e,其中 P 是行归一化的转移矩阵,e 是种子向量,α 取 0.85。
直接的矩阵求逆(I - αP)^{-1}在稀疏矩阵上其实也能算,但两万乘两万虽然稀疏、求逆之后会变稠密,内存吃不住。所以实际写代码时用固定点迭代,迭代 20–30 次基本收敛,误差在 1e-6 以下就不会影响下游了。更省事的方式是把它写成矩阵形式:把所有 GO 条目的种子拼成一个稀疏矩阵Y0(行是蛋白、列是条目),然后做Y_{t+1} = α P Y_t + (1-α) Y0,迭代 5–10 步就够。因为 GO 注释本身有长尾,迭代太多次反而会把信号抹平到全网。
import numpy as np import scipy.sparse as sp def rwr_propagate(A, Y0, alpha=0.85, steps=10): """ A: (n, n) 稀疏邻接矩阵(对称,未归一化) Y0: (n, m) 稀疏种子矩阵,m 为 GO 条目数 返回: (n, m) 稠密传播得分 """ deg = np.asarray(A.sum(axis=1)).ravel() deg[deg == 0] = 1.0 D_inv = sp.diags(1.0 / deg) P = D_inv @ A # 行归一化转移矩阵 Y = Y0.astype(np.float32).toarray() Y0d = Y.copy() for _ in range(steps): Y = alpha * (P @ Y) + (1 - alpha) * Y0d return Y有个细节值得说:种子矩阵里,用实验证据支撑的注释和电子推断的注释要分开加权。我给实验证据的权重是 1.0,电子推断的给 0.3–0.5。这样做是因为电子推断注释本身就是从序列相似性推出来的,如果把它和实验注释同等对待,等于把"序列相似"这个信号又偷偷塞回网络特征里,会削弱网络信息带来的独立增益。
3.3 网络嵌入:随机游走还是图神经网络
网络嵌入是另一条路,把节点映射到低维稠密向量。随机游走类方法(如把节点序列当作"句子"训练 skip-gram)实现简单、无需标签、在大图上很快,两万节点跑一遍也就十几分钟。它的问题是只编码结构,不编码节点的任何其他属性,而且对超参数(游走长度、窗口大小、每节点游走次数)比较敏感,不同随机种子跑出来的向量在下游任务上性能波动能到两三个百分点。
图神经网络则是端到端的,能把节点特征和邻接结构一起学。它的优势在于能针对具体任务优化表示,劣势是训练成本高、对超参敏感、而且有个著名的问题:层数一多就过平滑,所有节点的表示趋同,性能反而下降。我实测下来两到三层是最优区间,超过四层基本就没法看了。另外 GNN 在小图上(几千节点以下)相比手工特征 + 树模型没有明显优势,它的价值主要在中大规模图上体现。
我的实际配置是:小规模实验用随机游走 + 树模型快速验证特征方向,正式实验才上 GNN,两者并行跑,最后做集成。集成的时候不是简单平均,而是按节点度分箱,低度节点更信任传播类特征(因为它们邻居少,GNN 学不到东西),高度节点更信任 GNN 和拓扑特征。
3.4 序列与表达特征的插入位置
网络特征决定了下限,但真正把性能从"能用"推到"好用"的往往是序列特征。这里不必自己训模型,直接用现成的蛋白语言模型提嵌入就行——把每条序列过一遍,取最后一层的平均池化向量,通常是 1280 维或 1024 维。这个向量的好处是已经把远缘同源的结构信息编码进去了,比 k-mer 频率和 PSSM 强太多。
维度太高会拖慢训练,我一般用 PCA 降到 128 维再拼。另外两类值得加的特征:亚细胞定位预测的得分向量(对 CC 分支特别有用),以及组织表达谱的共表达相关性(对 BP 分支有增量价值,尤其是那些只在特定组织里起作用的蛋白)。融合方式上,我试过早期拼接和晚期融合两种,早期拼接更简单、在小样本条目上更稳,晚期融合在数据量大时略好一点,差距不大,不值得为此增加工程复杂度。
4. 模型训练:从树模型基线到图网络
4.1 先用梯度提升树立一个能对话的基线
我强烈建议任何网络相关项目都先建一个不使用图结构的基线,哪怕它性能很差。原因有二:一是它给你一个明确的"网络到底带来了多少增益"的参照;二是树模型训练快、可解释、调参直观,能让你在半天内摸清特征里哪些有用。
具体配置上,用 LightGBM 的话我一般这么起步:叶子数 63,学习率 0.05,树 800–1500 棵配早停(早停轮数 50),列采样 0.8,行采样 0.8。正负比不平衡特别严重时,用scale_pos_weight按实际比例设,但不要设得太激进——我试过直接按 1:100 设,结果召回上去了精确率掉得厉害,Fmax 反而更低。更稳的做法是把权重设在实际比例的一半左右,剩下靠阈值调整解决。
样本构造上,每个样本是"某个蛋白是否具有某个 GO 条目"。如果只对正样本训练,就变成了一分类问题。我的做法是给每个正样本配 3–5 个负样本,负样本从"同分支下其他条目"里随机抽——注意不能从"该蛋白已有注释的条目"里抽,这个约束必须严格检查,否则就是标签泄漏。
4.2 图神经网络的两个实操要点
如果用 GNN,第一个要点是邻居采样。全图训练在两万节点上勉强可行,但一旦加上序列特征的高维输入,显存立刻吃紧。用两层邻居采样(比如每层各采 15 和 10 个邻居)能把显存压下来,代价是每轮看到的邻居不一样,训练曲线会抖,用大一点的 batch 和稍低的学习率能缓解。
第二个要点是过平滑与残差。两层 GCN 之后,如果发现所有节点的表示余弦相似度都在 0.9 以上,那就是过平滑了。解决办法有几种:加残差连接(每层输出和输入相加)、用 Jumping Knowledge 把各层表示拼起来、或者干脆减少层数。我一般先试减少层数,最快见效。
import torch import torch.nn.functional as F from torch_geometric.nn import SAGEConv class ProteinSAGE(torch.nn.Module): def __init__(self, in_dim, hidden=256, n_class=1, dropout=0.4): super().__init__() self.conv1 = SAGEConv(in_dim, hidden) self.conv2 = SAGEConv(hidden, hidden) self.dropout = dropout self.head = torch.nn.Linear(hidden, n_class) def forward(self, x, edge_index): h1 = F.relu(self.conv1(x, edge_index)) h1 = F.dropout(h1, p=self.dropout, training=self.training) h2 = F.relu(self.conv2(h1, edge_index)) # 残差连接,缓解过平滑 h2 = h2 + h1 return self.head(h2)训练配置上,优化器用 AdamW,学习率从 1e-3 起步,权重衰减 1e-4,batch 每个 epoch 重新采样一次邻居。损失函数用带pos_weight的二元交叉熵,pos_weight按该条目正负比的对数设,比直接设成比例更平缓。我一般训 50–100 个 epoch,在验证集 Fmax 上做早停,耐心值 10。
4.3 阈值选择与 Fmax 的正确用法
多标签分类的输出是概率,必须选阈值才能算出精确率和召回率。这里最常见的错误是用测试集选阈值——选出来的 Fmax 是乐观偏差的,而且偏差有多大完全取决于测试集大小。正确做法是在验证集上按每个条目单独扫阈值(从 0.01 到 0.99 步长 0.01),取让该条目 F1 最大的那个,然后把这个阈值冻结,应用到测试集上。
如果条件允许,更好的做法是分层扫阈值:按条目的正样本数量分箱,每个箱内共用一组阈值。这么做是因为正样本少的条目,其最优阈值往往偏低;正样本多的条目最优阈值偏高。全局单一阈值在这两类条目上都会吃亏。
import numpy as np from sklearn.metrics import precision_recall_curve def fmax_scan(y_true, y_score, n_steps=99): y_true = np.asarray(y_true).ravel() y_score = np.asarray(y_score).ravel() best_f, best_t = 0.0, 0.5 for t in np.linspace(0.01, 0.99, n_steps): pred = (y_score >= t).astype(int) tp = np.sum((pred == 1) & (y_true == 1)) fp = np.sum((pred == 1) & (y_true == 0)) fn = np.sum((pred == 0) & (y_true == 1)) if tp == 0: continue prec = tp / (tp + fp) rec = tp / (tp + fn) f = 2 * prec * rec / (prec + rec) if f > best_f: best_f, best_t = f, t return best_f, best_t提示:Fmax 是个偏乐观的指标,因为它在所有阈值上取最大值,相当于让模型自己挑了一个最有利的操作点。报结果的时候一定要同时给出固定阈值(比如 0.5)下的精确率和召回率,让读者知道实际部署时的大致表现。
5. 评估这套东西到底靠不靠谱
5.1 随机划分是最容易自欺的一步
我最初的做法是标准的三折随机划分,Fmax 报出来很漂亮。后来做了个对照——把序列相似性超过 30% 的蛋白聚成一簇,按簇划分训练和测试——Fmax 直接掉了十几个百分点。这个差距就是同源性泄漏造成的:随机划分下,测试集里的蛋白在训练集里有大量近亲,模型基本是在做记忆题。
如果你的目标是"预测全新蛋白的功能",那必须做按序列相似性聚类划分,用 MMseqs2 这类工具在 30% 一致性阈值下聚类,然后整簇分到训练或测试。如果你的目标是"给已有部分注释的蛋白补全注释",那可以按时间划分:用某年之前的注释训练,预测之后的注释。两种划分测的是完全不同的能力,报结果时必须说清楚用的是哪种。
5.2 指标怎么读才不被忽悠
三个分支的性能差异要分开报,混在一起报一个平均数是没意义的。另外 micro 和 macro 这两个口径必须分清:micro-AUC 把所有"蛋白-条目"对混在一起算,正样本多的条目会主导结果;macro-AUC 是每个条目单独算 AUC 再平均,能反映模型在小样本条目上的能力。两个都要报,差距越大说明模型对条目频次的依赖越强。
我实测的经验区间是:CC 分支最容易,Fmax 能到 0.6 以上;MF 次之,0.5 上下;BP 最难,0.4 左右。如果你看到某篇报道 BP 的 Fmax 超过 0.6 且没有特别的说明,大概率是划分方式偏乐观,或者标签空间用了传播后的完整注释但没提。这不是说做不到,而是要看清前提。
| 指标 | 计算口径 | 反映什么 | 常见误读 |
|---|---|---|---|
| micro-AUC | 全部蛋白-条目对混算 | 整体排序能力 | 被高频条目主导 |
| macro-AUC | 每条目单独 AUC 再平均 | 各条目平均表现 | 小样本条目噪声大 |
| AUPR | 精确率-召回率曲线下面积 | 不平衡下的表现 | 基线随正样本率变,不能跨数据集比 |
| Fmax | 最优阈值下的 F1 最大值 | 最佳操作点 | 偏乐观,必须配固定阈值结果 |
| 覆盖率@K | 前 K 个预测中命中比例 | 实际可用性 | K 的取值需贴合实际使用场景 |
5.3 一个必须做的对照:度基线
这是我最想强调的一条。只用 PageRank 或节点度作为唯一特征,训一个同样的模型,看它能拿到多少分。我做过这个实验,在部分 GO 条目上,仅用度基线就能拿到相当高的 micro-AUC。原因在于"研究热度"这个混杂因素:被研究得多的蛋白,一方面互作数据多(度大),另一方面注释也全(正样本多),两个信号高度相关,模型学到的是"热度"而不是"功能"。
如果你不做这个对照,很可能把一个只学会了识别热门蛋白的模型当成有效模型。正确的报告方式是把这个基线数字也摆出来,主模型必须显著超过它,这个结论才站得住。更狠一点的做法是在特征里显式加入"该蛋白的注释条目数量"作为混杂控制,然后看模型在控制之后还剩多少增益——通常会掉,但掉完之后剩下的那部分才是真正的功能预测能力。
6. 实操中踩过的坑与调优经验
6.1 低度节点上的性能骤降与分层报告
新拿到的候选蛋白,在网络里的度往往是个位数,甚至为零。这类节点恰恰是最需要预测的,却也是模型表现最差的。我做过按度分箱的评估:度大于 50 的节点,CC 分支 Fmax 能到 0.7 以上;度在 1–5 之间的节点,掉到 0.35 左右;完全孤立的节点,模型基本退化成只用序列特征。
这个结果是正常的,不应该试图掩盖。我的处理有三条:一是在报告里按度分箱给出性能,让使用者知道在什么情况下该信这个预测;二是对孤立节点降级处理,直接输出"基于序列特征的预测,置信度低";三是在预测阶段回退到传播算法——传播算法对孤立节点的处理虽然也弱,但至少不会给出过度自信的错误答案,两者交叉验证可以筛掉一批明显不一致的预测。
6.2 标签泄漏的隐蔽来源
标签泄漏这个话题值得单独说,因为它的表现形式太隐蔽了。第一种是GO 注释本身来自互作证据,有些注释的证据码明确写了"由互作推断而来",如果你把它当正样本、又用互作网络做特征,那就是标准的循环论证。处理办法是过滤掉证据码里标明推断来源的注释,只用实验证据。
第二种是同源泄漏,前面讲过,靠按簇划分解决。第三种比较阴险:特征里混入了未来信息。比如你用某个数据库的"当前版本"来映射标识符,而这个版本的映射关系包含了测试集中蛋白的新命名,等于模型间接看到了答案。解决办法是把所有映射表、注释文件都按时间戳冻结,测试集只用训练时刻之前可得的信息。第四种是预处理阶段的全局泄漏,比如在划分之前就用全量数据做了 PCA 或标准化,测试集的分布信息漏进了训练。这类问题的修法很简单:任何 fit 操作都必须在训练集上做,然后 transform 测试集。
6.3 从实验代码到能跑起来的流程
实验做得再好,如果跑一次要人工改五个路径,那它就永远不会被用起来。我后来把整个流程拆成四个阶段:构建图与冻结映射、计算特征并缓存、训练与评估、推理与报告。每个阶段的输入输出都是磁盘上的文件,中间不共享内存。好处是可以单独重跑某一阶段——比如只更新了注释数据,就只重跑特征和训练,不用重新清洗整张图。
工程上有几个具体的优化点值得记下来。传播那步的矩阵乘法是性能瓶颈,用scipy.sparse的 CSR 格式比 CSC 快不少,但要注意P @ Y里 Y 是稠密的,这一步的内存占用是节点数乘以条目数再乘 4 字节,两万乘两千大约是 160MB,可以接受;如果条目数上万,就得按块处理。拓扑特征里的介数中心性用采样近似,我实测 500 个源点的结果和精确值在排序上的一致性已经足够。最后,所有随机过程必须固定种子并记录在配置里,包括负采样、邻居采样、初始化和数据划分,否则你永远无法复现一周前那个漂亮的数字。
注意:整个链路里最有价值的不是模型本身,而是那份冻结的映射表、缓存的特征文件和一键重跑的脚本。换个人接手时,这些东西决定了他是一天能上手还是一个月都在猜。
6.4 一个容易被忽略的后处理:层级一致性
GO 是 DAG 结构,预测结果必须满足 true path rule——如果模型预测某个蛋白具有某个具体条目,那它必然也应该具有该条目的所有祖先条目。原始预测经常违反这个约束,表现为一个节点的子条目概率高、父条目概率低。做一次自底向上的概率最大化传播(子节点的预测值向父节点传递取最大值),能把整体的 Fmax 抬 2–5 个百分点,几乎零成本。
这一步的意义不只是涨点。实际交付时,使用者看到的是一棵功能树,如果树上的结论自相矛盾,信任度立刻崩塌。我见过最尴尬的情况是模型预测某蛋白具有"某个非常具体的代谢过程",但没有"代谢过程"这个顶层条目——这种结果拿出去没法解释。加上层级一致性后处理,这类低级矛盾就消失了。
一些个人的体会
这套流程我前后迭代了大半年,最大的感受是:模型那部分其实是最不重要的。真正决定结果质量的是数据清洗的严谨程度、划分方式的诚实程度、以及对照实验做得够不够。我见过太多项目在模型结构上反复折腾,从树模型换到各种图网络,涨了一个百分点就兴高采烈,但划分方式还是随机划分,度基线从来不做,这样的数字没法复现也没法用。
如果让我给一个最实用的建议:先花两天把度基线跑出来,再花一天把按簇划分的评估框架搭好,这两件事做完,你对这套方法的能力边界就有了真实判断,后面所有的改进都有了参照系。至于用树模型还是图网络,用随机游走还是语言模型嵌入,那都是在这个参照系之上的边际优化,早做晚做都行。
另外一个常被忽略的点是结果的可解释性。预测一个新蛋白属于某个复合体,光给个概率不够,我在报告里会附上支撑这个预测的关键邻居——也就是传播得分贡献最大的那几个已知蛋白。这一步不增加任何计算成本,但对使用者来说,"它和这几个已知成员有互作"比"模型给出 0.87 的概率"要有说服力得多,也更容易被发现错误:如果这几个邻居本身注释就是错的,一眼就能看出来。