☰
蛋白质二级结构预测:CNN与Transformer实战与避坑指南
2026/10/1 12:35:30 网站建设 项目流程

简介:面向高校学生、生物信息学初学者及深度学习实践者的Python项目源码,基于CNN(卷积神经网络)或Transformer模型实现蛋白质二级结构预测。内容覆盖数据预处理、模型封装、训练、预测与结果保存等完整流程,适合作为期末大作业、课程设计或入门实战的可运行参考。压缩包共12个文件,主要包括5个Python脚本(对应模型定义、数据读取、训练与推理等环节)、2个预训练参数文件、3个文本说明/结果文件及若干缓存文件,总大小约43.76MB,目录按模型、参数、结果分区,结构清晰便于定位。项目已在本地编译通过,评审得分95分以上,难度适中;预置ResNet编码器与Transformer编码器两套参数,免去自行训练的时间,可直接运行测试,也可对比不同网络结构在二级结构预测任务上的效果。目前已有227人学习下载,对有课程设计或期末项目需求的读者来说,是一份完成度高、可直接复用的实战资料。

1. 蛋白质二级结构预测为什么选 CNN 或 Transformer:这个项目到底在做什么

蛋白质二级结构预测是典型的序列标注任务:输入一段氨基酸序列,输出每个残基对应的结构类别(常用三类 H/E/C,或 DSSP 八类)。用 Python 做这件事,方案基本收敛到两条技术路线——CNN 卷积神经网络和 Transformer 自注意力架构。做课设、毕设或入门深度学习在生物信息学上的落地,这个方向既能出可视化结果,又有明确的量化指标(Q3/Q8 准确率),很适合作为“高分实战可运行项目”来打磨。这篇笔记按我自己的实操顺序,把模型选型、数据构造、训练调参和答辩前必查的坑完整过一遍。

2. 先定模型再写代码:CNN 与 Transformer 的选型依据和直观差异

2.1 从序列到标签:二级结构预测是个逐残基分类问题

先把问题框死。一条蛋白质序列,比如长度为 100 的氨基酸链,输入到项目里后,期望的输出是长度同样为 100 的标签序列,每个位置标注这个残基属于 helix(螺旋)、sheet(折叠)还是 coil(无规卷曲)。这和你做过的一些 NLP 任务在结构上是一致的——命名实体识别也是这种逐 token 分类,只不过这里的“词”是氨基酸,上面没有文本语义,只有理化性质和进化信息。

这里有个新手容易绕进去的弯:既然输入是序列,为什么不直接把它当成文本,套个现成的 BERT 来用?答案是:标准蛋白质序列长度偏短(几百个残基),局部结构模式非常显著,而语言模型依赖的上下文语法规律在氨基酸序列上天然弱化。CNN 和 Transformer 都是合适的底座,但选型的逻辑不一样。

2.2 CNN 方案:一维卷积天然适合短窗口的局部模式

CNN 是二级结构预测里的经典选手。PSIPRED 这类早期工具本质上就是滑窗特征配合浅层神经网络,而现代做法是直接把序列转换成特征矩阵,用一维卷积堆叠来捕捉局部残基互作。一维卷积在长度维度上滑动,每个卷积核相当于一个短肽“模板”,比如窗口大小为 11 的卷积核看到的是一个 11 肽的局部模式,这对 alpha-helix 这种高度局部化的结构特别有效。

我一般会把 CNN 作为基线模型先跑通全流程。参数少、收敛快、对小数据集友好,在几百条蛋白质的训练集上也能收敛到不错的 Q3。如果项目时间紧,目标是快速打通端到端流程,CNN 几乎是唯一不会让你心态崩掉的选择。它的天花板在于长程依赖——一个残基的二级结构有时受 30 个残基之外的氨基酸影响,CNN 只能靠堆层数扩大感受野,效果和代价都上去了。

2.3 Transformer 方案:自注意力能抓到长程依赖但代价更大

