AI写基因组：用语言模型设计新型噬菌体

如果 AI 不仅能写代码、写文章，还能直接写出一段几十 kb 的病毒基因组，并且这段 DNA 被合成之后，在实验室里真的能组装成具有感染能力的噬菌体颗粒——你会怎么评估这件事的分量？

这正是斯坦福大学等机构发表在 Science 上的工作给生物技术圈带来的冲击。研究团队把过去五年在自然语言处理领域被反复验证的“语言模型”路线，直接搬到了基因组层面：先在大量噬菌体基因组上进行预训练，让模型学会 DNA 序列的“语法”，再让模型生成全新的噬菌体基因组，并通过合成生物学手段去验证这些序列是否真的成立。项目标题里那句“探索自然进化之外的生物设计”，点出了这件事最核心的野心：不再只靠筛选自然界的现成病毒，而是尝试用 AI 直接从序列空间里“写”出新的生命形态。

先给出我对这件事的判断：这项工作的意义不在于“AI 又生成了一段 DNA”这种新闻式表述，而在于它把生物设计从“理解规则后再改造”推进到了“从数据中隐式学习规则并直接生成”。对做 AI 的开发者来说，这是一个典型的大模型垂直落地样本；对做生物信息、合成生物学的工程师来说，它意味着一条新的“设计-合成-验证”闭环开始成型。

这篇文章会用尽量直白的语言拆解五件事：基因组语言模型到底学的是什么，为什么偏偏选噬菌体当试验场，论文的生成与验证流程可以怎么理解，你自己如何用开源工具跑通一个最小示例，以及真正落地时会遇到哪些工程问题、生物安全边界和下一步值得关注的方向。即使你不做生物信息，只看 AI 应用层，这篇文章也能帮你判断：当一个领域开始用“语言模型”重写原有流程时，真正的机会和坑分别在哪里。

1. 这篇文章真正要解决的问题

先说一个现实背景。抗生素耐药性已经成为全球公共卫生领域的头号难题，而噬菌体（能够特异性感染细菌的病毒）被认为是后抗生素时代最有潜力的替代方案之一。它不感染人类细胞，特异性高，而且在自然界里数量极其庞大。既然噬菌体这么好，为什么没有大规模进入临床和工业应用？核心障碍在于：传统的噬菌体筛选和改造太慢、太不可控。

传统路径有两种。一种是“筛”：从环境样本里分离野生噬菌体，测试它对目标病原菌的裂解能力。这个流程高度依赖运气，筛到的噬菌体往往裂解谱窄、携带不明功能基因，甚至带有溶原性风险。另一种是“改”：基于已知的受体结合蛋白、裂解机制等分子生物学知识，对噬菌体进行理性改造。问题在于，噬菌体的感染机制涉及多个基因和调控元件的协同，知识驱动的改造往往会顾此失彼，开发周期以年为单位。

基因组语言模型要解决的，正是这个“序列-功能”复杂性带来的设计瓶颈。它的思路和传统方法完全不同：不依赖人类专家把每一个规则都写清楚，而是让模型在大规模基因组数据里自己学习哪些序列模式是保守的、哪些区域之间存在相互作用，然后以生成的方式给出候选设计。

所以，谁最应该读这篇文章？

• 做 AI 算法和工程的人：你需要理解大模型如何从文本领域迁移到基因组领域，tokenization、模型结构、生成策略在哪里发生了改变。
• 做生物信息分析和合成生物学的人：你需要理解一条新的设计闭环，以及如何评估 AI 生成的基因组是否靠谱。
• 做药物研发和微生物检测产品的人：你需要判断这项技术什么时候能从论文走向管线。
• 纯技术爱好者：这是一个理解“大模型不是只能处理人类语言”的最佳案例。

2. 基础概念：基因组语言模型、噬菌体与“可设计的生物序列”

2.1 噬菌体为什么是一个好的设计对象

