上次我们用CNN重做了一遍RNN对电影评论分类 ,这次我们用RNN来重做一遍CNN识别DNA启动子 。
这两次"交换作业"是配套的:在《使用CNN对电影评论进行分类 》里,CNN抢了RNN的活;这一篇反过来,让LSTM去干CNN最擅长的活。只有把同一个任务交给两个模型各跑一遍,才能看清它们各自的适用边界在哪里。
CNN 用卷积核在序列上滑动,天然擅长 局部、位置无关 的模式。TATA box 正好是一个长度 6bp 的固定串,出现在哪个位置都算——这几乎是给卷积量身定做的任务。 LSTM 用循环单元一个碱基一个碱基地读过去,靠隐藏状态"记住"前面看过什么,它的强项是 长程依赖和变长序列 (见《循环神经网络RNN 》)。 所以这一篇的预期结果不是"RNN 更好",而是" RNN 也能做,但要付出更多代价 "。具体差多少,文末有对照。
import random import torch import torch.nn as nn import torch.optim as optim from torch.utils.data import DataLoader, TensorDataset 数据 还是一样地产生随机序列,这次motif加多一个,有两个变体,这样更真实一些。
def generate_data ( n= 3000 , length= 100 , motifs=( "TATAAA" , "TATACA" ) ): random.seed( 42 ) pos, neg = [], [] bases = "ACGT" for _ in range (n // 2 ): s = "" .join(random.choice(bases) for _ in
range (length)) m = random.choice(motifs) i = random.randint( 0 , length - len (m)) pos.append(s[:i] + m + s[i + len (m):]) s = "" .join(random.choice(bases) for _ in range (length)) while any (m in s for m in motifs): s= "" .join(random.choice(bases) for _ in range (length)) neg.append(s) X = seq_to_indices(pos + neg) y = torch.tensor([ 1 ]* len (pos) + [ 0 ]* len (neg), dtype=torch.long) perm = torch.randperm( len (y)) return X[perm], y[perm] 正样本是"替换",不是"插入"。 先造一条 100bp 的随机序列,再随机挑一个起点 i ,用 motif 覆盖掉原来的 6 个碱基( s[:i] + m + s[i+len(m):] )。因为是替换,序列长度仍然是 100。 i 的取值是 randint(0, length - len(m)) ,也就是 0~94,保证替换不会越界。
负样本的 while 循环不能省。 负样本要求"碰巧不含 motif",一旦发现就重新生成。随机序列碰巧含有某个变体的概率其实不算小:长度 100 有 95 个可能的起点,每个起点命中某个具体 6-mer 的概率是 ,两个变体合起来期望个数约
,也就是大约 4.5% 的随机序列会中招。如果不做这个过滤,1500 条负样本里会有 60 多条被贴上错误的标签 0,占全部样本的 2% 左右—— 准确率的上限就被钉死在 97.7% 附近了 ,后面无论怎么调模型都上不去。这类"标签噪声设了天花板"的现象,是排查"模型训不动"时经常被忽略的一环。
两个变体是刻意加的难度。 TATAAA 和 TATACA 只差第 5 位。真实的 TATA box 本身就是简并的(degenerate),不同基因的序列并不完全一样,所以加一个"错一个碱基"的变体更贴近真实情况。副作用是任务变难了:模型不能只记死一个字符串,得学到"前四位是 TATA、后两位是 A,中间可容忍一处不匹配"这种更抽象的模式。这是本文准确率明显低于 CNN 篇的原因之一——但注意,这属于 数据难度 的差异,不是模型能力的差异。
X[perm], y[perm] 用的是同一个 perm 。 这一点必须小心:如果特征和标签各用一次 randperm ,样本和标签就错位了,模型学到的是纯噪声,准确率会永远停在 50% 附近。
顺带一个可读性上的小问题: while any(m in s for m in motifs) 里的 m ,和上面 m = random.choice(motifs) 的 m 名字撞车了。好在生成器表达式有独立作用域,不会污染外层的 m ,所以代码是对的;但读起来容易误会,改成 for mt in motifs 会清爽很多。
另外 n=3000
正好是一半正一半负(1500 + 1500),类别平衡,所以后面用准确率(accuracy)当指标是合理的。如果类别不平衡,就该看 F1 或 AUC 了。
把序列变成索引 然后是把序列变成数字,这是为了满足后面转成Embedding向量。
def seq_to_indices ( sequences ): mapping = { "A" : 0 , "C" : 1 , "G" : 2 , "T" : 3 } x = [[mapping[base] for base in seq] for seq in sequences] return torch.tensor(x) 输出是一个形状为 (N, L) 、类型为 torch.long 的整数张量—— A/C/G/T 各自对应 0/1/2/3。
为什么是整数而不是独热编码? 因为下一步用的是 nn.Embedding 。 nn.Embedding 只接受索引,它内部维护一张形状为 (4, 32) 的查表矩阵,输入编号 k 就取出第 k 行那个 32 维向量。换句话说, Embedding 层本身就是"可学习的独热编码" :独热编码是固定死的(某一维为 1,其余为 0),而 Embedding 让模型自己决定每个碱基该用什么向量表示,顺便把维度从 4 升到了 32。
对照 CNN 那一篇:那边我们手工做了 (N, 4, L) 的独热编码,因为 nn.Conv1d 要的是"通道 × 长度"这种像图片一样的输入;RNN 这边要的是"一个位置一个 token"的序列,所以用索引 + Embedding。 同一份序列数据,喂给不同模型时的格式并不一样 。
( mapping 里谁编 0 谁编 1 其实无所谓,只要全局一致、并且和 Embedding 的词表大小 4 对得上就行。)
模型 class LSTMClassifier (nn.Module): def __init__ ( self ): super ().__init__() self.embedding = nn.Embedding( 4 , 32 ) self.lstm = nn.LSTM( 32 , 64 , batch_first= True ) self.fc = nn.Linear( 64 , 2 ) def forward ( self, x ): x = self.embedding(x) out, _ = self.lstm(x) x = out.mean(dim= 1 ) x = self.fc(x) return x
词表的大小是4,转成Embedding,变成32维向量。输入LSTM是32维,输出隐藏64维。全连接层吃64维,出来2维。
__init__ 里三层各自的职责:
nn.Embedding(4, 32) :词表大小 4,嵌入维度 32。权重矩阵形状 (4, 32) ,一共只有 128 个参数,是个非常小的查表。 nn.LSTM(32, 64, batch_first=True) : input_size=32 接 Embedding 的输出, hidden_size=64 。 batch_first=True 是个高频坑 ——它表示输入形状是 (N, L, 32) (批次在前),而不是 PyTorch 默认的 (L, N, 32) 。忘了写这个参数代码不会报错,但会把"批次"当成"序列长度",训练结果完全错乱。 nn.Linear(64, 2) :把 64 维的序列表示映射成 2 个类别的分数(logits),交给 CrossEntropyLoss 。 forward 里的形状变化,一路写清楚:
x (N, 100) self.embedding(x) (N, 100, 32) self.lstm(x) (N, 100, 64) out.mean(dim=1) (N, 64) self.fc(x) (N, 2)
LSTM的输出是 (out, (h, c)) ,out是 (N, L, 64) ,每一个位置(碱基)都有64维的记忆状态。这里把L(长度)那一维的元素求平均,这就是所谓的平均池化。这就是 out.mean(dim=1) 这一句代码干的事情。
展开说一下这三个返回值:
out :形状 (N, L, 64) , 每一个时间步 的隐藏状态都保留着。L=100,所以每个碱基位置都有一个 64 维的"记忆"。 h_n :形状 (1, N, 64) ,只有 最后一个时间步 的隐藏状态(因为 num_layers=1 且是单向,第一维是 1)。 c_n :LSTM 特有的记忆单元(cell state)的最终值,形状同上。 这里只用了 out ,另外两个被 _ 丢掉了,所以写成 out, _ = self.lstm(x) 。
那为什么用 out.mean(dim=1) 而不是 out[:, -1, :] ? 三种做法都常见,但适用场景不同:
取最后一个时间步 ( out[:, -1, :] ):假设"读完整条序列后,最终的记忆包含了全部信息"。这对情感分类那种"读完整句话下判断"的任务很自然,但对"序列里某处藏了个 motif"就有点别扭——motif 要是在前半段,它的信息得一路经过几十个时间步才能传到末尾,中间还有被遗忘门丢掉的风险。 对所有时间步求平均 (本代码的做法):相当于"把整条序列的记忆平均一下"。好处是每个位置的梯度都能直接回传,motif 无论出现在哪里贡献都差不多;坏处是 motif 只占 100 个位置里的 6 个,信号被 95% 的背景稀释了。 最大池化 ( out.max(dim=1).values ):只看"哪个位置反应最强烈"。对"检测某个局部模式是否存在"这类问题,这才是最对路的池化方式——它在问"最像 motif 的那个位置有多像",而平均池化在问"整条序列平均下来长什么样"。 顺带一提,这里的 LSTM 是 单向 的,读到第 i 个碱基时只知道它左边发生了什么。motif 检测并不要求因果性(我们不是预测下一个碱基,而是读完一整条再判断),所以换成 bidirectional=True 会更自然,每个位置能同时看到左右两侧的上下文——代价是 fc 的输入维度要改成 64 × 2 = 128。
训练与评估
def train_one_epoch ( model, loader, optimizer, criterion, device ): model.train() total_loss = 0.0 for x, y in loader: x, y = x.to(device), y.to(device) optimizer.zero_grad() outputs = model(x) loss = criterion(outputs, y) loss.backward() optimizer.step() total_loss += loss.item() return total_loss / len (loader) def evaluate ( model, loader, device ): model. eval () correct, total = 0 , 0 with torch.no_grad(): for x, y in loader: x, y = x.to(device), y.to(device) outputs = model(x) correct += (outputs.argmax(dim= 1 ) == y). sum ().item() total += y.size( 0 ) return correct / total * 100 train_one_epoch 就是标准的四步循环:清零梯度 → 前向 → 算损失 → 反向 + 更新。开头的 model.train() 把模型切到训练模式,对本文这个模型其实没有影响(它没有 Dropout 和 BatchNorm),但这是必须养成的习惯——哪天加了 Dropout,不写这句就等于没生效。
损失累加那里有个小细节: total_loss += loss.item() 累加的是每个 batch 的 平均 损失,最后除以 len(loader) (batch 个数)。严格来说,只有当每个 batch 大小相同时,这个值才等于整个 epoch 的样本平均损失。本文训练集 2400 条、 batch_size=64 ,最后一个 batch 只有 32 条,所以算出来的损失有极轻微的偏差。影响可以忽略,但知道这个口径有助于以后排查"loss 看着不对劲"的问题;要精确的话可以累加 loss.item() * x.size(0) ,最后再除以样本总数。
evaluate 里有几个值得指出的点:
model.eval() 和 torch.no_grad() 是两件事。前者关闭训练专属行为,后者不建计算图,省显存也更快。评估时忘了 no_grad 不会算错,只是白白多占一份显存。 outputs.argmax(dim=1) 在 2 维的类别分数上取较大者,得到预测类别; (outputs.argmax(dim=1) == y) 得到的是布尔张量, .sum().item() 会先把它转成 0/1 再求和——这是 PyTorch 里统计正确率最常用的写法。 运行 def main (): device = torch.device( "cuda" if torch.cuda.is_available() else "cpu" ) print ( f"使用设备: {device} " ) X, y = generate_data( 3000 ) train_loader = DataLoader(TensorDataset(X[: 2400 ], y[: 2400 ]), batch_size= 64 , shuffle= True ) test_loader = DataLoader(TensorDataset(X[ 2400 :], y[ 2400 :]), batch_size= 128 ) model = LSTMClassifier().to(device) criterion = nn.CrossEntropyLoss() optimizer = optim.Adam(model.parameters(), lr= 1e-3 ) for epoch in range ( 1 , 21 ): loss = train_one_epoch(model, train_loader, optimizer, criterion, device) acc = evaluate(model, test_loader, device) if epoch% 5 == 0 or epoch== 1 : print ( f"Epoch {epoch:2d} /20 | 损失: {loss: .4 f} | 准确率: {acc: .2 f} %" ) print ( f"\n完成!准确率 {acc: .2 f} %(预期 95%+)" ) print ( "→ 结论:局部模式检测用 CNN;长程依赖/变长序列才需要 RNN/Transformer" ) 几处说明:
数据划分是 8:2 (2400 训练 / 600 测试),没有单独切验证集,是直接拿测试集来观察准确率的。严格说,用测试集反复挑 epoch、调参数,就等于把测试集的信息泄漏进了模型选择,报出来的数字会偏乐观。这个任务简单、也没有做 early stopping,所以影响不大;但正式实验里应该从训练集再切 10% 出来当验证集。 shuffle=True 只给训练集开 :训练时打乱能破坏样本顺序带来的相关性;测试集不需要打乱,评估结果和顺序无关。 只在第 1 轮和每 5 轮打印 ,是为了日志清爽。想看收敛细节,把 epoch % 5 == 0 这个条件去掉即可。 运行 main() 的输出是:
使用设备: cuda Epoch 1 / 20 | 损失: 0 .6841 | 准确率: 66 .50 % Epoch 5 / 20 | 损失: 0 .5291 | 准确率: 73 .00 % Epoch 10 / 20 | 损失: 0 .2018 | 准确率: 94 .50 % Epoch 15 / 20 | 损失: 0 .1182 | 准确率: 95 .67 % Epoch 20 / 20 | 损失: 0 .0951 | 准确率: 97 .17 % 完成!准确率 97.17 %(预期 95 %+) → 结论:局部模式检测用 CNN;长程依赖/变长序列才需要 RNN/Transformer 结果解读 把日志和 CNN 那一篇放在一起看,有几件事值得注意。
一、模型其实还没训完。 损失从 0.0951 还在往下走,准确率 94.50 → 95.67 → 97.17 一路在涨,两条曲线都没有走平的迹象。对比 CNN 版:第 5 个 epoch 就到 98.88%,第 7 个 epoch 损失降到 0.0052,第 10 个 epoch 99.88%—— CNN 早就收敛了,LSTM 到第 20 轮还差得远 。这正是 RNN 的固有短板:序列必须一个时间步一个时间步顺序算,没法并行;梯度又要沿 100 个时间步连乘回传(BPTT),收敛自然慢。所以 97.17% 是 被 epoch 数限制住的下限,不是 LSTM 能力的上限 。要公平比较,得先把 epoch 加到 50~100,或者用早停。
二、这个任务其实有个"作弊解"。 只要判断序列里有没有 TATAAA 或 TATACA ,就能拿到接近 100% 的准确率(负样本已经被刻意剔干净了)。CNN 的 99.88% 基本摸到了这个上界;LSTM 的 97.17% 说明它还没学会这条规则,漏掉的那部分样本,多半是 motif 落在序列两端、或者被平均池化稀释掉的那些。
三、池化方式的影响,可能比"用 CNN 还是 RNN"更大。 CNN 那篇最后用的是 MaxPool1d ,等价于"哪个位置反应最强就取哪个";这里的 LSTM 用的是 out.mean(dim=1) ,是"把所有位置平均"。对"序列里有没有出现某个局部模式"这类问题, 最大池化天然更对路 :motif 的 6 个位置被 94 个背景位置一平均,信噪比就掉下来了。
这是这篇笔记最值得动手验证的一步:把 out.mean(dim=1) 换成 out.max(dim=1).values ,其他都不动,再跑一遍看看准确率能到多少。按上面的分析,应该会明显高于 97.17%。
四、两个 motif 变体确实让任务变难了。 CNN 篇只有一个 TATAAA ,正样本里全是完全匹配;这里多了 TATACA ,模型必须容忍一处错配。不过这是 数据难度 的差异,不能直接拿 99.88% 和 97.17% 去说"CNN 比 LSTM 强"——要比就得在同一份数据、同样的 epoch 数下比。
小结:CNN vs LSTM 做同一件事 索引 (N, L) + Embedding,像一个句子
MaxPool1d mean over L
本文里 LSTM 只跑了 20 个 epoch、用的是平均池化、还是单向的,这三点都还没优化,所以"97% vs 99.88%"只能说明 CNN 在这个任务上更省事 ,不能说明 LSTM 学不会。真想把 LSTM 的潜力压出来,至少可以做三件事:① 池化换成 max;② 换成双向 LSTM;③ epoch 加到 50 以上或加早停。
参考资料 Hochreiter, S., & Schmidhuber, J. (1997). Long Short-Term Memory. Neural Computation, 9(8), 1735–1780. —— LSTM 的原始论文 Alipanahi, B., Delong, A., Weirauch, M. T., & Frey, B. J. (2015). Predicting the sequence specificities of DNA- and RNA-binding proteins by deep learning. Nature Biotechnology, 33(8), 831–838. —— DeepBind,用 CNN 做 motif / 结合位点预测的代表工作,也是"这类任务上 CNN 是主流"的依据 PyTorch 文档: nn.LSTM https://pytorch.org/docs/stable/generated/torch.nn.LSTM.html 、 nn.Embedding https://pytorch.org/docs/stable/generated/torch.nn.Embedding.html