Transformer 自注意力机制的优势正好对应在 CNN 最吃力的地方。注意力矩阵是全局的,第 i 个残基可以和序列内任何位置的残基直接建立关系,不需要像卷积那样一层层传递信息。这对那些二硫键约束、远端互作形成的结构片段尤其重要。代价也明显:参数量大、需要更多数据、在小训练集上很容易过拟合,训练时 loss 稳定下来的速度比 CNN 慢得多。

如果项目数据量充足(比如公开数据集用于训练的蛋白质超过一万条),Transformer 编码器结构就值得认真调。通常我只会用 Encoder 部分,把 d_model 设成 64 或 128,加上 4 到 8 层编码器,输出端接一个线性层做逐位置分类。位置编码不能省,氨基酸顺序信息全靠它;训练时还要加 dropout,不然训练集准确率会冲到 95% 以上而验证集纹丝不动。

对比项CNN 基线Transformer 编码器
训练速度快慢
参数量小,适合几百条蛋白质的小数据集大,需要数据量支撑
长程依赖弱,靠层数堆感受野强,自注意力直接全局建模
过拟合风险低高,必须配 dropout 和早停
可解释性卷积核模式较难直观表述注意力权重可直接可视化

3. 数据准备是 80% 的功夫:PSSM 特征、DSSP 标签与数据集划分

3.1 第一步:把氨基酸序列变成模型能吃的特征矩阵

不要一上来就搭模型,数据特征没做好,后面所有分数都是虚的。二级结构预测的特征矩阵常见做法由三部分拼接而成:PSSM 位置特异性打分矩阵、氨基酸独热编码、以及可选的理化性质特征。PSSM 是关键项,它记录了每个位置上的氨基酸在进化上被替换成其他氨基酸的倾向性,相当于把同源序列的进化信息压缩进了特征里,对预测效果的提升非常明显。

PSSM 生成依赖 PSI-BLAST 这类序列比对工具,比对到蛋白质数据库后输出一个 L×20 的矩阵。纯 Python 不负责生成 PSSM,而是负责把工具输出解析成干净的特征文件。这一步较慢,通常要跑一段时间,是项目里最容易被忽视的“时间黑洞”。模型训练反而很快。

import numpy as np def parse_pssm_to_matrix(pssm_path, seq_len): """ 把 PSI-BLAST 输出的 PSSM 文本解析为 L x 20 的 float 矩阵 这类工具输出的原始文本带表头表尾,需要跳过无关行。 """ matrix = np.zeros((seq_len, 20), dtype=np.float32) with open(pssm_path, "r", encoding="utf-8", errors="ignore") as f: lines = f.readlines() row_idx = 0 for line in lines: parts = line.split() # 判断是不是数据行:开头是残基序号,第二列是氨基酸字母 if len(parts) >= 22 and parts[0].strip().isdigit(): try: scores = [float(x) for x in parts[2:22]] matrix[row_idx, :] = scores row_idx += 1 except ValueError: continue # 归一化:PSSM 原始数值浮动范围大,按行做 z-score 后模型更稳定 mean = matrix.mean(axis=1, keepdims=True) std = matrix.std(axis=1, keepdims=True) + 1e-6 matrix = (matrix - mean) / std return matrix

逻辑说明:PSSM 输出文件的行格式在不同版本工具中略有差异,但核心规律一致——中间有 20 列数值,对应 20 种标准氨基酸的打分。解析时先判断行首是否为数字,避开表头表尾;取第 3 到第 22 列是通用做法。归一化这步特别重要,直接喂原始分数会让模型的权重初始化失去意义,收敛速度显著变慢。

参数说明:seq_len必须和实际序列长度一致,建议在预处理阶段维护一个{protein_id: length}字典,避免矩阵维度对不上。+1e-6是防止某行全为零时除零报错,这种行在序列两端常出现。

3.2 第二步:用 DSSP 生成逐残基标签并处理类别不平衡

标签的标准来源是 DSSP 工具,它根据蛋白质三维结构(PDB 文件)计算每个残基的二级结构归属,通常输出八类别(3_10-helix、alpha-helix、beta-sheet 等)。训练时既可以精调到 Q8(八类准确率),也可以合并成 Q3(三类准确率)方便和论文对比。我的建议是两者都算,答辩时 Q8 更能说明模型的判别粒度。