噬菌体本质上是一类专门感染细菌的病毒。它的基因组大小通常在 5 kb 到 200 kb 之间，相比动辄上百万甚至上亿碱基的细菌和真核基因组，要小好几个数量级。正因为基因组小，它成为合成生物学里最容易被“整段合成、整段测试”的生命对象。

现代 DNA 合成技术可以按设计好的碱基序列直接合成 DNA 片段，因此只要你有了一段完整的噬菌体基因组序列，理论上就能把它合成出来，导入宿主细菌后让它重新组装成有活性的病毒颗粒。这就是基因组语言模型的输出能“闭环验证”的关键前提。你在模型侧生成序列的成本很低，但在湿实验侧验证一条序列是否成立的成本依然很高，所以噬菌体这种“小基因组、可合成、表型清晰”的对象，就成了最理想的试验场。

2.2 从语言模型到基因组语言模型

我们熟悉的 ChatGPT 类模型，处理的对象是人类语言，基本单位是“词”或“子词”。模型通过预测下一个词的方式，从海量文本中学习语言的内在规律。基因组语言模型做的事情几乎一样，只是把对象换成了 DNA 序列。

DNA 可以看作由 A、T、G、C 四种碱基组成的字符串。一个噬菌体基因组，就是一条由数万个碱基组成的“句子”。模型读入大量这样的“句子”后，会学习到一系列规律：哪些位置容易出现编码蛋白的开放阅读框，哪些位置是启动子，哪些碱基组合在进化压力下被强烈保守，哪些区域可以有更高的突变自由度。

需要强调的一点是：模型并没有被显式告诉“这里是基因”或“这里是启动子”，它是通过序列上下文自己推断出这些统计规律的。这是基因组语言模型与基于规则的生物信息工具的本质区别——前者学的是隐式“语法”，后者用的是显式“规则”。

2.3 与蛋白质语言模型和传统生物信息方法的区别

很多人会把基因组语言模型和蛋白质语言模型混为一谈，但二者在输入单位、建模范围和输出形态上都有明显差异。

蛋白质语言模型擅长在“单个蛋白质”的尺度上理解序列-结构-功能关系，而基因组语言模型的野心更大：它试图理解整条染色体或整个基因组内部的“上下文关系”。比如一个基因的表达强弱可能取决于它周围几 kb 的调控序列，一个噬菌体的感染能力可能由多个模块协同决定。这些跨区域、长距离依赖，恰恰是基因组语言模型想要捕捉的东西。

3. 为什么选噬菌体：最小的“可验证、可闭环”基因组对象

如果把基因组设计比作写文章，那么设计一个细菌基因组相当于写一部长篇小说，设计一个噬菌体基因组则相当于写一篇短篇散文。篇幅越小，模型的出错率就越可控，湿实验验证的周期也越短。

噬菌体基因组有三个非常适合作为生成任务验证对象的特征。

第一个特征是模块化。噬菌体基因组通常存在清晰的模块结构：吸附模块负责识别并结合宿主表面的受体，注入模块负责把 DNA 注射进细菌，复制模块负责劫持宿主机器合成子代，裂解模块负责在组装完成后裂解宿主细胞释放新颗粒。这种模块化结构意味着模型生成的序列如果出问题，研究者可以通过比对模块来判断问题出在哪一段。

第二个特征是表型可观察。把合成的噬菌体基因组转入大肠杆菌等宿主后，如果序列本身是完备的、功能是正常的，培养皿上会出现清晰的噬菌斑。这个表型读数是二进制式的：有活性还是没有活性，不需要复杂的生化检测就能判断。

第三个特征是进化压力明确。自然界中的噬菌体经过了漫长的宿主-寄生军备竞赛，序列里保留了大量与功能和适应性相关的信号。对模型来说，这些信号就是训练时最强的“监督信息”；对研究者来说，生成的序列是否“像真的”，也可以通过序列统计特征快速判断。

