2. 基因组组装质量对注释的影响：N50、contig完整性、scaffold连接

做基因组注释这么多年，我踩过的坑一半以上都出在「组装质量」这关。说白了，上游组装要是豆腐渣，下游注释做得再花哨也是白搭。你想想看，基因预测依赖于序列的连续性和准确性——断裂的contig会切碎基因结构，错误的scaffold连接会把不同染色体片段拼在一起，导致预测出一堆根本不存在于真实基因组里的「嵌合基因」。

2.1 为什么组装质量会直接影响注释

注释的本质是在组装序列上寻找功能性特征。如果组装本身就有问题，那么：

• 断裂 → 基因被切成几段，预测软件要么漏掉，要么输出一堆残缺的外显子集。
• 错误连接 → 两条不同染色体的序列被拼到一起，结构基因两侧的调控元件全部错位。
• 冗余 → 重复区域被过度合并或错误拆分，导致基因拷贝数偏离真实情况。

我记得有个项目，用三代测序组装的基因组N50挺漂亮，但BUSCO完整性只有65%。一查发现，组装器把高重复区域全压成了N，结果注释出来的基因集缺少了一大堆保守的单拷贝基因。最后整个注释结果没法用，只能重新调参组装。

下面这张流程图清晰展示了组装质量指标与注释各环节之间的牵连关系：

2.2 N50：连续性指标如何影响基因预测

N50常被拿来吹嘘组装质量，但它和注释的关系可没那么简单。绝大多数基因的跨度在2～20 kb之间，如果contig N50低于基因中位数长度，那么超过一半的基因会被切成两段甚至更多段。

实战经验：我习惯在注释前先跑一个简单脚本，统计组装序列中长度大于等于基因中位数的contig数量占比。如果这个比例低于70%，就要警惕了——即使BUSCO分数还行，基因预测的完整性也会大打折扣。

来看一个真实的对比：

你想想看，B版本里每6个预测基因就有1个可能是假阳性——后续实验验证的成本谁受得了？

2.3 contig完整性：BUSCO 不是万能的，但少它不行

N50只告诉你长度，不告诉你内容。contig完整性指的就是组装序列中有多少「真实基因组」的基因区域被覆盖到。我们常用的工具是BUSCO，会用上千个保守单拷贝直系同源基因来评估。我个人认为，对于真核生物，BUSCO完整比例低于80%的组装，注释结果基本没办法直接发表。

我曾经接手过一个哺乳动物基因组项目，组装BUSCO完整度只有78%，但N50高达60 Mb。很多人会被N50迷惑，觉得「这么长肯定没问题」。结果一注释，发现大量看家基因只有部分外显子，有些甚至完全缺失。最后花了两个月手工补全长读长覆盖率才挽回局面。你说冤不冤？

评估命令其实很简单，用BUSCO跑一下：

busco -i genome.fasta -l eukaryota_odb10 -o busco_out -m genome

重点关注 Complete 和 Fragmented 两列。Fragmented比例高说明基因被截断，往往直接对应低质量注释。

2.4 scaffold连接：gap 和 chimeric 的陷阱

Scaffold的N50通常比contig N50大得多，但这里藏着最多的坑。scaffold连接依赖pair-end或Hi-C数据，如果连接错误（chimeric scaffold），就会把不同染色体上的序列粘在一起。注释软件拿到这样的序列，会预测出「跨染色体的融合基因」——这在真实基因组里根本不可能存在。

我记得有一次，某个真菌基因组scaffold N50达到整个染色体的80%，但注释出来一堆跨着丝粒的基因模型。一查发现，那个scaffold内部有大量N碱基填充区（gap），而gap两侧的contig其实不属于同一条染色体。这就是典型的「连接假象」。

怎么评估scaffold质量？我常用的组合拳：

1. 统计 gap数量 / 总scaffold长度 —— 如果每Mb超过100个gap，就要注意了。
2. 看 gap N50（即gap长度的中位数）—— 越长说明连接越模糊。
3. 比对到近缘参考基因组（如果有），检查scaffold边界处是否存在结构重排。

下面是我自己写的一个小脚本，用来快速评估scaffold连接质量：

#!/usr/bin/env python3
import sys
from collections import Counter

gap_count = 0
total_n = 0
scaf_lengths = []

for record in SeqIO.parse(sys.argv[1], "fasta"):
    seq = str(record.seq).upper()
    scaf_lengths.append(len(seq))
    gaps = seq.count('N' * 10)  # 长度≥10的gap串
    gap_count += gaps
    total_n += seq.count('N')

print(f"Scaffolds: {len(scaf_lengths)}")
print(f"Total gap count: {gap_count}")
print(f"Total N bases: {total_n}")
print(f"Gap density (per Mb): {total_n / (sum(scaf_lengths)/1e6):.1f}")

2.5 实践建议：组装质量审计三步走

开始注释前，我至少会花一个下午做质量审计。不是小题大做——后续所有工作都依赖这份审计结果。以下是我的三步流程：

• 第一步： 用QUAST跑一份报告，重点看 Largest contig、GC content 和 # mismatches per 100 kbp（如果有参考）。
• 第二步： 跑BUSCO，用与该类群最匹配的谱系数据库。真核生物用 eukaryota_odb10，昆虫用 insecta_odb10，以此类推。
• 第三步： 手动检查几个大scaffold的末端结构。我会选N50前后的scaffold，提取两侧各10 kb，与已知蛋白数据库比对。如果两侧的蛋白同源物属于不同物种或不同染色体，那就要重新审视整个组装策略了。

做完这三步，你心里就有底了。注释前明确知道哪些区域可靠、哪些区域需要特殊处理（比如只用证据较少的从头预测）。这样最终输出的注释集，投稿时才有底气。