def dssp_8class_to_3class(mapping_dict): """ DSSP 八类标签合并为三类,用于 Q3 评测。 合并规则是常见做法:G/H/I 归为 H,E/B 归为 E,其余归为 C。 """ label_8_to_3 = { "H": "H", "G": "H", "I": "H", "E": "E", "B": "E", "C": "C", "T": "C", "S": "C" } seq_labels = [] for aa in mapping_dict: seq_labels.append(label_8_to_3.get(aa, "C")) return "".join(seq_labels)

逻辑说明:DSSP 的八类标签中,H 是 alpha 螺旋、G 是 3_10 螺旋、I 是 pi 螺旋,三者结构上都算螺旋族;E 是 beta 折叠、B 是孤立 beta 桥,合并为 sheet;其余 loop、turn、bend 全部归为 coil。这个映射是学术界的标准做法,不是拍脑袋定的。

参数说明:如果模型输出八类,则最后一层全连接输出 8 个神经元;三类则输出 3 个神经元。两者对应不同的损失函数配置,训练前一定要确认标签编码和模型输出维度一致,不然训练循环第一轮就报 mismatch 错误。

3.3 第三步:按蛋白质划分数据集,避免序列泄漏导致虚高分数

这是整个项目里最容易拿高分又最容易翻车的环节。很多教程会把所有蛋白质的滑窗片段混合在一起,随机按 8:2 划分训练集和测试集——这样做出来的验证准确率可以虚高 10 个百分点以上,因为同一条蛋白质和它的同源序列片段会同时出现在训练集和测试集里,模型实际上在“背答案”,而不是在泛化。

正确做法是以蛋白质为单位划分。每一条蛋白质的全部滑窗数据要么全在训练集,要么全在测试集,不允许交叉。如果想让分数更有说服力,可以按序列相似性聚类后再划分,但把这层逻辑写进课程项目里已经足够赢得评审认可。

from sklearn.model_selection import train_test_split # protein_ids 是所有蛋白质 ID 列表 # 按蛋白质粒度切分,而不是按滑窗样本切分 protein_train, protein_test = train_test_split( protein_ids, test_size=0.2, random_state=42 ) train_samples = [s for s in all_samples if s["protein_id"] in protein_train] test_samples = [s for s in all_samples if s["protein_id"] in protein_test]

逻辑说明:all_samples是滑窗化之后的数据列表,每个元素记录了它来自哪条蛋白质。切分时先对蛋白质 ID 集合做随机划分,再用集合成员关系过滤样本。random_state=42固定随机种子,保证每次复现结果一致,这是答辩时“可复现实验记录”的必要条件。

参数说明:test_size=0.2是常用比例,但如果总蛋白质数少于 100 条,建议降到 0.15 甚至 0.1,否则测试集蛋白质太少,单条异常序列会把整体准确率拉低好几个点。

4. 从零搭起最小可用项目:CNN 与 Transformer 的 PyTorch 实现要点

4.1 搭建 CNN 基线的 60 行代码和 5 个关键参数

CNN 基线的目标不是刷分,而是验证全流程能跑通。我用 PyTorch 写了一个最小结构:嵌入层不做,直接用特征矩阵作为输入;一层一维卷积加 BatchNorm 加 ReLU;接一个全局最大池化(在长度维度上);最后过全连接层输出逐位置分类。这里不做滑窗,而是把整条序列的特征矩阵直接输入,让卷积核跨位置滑动,输出同样长度的预测序列。