换句话说，噬菌体之于基因组设计，有点像 MNIST 之于深度学习：规模小、语义清楚、验证便宜。但它又不是一个简单玩具，因为即使是一个 50 kb 的基因组，背后也包含了几十个基因、多个调控元件和一套完整的生命周期逻辑。能在这个小基因组上跑通“生成-合成-验证”闭环，才有底气去挑战更大的基因组。

4. 核心方法拆解：一个基因组语言模型的训练与生成流程

从公开信息和这类工作的一般流程来看，用基因组语言模型生成新型噬菌体，并不只是“把序列丢进 GPT 再采样”这么简单。整个流程可以拆成五步，每一步都有自己的技术难点。

4.1 数据准备：把噬菌体基因组整理成“语料库”

语言模型的质量上限由语料决定，基因组语言模型同理。训练数据的核心来源是公共数据库中的噬菌体基因组序列，包括完整的基因组、注释信息以及宿主信息。数据准备阶段要做的事情包括：去冗余、过滤低质量序列、统一正负链方向、剔除混入的质粒或宿主 DNA 污染。

这一步决定了模型能学到什么。如果数据里混入了大量非噬菌体序列，模型生成的“噬菌体”可能实际更像质粒；如果只用了少数几个科的噬菌体，模型就无法覆盖噬菌体基因组的多样性。

4.2 Tokenization：基因组如何变成 token 序列

自然语言处理中的 tokenizer 负责把文本切成模型能处理的 token。DNA 序列只有四个字母，看起来很简单，但 tokenization 策略并不 trivial。

最简单的做法是单碱基 token，也就是把 A、T、G、C 分别映射为四个 token。但这种做法会让模型看到的每个位置信息量太少，很难捕捉密码子级别的偏好。常用的替代方案是 k-mer tokenization，即把序列切成长度为 k 的片段，k 通常取 3 到 6。k=3 时正好对应一个密码子，能帮模型捕捉编码区统计规律；k=6 时能覆盖部分短 motif，但词表会膨胀到 4 的六次方即 4096，尚在可控范围。

更高级的方案是直接在 DNA 序列上做类似 BPE（Byte Pair Encoding）的子词切分，让模型自己学习哪些片段应该合并成一个 token。不同切分方式会影响模型能学到的规律，这是基因组语言模型工程里最容易踩坑、也最容易被忽视的环节。

4.3 模型结构：生成式预训练模型如何学习基因组语法

从技术选型看，这类工作通常采用生成式预训练 Transformer，训练目标就是经典的自回归下一个 token 预测：给定前面一段序列，预测下一个 token。这个目标不需要人工标注，正好匹配基因组数据“量大但标注稀少的”特点。

模型在训练中会隐式学到多种规律：密码子偏好、基因边界信号、启动子和终止子的序列模式、不同功能模块之间的共现关系。随着模型规模增大和训练数据增多，它能捕捉的依赖距离也会变长，从局部密码子偏好逐步扩展到跨基因的关联。

4.4 生成策略：怎么从模型里采样“像噬菌体又不一样”的序列

训练完成后，生成策略决定了产出序列的性质。如果直接贪心解码，模型大概率会输出训练集中常见序列的“平均态”，缺乏多样性。如果完全随机采样，又会产生大量不符合生物学约束的序列。

实际工作中通常会在采样时引入温度参数和 top-k/top-p 截断，控制生成序列的多样性与质量之间的平衡。更进一步的策略是约束生成，比如在采样时强制保持 GC 含量在合理区间，或者把某些关键功能域序列作为锚定片段，让模型围绕锚定片段生成周边序列。从论文报道的性质看，“生成新型噬菌体”并不是让模型随机乱写，而是在多样性、新颖性和功能可行性之间做精细调衡。

4.5 湿实验验证：生成序列如何被合成并测试

生成序列只是起点。论文的完整证据链必须包含湿实验验证：将模型生成的基因组序列通过 DNA 合成手段拼接起来，转入宿主细菌，观察是否能够产生具有感染活性的噬菌体颗粒。这一步才是“自然进化之外的生物设计”真正成立的关键。

