Py学习  »  机器学习算法

深度学习连载10:用RNN识别DNA启动子

Y叔 余光创 • 1 月前 • 305 次点击  

上次我们用CNN重做了一遍RNN对电影评论分类,这次我们用RNN来重做一遍CNN识别DNA启动子。

这两次"交换作业"是配套的:在《使用CNN对电影评论进行分类》里,CNN抢了RNN的活;这一篇反过来,让LSTM去干CNN最擅长的活。只有把同一个任务交给两个模型各跑一遍,才能看清它们各自的适用边界在哪里。

  • CNN 用卷积核在序列上滑动,天然擅长局部、位置无关的模式。TATA box 正好是一个长度 6bp 的固定串,出现在哪个位置都算——这几乎是给卷积量身定做的任务。
  • LSTM 用循环单元一个碱基一个碱基地读过去,靠隐藏状态"记住"前面看过什么,它的强项是长程依赖和变长序列(见《循环神经网络RNN》)。

所以这一篇的预期结果不是"RNN 更好",而是"RNN 也能做,但要付出更多代价"。具体差多少,文末有对照。

  • 1
  • 2
  • 3
  • 4
  • 5
import randomimport torchimport torch.nn as nnimport torch.optim as optimfrom torch.utils.data import DataLoader, TensorDataset

数据

还是一样地产生随机序列,这次motif加多一个,有两个变体,这样更真实一些。

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
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向量。

  • 1
  • 2
  • 3
  • 4
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 对得上就行。)

模型

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
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__ 里三层各自的职责:

  1. nn.Embedding(4, 32):词表大小 4,嵌入维度 32。权重矩阵形状 (4, 32),一共只有 128 个参数,是个非常小的查表。
  2. nn.LSTM(32, 64, batch_first=True):input_size=32 接 Embedding 的输出,hidden_size=64。batch_first=True 是个高频坑——它表示输入形状是 (N, L, 32)(批次在前),而不是 PyTorch 默认的 (L, N, 32)。忘了写这个参数代码不会报错,但会把"批次"当成"序列长度",训练结果完全错乱。
  3. nn.Linear(64, 2):把 64 维的序列表示映射成 2 个类别的分数(logits),交给 CrossEntropyLoss。

forward 里的形状变化,一路写清楚:

步骤
代码
输出形状
含义
输入
x(N, 100)
碱基索引
嵌入
self.embedding(x)(N, 100, 32)
每个碱基变成一个 32 维向量
循环
self.lstm(x)
 的 out
(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。

训练与评估

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
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 里统计正确率最常用的写法。

运行

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
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)   ifepoch%5==0orepoch== 1:            print(f"Epoch {epoch:2d}/20 | 损失: {loss:.4f} | 准确率: {acc:.2f}%")
    print(f"\n完成!准确率 {acc:.2f}%(预期 95%+)")    print("→ 结论:局部模式检测用 CNN;长程依赖/变长序列才需要 RNN/Transformer")

几处说明:

  • 数据划分是 8:2(2400 训练 / 600 测试),没有单独切验证集,是直接拿测试集来观察准确率的。严格说,用测试集反复挑 epoch、调参数,就等于把测试集的信息泄漏进了模型选择,报出来的数字会偏乐观。这个任务简单、也没有做 early stopping,所以影响不大;但正式实验里应该从训练集再切 10% 出来当验证集。
  • shuffle=True 只给训练集开:训练时打乱能破坏样本顺序带来的相关性;测试集不需要打乱,评估结果和顺序无关。
  • 只在第 1 轮和每 5 轮打印,是为了日志清爽。想看收敛细节,把 epoch % 5 == 0 这个条件去掉即可。

运行main()的输出是:

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
使用设备: cudaEpoch  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 做同一件事


用 CNN(05 篇)
用 RNN/LSTM(本篇)
输入编码
独热 (N, 4, L),像一张图
索引 (N, L) + Embedding,像一个句子
核心算子
卷积核滑窗(局部、权重共享)
循环单元(顺序、隐藏状态)
池化
MaxPool1d
(取最强响应)
mean over L
(取平均)
数据
4000 条,1 个 motif
3000 条,2 个 motif 变体
epoch
10
20
准确率
99.88%(第 7 轮就收敛)
97.17%(仍在上升)
能否并行
能
不能,必须顺序
适合
局部、位置无关的模式
长程依赖、变长序列

本文里 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

Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/201360