import torch.nn as nn class CNNBaseline(nn.Module): def __init__(self, in_channels=22, num_classes=3): super().__init__() # in_channels = PSSM 的 20 维 + 独热编码的 20 维 + 理化性质 2 维(示例) self.conv1 = nn.Conv1d(in_channels, 64, kernel_size=11, padding=5) self.bn1 = nn.BatchNorm1d(64) self.conv2 = nn.Conv1d(64, 128, kernel_size=7, padding=3) self.bn2 = nn.BatchNorm1d(128) self.dropout = nn.Dropout(0.3) # 逐位置分类头,kernel_size=1 等价于在序列长度维度上共享权重 self.classifier = nn.Conv1d(128, num_classes, kernel_size=1) def forward(self, x): # x 形状: (batch, seq_len, in_channels) -> 转为 (batch, in_channels, seq_len) x = x.transpose(1, 2) x = torch.relu(self.bn1(self.conv1(x))) x = self.dropout(x) x = torch.relu(self.bn2(self.conv2(x))) x = self.dropout(x) # 输出形状: (batch, num_classes, seq_len) -> 转回 (batch, seq_len, num_classes) return self.classifier(x).transpose(1, 2)

逻辑说明:Conv1d 默认在最后一维上滑动,所以输入先转成通道在前的格式。padding=5配合kernel_size=11是为了保持输出长度和输入长度一致,这是序列标注任务的关键——预测序列必须和输入序列一样长。最后的分类头用kernel_size=1的卷积,等价于在序列每个位置上做一次共享的全连接变换,既减少参数量又避免序列长度对全连接层输入维度的限制。

参数说明:5 个关键参数分别是in_channels(特征维度)、kernel_size(感受野)、hidden_dim(通道数 64/128)、dropout、num_classes。kernel_size在二级结构预测里常用 9~15,太小时感受野不足,太大时参数量上涨且容易过拟合;dropout=0.3是 CNN 基线的稳妥值,数据量越小越要往上调。

4.2 把 CNN 换成 Transformer 编码器:位置编码与注意力掩码怎么写

Transformer 方案我直接用了torch.nn.TransformerEncoderLayer作为组件,没有自己手写多头注意力,但位置编码必须自己加。标准做法是正弦位置编码,核心思想是让每个位置有一个唯一的编码向量,模型才能区分“第 3 个残基”和“第 7 个残基”。这里要特别提醒:序列长度不一,位置编码必须能泛化到训练时未见过的长度,正弦编码满足这个条件。

import torch import torch.nn as nn import math class PositionalEncoding(nn.Module): def __init__(self, d_model, max_len=512): super().__init__() pe = torch.zeros(max_len, d_model) position = torch.arange(0, max_len, dtype=torch.float).unsqueeze(1) div_term = torch.exp(torch.arange(0, d_model, 2).float() * (-math.log(10000.0) / d_model)) pe[:, 0::2] = torch.sin(position * div_term) pe[:, 1::2] = torch.cos(position * div_term) pe = pe.unsqueeze(0) # (1, max_len, d_model) self.register_buffer("pe", pe) def forward(self, x): return x + self.pe[:, :x.size(1), :] class TransformerBaseline(nn.Module): def __init__(self, in_channels=22, d_model=128, num_heads=8, num_layers=4, num_classes=3): super().__init__() self.input_proj = nn.Linear(in_channels, d_model) self.pos_enc = PositionalEncoding(d_model, max_len=512) encoder_layer = nn.TransformerEncoderLayer( d_model=d_model, nhead=num_heads, dim_feedforward=512, dropout=0.3, batch_first=True ) self.encoder = nn.TransformerEncoder(encoder_layer, num_layers=num_layers) self.classifier = nn.Linear(d_model, num_classes) def forward(self, x): # x: (batch, seq_len, in_channels) x = self.input_proj(x) x = self.pos_enc(x) x = self.encoder(x) return self.classifier(x)

逻辑说明:TransformerEncoderLayer 在批量训练时要求同批次序列长度一致,所以数据加载器里要做 padding,用一个占位符向量填充短序列到统一长度。上面的代码没有显式传入 padding mask,实际训练时必须把src_key_padding_mask传进去,否则模型会把 padding 位置当作有效残基来算注意力,预测结果里会出现一排无意义的边缘分数。

参数说明:d_model=128是序列特征映射到注意力空间的维度;num_heads=8是常规配置,但如果数据集小,4 个头更不容易过拟合;num_layers=4对二级结构预测基本够用,堆到 8 层以上在小数据集上纯粹是折磨优化器。max_len=512表示最长支持 512 个残基,长于这个长度的蛋白质在推理时要拆段。