从工程视角看，整个流程是一个典型的闭环：模型生成候选序列，计算层过滤掉明显不合理的候选，再进入湿实验验证，验证结果反过来指导模型迭代。这个闭环的高效运转，才是基因组语言模型真正的价值所在。

5. 环境准备与前置条件

下面进入可操作环节。需要先说明：我不可能在这里复现 Science 论文的完整实验，那需要大规模 GPU 集群和完整的合成生物学实验室。我会带你跑通一个“最小教学版”，目标是让你理解基因组语言模型的核心原理，并能够处理真实的噬菌体基因组数据。

推荐环境如下：

• 操作系统：Linux（Ubuntu 20.04 或 22.04 最佳），macOS 也基本兼容。
• Python 版本：3.9 以上。
• 基础依赖：BioPython（用于读取 FASTA 基因组文件）、PyTorch（用于模型训练和生成）、pandas（用于数据处理）。
• 硬件要求：教学示例只需要 CPU 即可跑通；如果想训练更大的模型，建议使用 16 GB 以上显存的 GPU。

安装命令如下：

# 创建虚拟环境（推荐） python3 -m venv phage_env source phage_env/bin/activate # 安装基础依赖 pip install --upgrade pip pip install biopython pandas torch

# 验证安装 python3 -c "import Bio; print('BioPython', Bio.__version__)" python3 -c "import torch; print('PyTorch', torch.__version__)"

需要说明的是，真实论文使用的数据规模、模型参数规模和训练算力都会远超教学示例。这里的代码价值在于帮助你建立从“原始基因组数据”到“模型训练”再到“序列生成”的完整技术直觉，而不是替代生产环境。

6. 从零到一：一个可运行的基因组语言模型最小示例

这一节的代码全部是教学级实现，我会给出完整的文件路径和说明。你可以照着把流程跑通，然后用自己的数据做替换实验。

6.1 读取噬菌体基因组 FASTA

第一步是读取基因组序列。我们使用 BioPython 解析 FASTA 文件，并做最小预处理。建立一个examples目录，把下面的代码保存为examples/phage_loader.py。

# 文件路径：examples/phage_loader.py from pathlib import Path from Bio import SeqIO def load_genomes(fasta_path, min_len=5000): """ 读取 FASTA 文件中的基因组序列。 返回一个列表，每个元素是 {'id': 序列ID, 'seq': 大写DNA序列}。 """ sequences = [] for record in SeqIO.parse(fasta_path, "fasta"): seq = str(record.seq).upper().replace("U", "T") if len(seq) >= min_len: sequences.append({"id": record.id, "seq": seq}) return sequences def load_genomes_from_dir(fasta_dir, min_len=5000): """批量读取目录下的多个 FASTA 文件。""" fasta_dir = Path(fasta_dir) patterns = ["*.fna", "*.fasta", "*.fa"] result = [] for pattern in patterns: for fasta_file in sorted(fasta_dir.glob(pattern)): result.extend(load_genomes(fasta_file, min_len=min_len)) return result

这段代码的作用很直接：把基因组文件解析成 Python 字典列表，同时做两个关键预处理——统一大写、把 RNA 序列中的 U 替换为 T。min_len参数可以帮助过滤掉碎片化的短序列。

6.2 用 k-mer 构造 token 序列

第二步实现一个简化的 k-mer tokenizer。注意，这只是教学用的极简版本，真实基因组语言模型的 tokenizer 要复杂得多。保存为examples/phage_tokenizer.py。

# 文件路径：examples/phage_tokenizer.py from collections import Counter class KMerTokenizer: """ 简化的 k-mer tokenizer。 把 DNA 序列切分为长度为 k 的片段，并为每个片段分配一个整数 ID。 """ def __init__(self, k=6): self.k = k self.vocab = {} self._id2kmer = {} def fit(self, sequences): """根据训练序列构建词表。""" counter = Counter() for seq in sequences: seq = seq.upper().replace("U", "T") if len(seq) < self.k: continue for i in range(len(seq) - self.k + 1): counter[seq[i:i + self.k]] += 1 # 0 号 ID 留给未知 token self.vocab = {"<UNK>": 0} for idx, (kmer, _) in enumerate(counter.items(), start=1): self.vocab[kmer] = idx self._id2kmer = {idx: kmer for kmer, idx in self.vocab.items()} return self def encode(self, seq): """DNA 序列 -> token ID 列表。""" seq = seq.upper().replace("U", "T") tokens = [] for i in range(0, len(seq) - self.k + 1, self.k): kmer = seq[i:i + self.k] tokens.append(self.vocab.get(kmer, 0)) return tokens def decode(self, tokens): """token ID 列表 -> DNA 序列（仅用于教学演示）。""" return "".join(self._id2kmer.get(t, "N") for t in tokens) @property def vocab_size(self): return len(self.vocab)

你可以看到，这个 tokenizer 的核心逻辑就是滑动窗口切 k-mer，然后查表转 ID。decode方法只是把 token 对应的 k-mer 拼接起来，用来观察生成结果；由于我们切分时步长等于 k，理论上可以无缝拼接回去。

6.3 训练一个最小 GPT 模型

第三步是模型本身。这里我用 PyTorch 实现一个带因果掩码的最小 Transformer，命名叫MiniPhageGPT，方便你理解“生成式预训练”的含义。保存为examples/mini_phage_gpt.py。

# 文件路径：examples/mini_phage_gpt.py import torch import torch.nn as nn import torch.nn.functional as F class MiniPhageGPT(nn.Module): """ 极简基因组语言模型。 只用来说明核心原理，不追求性能。 """ def __init__(self, vocab_size, d_model=128, nhead=4, num_layers=2, max_len=256, dropout=0.1): super().__init__() self.token_embedding = nn.Embedding(vocab_size, d_model) self.position_embedding = nn.Embedding(max_len, d_model) self.max_len = max_len encoder_layer = nn.TransformerEncoderLayer( d_model=d_model, nhead=nhead, dim_feedforward=d_model * 4, dropout=dropout, batch_first=True, ) self.blocks = nn.TransformerEncoder(encoder_layer, num_layers=num_layers) self.lm_head = nn.Linear(d_model, vocab_size) def forward(self, x): """ x: (batch, seq_len) 的 token ID 返回: (batch, seq_len, vocab_size) 的 logits """ batch, seq_len = x.shape positions = torch.arange(seq_len, device=x.device).unsqueeze(0) hidden = self.token_embedding(x) + self.position_embedding(positions) # 因果掩码：上三角设为 -inf，使每个位置只能看到当前位置及之前 causal_mask = torch.triu( torch.full((seq_len, seq_len), float("-inf"), device=x.device), diagonal=1, ) hidden = self.blocks(hidden, mask=causal_mask) logits = self.lm_head(hidden) return logits @torch.no_grad() def generate(self, prompt_tokens, max_new_tokens=100, temperature=0.8, top_k=50): """ 自回归生成：根据 prompt 逐个预测下一个 token。 """ self.eval() tokens = list(prompt_tokens) for _ in range(max_new_tokens): input_ids = torch.tensor( [tokens[-self.max_len:]], dtype=torch.long ) logits = self(input_ids)[0, -1, :] / temperature if top_k is not None: top_values, _ = torch.topk(logits, top_k) logits[logits < top_values[-1]] = float("-inf") probs = F.softmax(logits, dim=-1) next_token = torch.multinomial(probs, num_samples=1).item() tokens.append(next_token) return tokens

代码里最关键的是forward中的因果掩码。TransformerEncoder 默认会看到整个序列，但生成式模型要求当前 token 只能看到它之前的 token。我把上三角矩阵设为-inf，这样自注意力在计算时，未来位置的信息会被完全屏蔽。generate方法实现了标准的自回归生成：把当前 token 序列喂给模型，取最后一个位置的 logits，经过温度和 top-k 控制后采样下一个 token。

6.4 训练与生成主流程