4.3 训练循环里必须监控的四个指标

训练循环本身不复杂,难的是懂得看哪些数字。除了常规的 train loss,四个指标必须同步打印:验证集 Q3 准确率、验证集平均损失、类别维度的召回率、以及训练 loss 和验证 loss 的差值。这四个数字能直接告诉你模型是在正常学习还是已经过拟合。

def train_one_epoch(model, dataloader, optimizer, criterion, device): model.train() total_loss = 0.0 for batch in dataloader: x = batch["feature"].to(device) y = batch["label"].to(device) mask = batch["mask"].to(device) logits = model(x) # 只计算真实残基位置的损失,忽略 padding 位置 loss = criterion(logits.view(-1, logits.size(-1)), y.view(-1)) loss = loss * mask.view(-1).float() loss = loss.sum() / mask.sum() optimizer.zero_grad() loss.backward() optimizer.step() total_loss += loss.item() * mask.sum().item() return total_loss / len(dataloader.dataset)

逻辑说明:mask是数据加载器里同步生成的布尔张量,标记哪些位置是真实残基、哪些是 padding。计算损失时先把预测和目标都展平,逐位置算 CrossEntropy,然后用 mask 把 padding 位置的损失清零,最后除以真实位置的数量做平均。如果不加 mask,padding 位置会被模型强行预测成某个类别,这部分损失会干扰梯度,导致验证准确率偏低。

参数说明:优化器我一般用 AdamW,学习率 1e-3 起步,配 CosineAnnealingLR 调度器;如果 loss 震荡严重,降到 5e-4 再做 warmup。batch size 在序列任务中不用开太大,16 到 32 之间通常足够,太大反而让 BatchNorm 在短序列上的统计量不稳定。

5. 避坑指南:从读不懂损失到测试集分数虚高,5 次翻车记录

5.1 加了类别权重后,验证指标反而下降

现象:C 类(coil)样本数量多,模型倾向于把所有位置都预测成 C,导致 H 类和 E 类的召回率极低。我一开始给 CrossEntropyLoss 加了类别权重,把 H 和 E 的权重调大,结果 Q3 不升反降。

原因:类别权重确实提高了少数类的召回率,但副作用是模型把大量模糊的 C 类样本误判成了 H 或 E,整体 Q3 被拉低。Q3 是所有位置的平均准确率,少数类权重过高会牺牲多数类的精度。

解决:改成先不加权重跑一版,看混淆矩阵再决定。如果 H/E 召回率实在不能看,把权重控制在 1.2 倍以内,同时调整分类阈值而不是直接调损失权重。在答辩项目里,给出“类别不平衡如何影响 Q3”的分析比盲目上权重更有说服力。

5.2 PSSM 生成在 Windows 上跑不起来,或比对速度慢到让人想放弃

现象:在 Windows 环境下执行 PSI-BLAST 命令提示找不到可执行文件;换成 Linux 后单条序列比对耗时达几十秒,几百条蛋白质根本不敢想。

原因:PSI-BLAST 依赖本地序列数据库,Windows 下需要配置环境变量且数据库格式必须是makeblastdb处理过的;耗时长的原因通常是数据库太大或未限制迭代次数。

解决:生成 PSSM 时统一用-num_iterations 3限制迭代次数,这是二级结构预测的常用值,再多迭代对特征提升有限而耗时翻倍。数据库用小型非冗余库(比如 NR 的子集)即可满足教学项目需求。另外把这个步骤做成离线前置任务,批量生成后保存为.npy文件,训练时直接加载,不要再在训练脚本里实时调用比对工具。

5.3 Transformer 训练到中途 loss 变成 NaN

现象:训练到第 20 个 epoch 左右,loss 突然变成 NaN,之前一切正常。

原因:自注意力机制里点积结果的数值范围随序列长度增长而变大,softmax 之后出现梯度爆炸;d_model 较小而学习率偏高时尤其容易触发。