第四步把前面的模块串起来。保存为examples/train_mini_phage_gpt.py，用法是传入一个 FASTA 文件和训练轮数。

# 文件路径：examples/train_mini_phage_gpt.py import argparse import torch import torch.nn.functional as F from torch.utils.data import DataLoader, Dataset from phage_loader import load_genomes from phage_tokenizer import KMerTokenizer from mini_phage_gpt import MiniPhageGPT class SeqDataset(Dataset): """把 token 序列切成固定长度的训练样本。""" def __init__(self, tokenized_seqs, block_size=128): self.samples = [] for tokens in tokenized_seqs: if len(tokens) < block_size: continue for start in range(0, len(tokens) - block_size, block_size // 2): chunk = tokens[start:start + block_size] self.samples.append(chunk) def __len__(self): return len(self.samples) def __getitem__(self, idx): chunk = self.samples[idx] x = torch.tensor(chunk[:-1], dtype=torch.long) y = torch.tensor(chunk[1:], dtype=torch.long) return x, y def collate_fn(batch): """按 batch 内最大长度做 padding，target 部分用 -100 屏蔽。""" max_len = max(x.size(0) for x, _ in batch) xs, ys = [], [] for x, y in batch: pad_len = max_len - x.size(0) xs.append(torch.cat([x, torch.zeros(pad_len, dtype=torch.long)])) ys.append(torch.cat([y, torch.full((pad_len,), -100, dtype=torch.long)])) return torch.stack(xs), torch.stack(ys) def train(model, tokenizer, sequences, epochs=10, batch_size=8, lr=3e-4): tokenized = [tokenizer.encode(s) for s in sequences] dataset = SeqDataset(tokenized, block_size=128) loader = DataLoader(dataset, batch_size=batch_size, shuffle=True, collate_fn=collate_fn) optimizer = torch.optim.AdamW(model.parameters(), lr=lr) model.train() for epoch in range(epochs): total_loss = 0.0 for x, y in loader: logits = model(x) loss = F.cross_entropy( logits.view(-1, logits.size(-1)), y.view(-1), ignore_index=-100, ) optimizer.zero_grad() loss.backward() torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0) optimizer.step() total_loss += loss.item() print(f"epoch {epoch + 1}/{epochs}, loss={total_loss / len(loader):.4f}") def main(): parser = argparse.ArgumentParser(description="Mini Phage GPT 教学示例") parser.add_argument("--fasta", required=True, help="噬菌体基因组 FASTA 文件") parser.add_argument("--epochs", type=int, default=10) parser.add_argument("--k", type=int, default=6) args = parser.parse_args() genomes = load_genomes(args.fasta, min_len=5000) if not genomes: print("未读取到足够长的基因组序列，请检查 FASTA 文件。") return sequences = [g["seq"] for g in genomes] print(f"读取到 {len(sequences)} 条基因组序列") tokenizer = KMerTokenizer(k=args.k).fit(sequences) print(f"词表大小: {tokenizer.vocab_size}") model = MiniPhageGPT(vocab_size=tokenizer.vocab_size) train(model, tokenizer, sequences, epochs=args.epochs) # 用第一条序列的前 100 个 token 作为 prompt，生成后续序列 prompt_tokens = tokenizer.encode(sequences[0])[:100] generated_tokens = model.generate(prompt_tokens, max_new_tokens=200) generated_seq = tokenizer.decode(generated_tokens) print("生成序列（片段）:") print(generated_seq) if __name__ == "__main__": main()

这段代码把所有环节串成了一条线：读取 FASTA、拟合 tokenizer、创建数据集、训练模型、自回归生成。训练时用ignore_index=-100屏蔽 padding 位置，避免它们参与 loss 计算。生成时以一条真实序列的前 100 个 token 作为 prompt，让模型续写后面的内容。

6.5 运行与验证

假设你把噬菌体基因组 FASTA 文件放在data/phages.fasta，运行命令如下：

cd examples python train_mini_phage_gpt.py --fasta ../data/phages.fasta --epochs 10 --k 6

预期输出类似：

读取到 50 条基因组序列 词表大小: 4096 epoch 1/10, loss=5.9841 epoch 2/10, loss=5.6103 ... epoch 10/10, loss=4.8720 生成序列（片段）: ATGCGT...（一段模型生成的新序列）

如何判断训练是否成功？主要看两点：第一，loss 是否整体呈下降趋势；第二，生成序列是否仍然由合法的 k-mer 组成，而不是停在中途或出现大量未知 token。如果 loss 一直不降，优先检查数据量是否过少、学习率是否过高、block_size 是否设置合理。

7. 生成结果怎么看：in silico 与湿实验两类验证

教学示例跑通后，下一个核心问题是：怎么判断一段生成序列“像不像”一个真实的噬菌体基因组？论文工作使用了完整的湿实验验证，但我们在日常工程中会先做一系列计算层面的筛选，把明显不靠谱的候选过滤掉，再决定哪些值得进入合成环节。

7.1 序列层面的评价指标

下面给出一个简单的分析脚本，用几个基础统计量来比较生成序列与训练集序列的差异。保存为examples/analyze_phage_seq.py。

# 文件路径：examples/analyze_phage_seq.py from collections import Counter import argparse def gc_content(seq): """计算 GC 含量。""" seq = seq.upper() if not seq: return 0.0 return (seq.count("G") + seq.count("C")) / len(seq) def kmer_frequency(seq, k=3): """统计 k-mer 频率。""" seq = seq.upper() if len(seq) < k: return Counter() return Counter(seq[i:i + k] for i in range(len(seq) - k + 1)) def jaccard_kmers(counter_a, counter_b): """计算两个 k-mer 集合的 Jaccard 相似度。""" set_a, set_b = set(counter_a), set(counter_b) if not set_a or not set_b: return 0.0 return len(set_a & set_b) / len(set_a | set_b) def codon_density(seq): """粗略统计起始/终止密码子密度，仅供教学演示。""" seq = seq.upper() starts = sum( 1 for i in range(len(seq) - 2) if seq[i:i + 3] == "ATG" ) stops = 0 for codon in ("TAA", "TAG", "TGA"): stops += sum( 1 for i in range(len(seq) - 2) if seq[i:i + 3] == codon ) total_codons = max(len(seq) // 3, 1) return starts / total_codons, stops / total_codons def main(): parser = argparse.ArgumentParser(description="分析噬菌体序列基础统计") parser.add_argument("--seq", required=True, help="待分析的序列文件") parser.add_argument("--ref", help="参考序列文件，计算 Jaccard 相似度") args = parser.parse_args() with open(args.seq, "r") as f: seq = f.read().strip().upper() print(f"序列长度: {len(seq)}") print(f"GC 含量: {gc_content(seq):.4f}") start_density, stop_density = codon_density(seq) print(f"起始密码子密度: {start_density:.4f}") print(f"终止密码子密度: {stop_density:.4f}") seq_kmers = kmer_frequency(seq, k=3) if args.ref: with open(args.ref, "r") as f: ref_seq = f.read().strip().upper() ref_kmers = kmer_frequency(ref_seq, k=3) sim = jaccard_kmers(seq_kmers, ref_kmers) print(f"与参考序列的 3-mer Jaccard 相似度: {sim:.4f}") if __name__ == "__main__": main()

运行方式：

python examples/analyze_phage_seq.py --seq generated_seq.txt --ref natural_phage_seq.txt

这类指标的意义在于快速筛查。比如真实噬菌体基因组的 GC 含量通常有比较固定的范围，如果生成序列的 GC 含量偏离训练集统计分布太多，说明模型没有学好基本的碱基组成规律。k-mer 相似度则反映生成序列在短片段模式上与自然序列的重合程度，相似度太高说明模型只是在复读训练数据，太低则说明生成结果可能是无意义的随机串。

7.2 功能层面的验证

计算指标只是第一层筛选。真正能让一段生成序列“成立”的，是湿实验验证。从这类工作的一般流程看，