解决:先降低学习率到原来的十分之一重跑;如果还复现,检查dim_feedforward是否过大导致的激活值饱和。绕不开的一步是在优化器更新后用torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)做梯度裁剪,这对 Transformer 是保命操作。我后面所有训练脚本都默认带梯度裁剪,损失函数再也没出现过 NaN。

5.4 评估指标算出来 88%,但换一拨序列就打回原形

现象:项目自带的测试集上 Q3 达到 88%,同学拿一条新的蛋白质序列来测,预测结果和实际结构对不上,分数掉到 70% 不到。

原因:最典型的是数据泄漏——切分训练集和测试集时没有按蛋白质划分,同源序列同时出现在两边。另一种可能是测试集蛋白质和训练集蛋白质来源相同(比如同一个蛋白家族),模型记住了家族特征而不是通用结构规则。

解决:按第 3.3 节的做法做蛋白质级别切分,并检查测试集蛋白质 ID 与训练集有无重叠。如果条件允许,用一组外部蛋白质序列做二次验证,报告里写清楚这个验证结果,反而比单一测试集分数更经得起评审追问。

5.5 直接用 FASTA 序列做预测,效果和论文结果差 20 个百分点

现象:模型在带 PSSM 特征的测试集上 Q3 有 80% 以上,但用一段新 FASTA 序列直接用独热编码喂给模型,准确率掉得没法看。

原因:模型训练时输入包含了 PSSM 特征,这些进化信息承载了大部分预测能力。新序列没有生成 PSSM 就直接输入,相当于把特征从 20 维进化打分变成 20 维独热编码,分布完全变了。

解决:推理阶段必须和训练阶段走同一条特征流水线。新序列先用 PSI-BLAST 生成 PSSM,再做同样的归一化,最后拼上独热编码才能输入模型。我自己的项目脚本里把“FASTA -> PSSM -> 特征拼接 -> 模型推理”封装成了一条命令,避免手动步骤遗漏。

6. 跑通只是起点:用 Q3/Q8 指标做评审,把注意力图做成可解释结果

6.1 评审前用这几张图表汇报:混淆矩阵、逐位置准确率、注意力热度图

课程项目拿高分的关键不在模型多复杂,而在验证是否完整。Q3/Q8 数字只是第一步,我建议额外输出三张图。第一张是混淆矩阵,展示 H/E/C 三类互相混淆的情况,评审一眼就能看出模型在哪些结构上区分度不足。第二张是逐位置准确率曲线,统计不同序列位置上预测正确率的变化,通常序列两端准确率偏低,因为卷积或注意力在这些位置能利用的上下文更少,这个现象可以提前在答辩里讲清楚。

Transformer 分支可以增加第三张图:注意力热度图。把某个查询位置对所有位置的注意力权重画成热力图,能直观显示第 i 个残基主要关注了哪些远端残基。一级章里提到的“自注意力能抓长程依赖”,在这张图上得到直接可视化验证,比任何口头解释都有力。

6.2 序列窗口滑动推理:长序列不截断的批处理写法

推理阶段会遇到长序列超出模型输入长度限制的情况,直接把序列截断会丢掉远端依赖。常见做法是用滑动窗口配合重叠拼接:窗口大小为 256,步长为 128,两个相邻窗口重叠部分取两次预测的平均值。这样既控制单次计算量,又让重叠区域的预测利用了两侧上下文,避免窗口边界处准确率下降。

6.3 保存 checkpoint 与可复现实验记录

训练结束时保存完整的信息,不只是模型权重,包括特征归一化参数、类别编号映射、PSSM 解析参数和数据集划分的随机种子。这几个因素里任何一个对不上,复现出来的分数都可能有明显出入。我自己的习惯是把训练命令和每次实验的指标写进一个实验记录文件,哪怕后来只改了一个 dropout,也会单独记一行。这个习惯在答辩时帮了大忙,评审问到任何参数都能当场说清。

这个项目做完之后让我最深的体会是:二级结构预测的模型选型其实只决定分数的上限,数据特征和评测方式才决定这个上限能不能兑现。如果你打算在这个方向上投入,先花时间把 PSSM 特征和按蛋白质划分数据集这两件事做好,再回头调模型结构,路径会顺很多。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询