做基因组组装的人大概都经历过这个场景contig N50 一路刷到了几十 Mb组装总长度跟流式估计的基因组大小也对得上染色体挂载也完成了满心以为这个基因组稳了。可等到做基因注释、做比较基因组学才发现一堆基因缺胳膊少腿某个基因家族莫名其妙少了好几个成员。这个时候再回头看组装N50 再漂亮也只能说明 contig 拼得长根本没法回答一个更关键的问题这个基因组把该有的基因完整保留下来了吗衡量这个问题BUSCO 基本是行业里默认的那把尺子。BUSCO 全称 Benchmarking Universal Single-Copy Orthologs通过一套覆盖不同谱系的单拷贝直系同源基因库去检查你的基因组、蛋白组或转录组里能不能找齐这些“标配”基因从而给出一个完整性百分比。现在基因组文章里如果缺了 BUSCO 结果评审基本直接拒稿。这篇文章我就把 BUSCO 的原理、实操、结果解读和那些最容易翻车的细节完整梳理一遍。1. 基因组完整性的需求拐点N50回答不了的问题1.1 传统指标为什么“看着够用”先聊个实际感受。组装指标里大家最常晒的就是 contig N50、scaffold N50、L50、组装总长、GC 含量这些。N50 翻译成人话就是把组装出的序列按长度从长到短排序累计长度达到总长一半时那条序列的长度。它衡量的是连续性——你拼出来的片段有多长可以直观理解为“这堆乐高拼得有多整”。但连续性和完整性是两个维度。我见过一个极端例子一条超长 contig 里其实插了一段寄生序列或者某条染色体臂整个丢失了但 N50 毫发无伤反过来一个碎片化严重、N50 只有几百 kb 的组装也可能在基因层面相当完整。这就好比衡量一栋房子N50 告诉你墙面有多完整、梁柱连得多顺但它不告诉你客厅和卧室是不是漏掉了一间。漏掉的那间房恰恰是对基因注释、比较基因组学、系统发育分析最伤的东西——因为后续分析默认假设你能拿到全套基因。组装总长度虽然可以和流式细胞术估算的基因组大小互相对照但同样有局限。基因组里有大量重复序列和高 AT 区域这些区域不容易组得全总长度偏小几百 Mb 也不奇怪。如果只盯着“长度够不够”很容易掩盖局部基因缺失的问题等往下游走才发现满盘皆输。1.2 BUSCO补的是哪块短板BUSCO 的思路完全绕开了连续性和总长度直接问你在进化上被广泛保留的基因你的组装里到底找齐了没有思路有点像装修验收时请第三方监理——你说墙面多平、板材多贵都没用我只抽查标准清单里的必检项目水电点位、防水高度、消防设施逐项核对。BUSCO 抽取的就是不同谱系里几乎每个物种都该有的单拷贝基因清单逐项核对你的基因组货架上有货没货。这也是为什么现在主流基因组文章都把 BUSCO 结果当作硬指标。除了期刊要求主流公共数据库提交基因组时也会要求提供组装质量评估信息BUSCO 就是最常见的一项。尤其是做了 scaffold 到 chromosome 级别提升的项目评审会非常在意挂载后有没有引入错误、有没有把不同单倍型的序列错误拼到一起一个靠谱的 BUSCO 结果能省掉很多来回补实验的麻烦。把 BUSCO 和传统指标放在一起看盲区差距一目了然指标衡量的是什么盲区contig/scaffold N50连续性不反映基因缺失组装总长度基因组大小近似重复区难组装掩盖局部缺失GC 含量碱基组成与完整性没有直接关系BUSCO C/D/F/M保守单拷贝基因完整度依赖谱系与数据库版本2. 为什么偏偏选中“单拷贝直系同源基因”当标尺2.1 单拷贝直系同源基因的筛选逻辑BUSCO 里每个字母都有讲究核心资产是 OrthoDB 数据库里那些符合“单拷贝、直系同源、谱系内广泛保守”的基因家族。先理解为什么强调单拷贝。如果一个基因在大多数物种里本身就存在多个拷贝比如某些抗病基因家族那么在评估基因组时就算组装里出现了 5 个拷贝你也很难判断这是天然如此、还是组装重复造成的假象。而选择在谱系里通常只有一个拷贝的基因一旦你的组装里出现两个及以上的 hits基本可以断定组装里存在冗余片段或未处理的单倍型。这种干净可判的性质让单拷贝基因成为天然质检探针。另一个容易被忽略的点BUSCO 要求这些基因在进化上足够保守但又不能过于保守。太保守的基因容易在重复区和假基因区域产生假阳性命中保守度适中的基因既能稳定存在又有足够的序列差异让比对软件能定位到正确的基因座。所以 OrthoDB 的筛选实际上是用大量参考基因组做交叉验证把那些在多物种中呈单拷贝、序列保守的基因家族挑出来再按谱系打包成数据库。你下载到的 arthropoda_odb10 这类文件本质上就是一张经过了重重筛选的“必备基因清单”。2.2 四种状态是怎么判出来的BUSCO 把每个基准基因在你输入文件中的状态归为四类CompleteC基因完整存在通常要求比对覆盖度和序列得分都超过阈值。DuplicatedD在 Complete 的基础上检测到 2 个及以上拷贝。它是 C 的子集。FragmentedF只找到了基因的一部分达不到完整长度的门槛。MissingM在你的输入里完全找不到对应序列。这里要特别说明在 genome 模式下这个判定不是简单拿序列去 BLAST 然后数数。BUSCO 会把输入基因组先做一轮基因结构预测——原核生物用 Prodigal真核生物会用 Augustus 或 MetaEuk 这类工具先把候选基因区域预测出来再把这些候选蛋白序列与对应的 BUSCO profile 用 HMMER 做比对。比对结果还要经过筛选、聚类、坐标合并最后才给出状态。proteins 模式和 transcriptome 模式则跳过了预测或使用不同流程所以速度快但也更受输入文件质量影响。一句话总结BUSCO 相当于拿一张清单去仓库盘点但它的目的不是把仓库所有货品都数一遍而是只抽查清单上那些“必备货号”。每个货号找到完整商品就记 C找到两个以上完整商品就记 D只找到一个零件就记 F压根没找到就记 M。3. 从零跑通一次BUSCO评估安装、选库、启动命令3.1 安装conda一条命令还是源码编译我自己的经验90% 的场景直接 conda 就能解决。建议新建一个独立环境避免跟其他生物信息软件互相污染依赖conda create -n busco_env -c conda-forge -c bioconda busco5.4.3 conda activate busco_env busco --version这里有一个版本细节BUSCO 5.x 和 4.x 的结果不能直接混比这一点在后面的踩坑章节我会重点展开。所以用 conda 创建环境时锁定版本号是最稳妥的比如busco5.4.3不要装成最新版然后跟旧结果对比那样得到的完整性变化可能就是版本噪声。没有 conda 或者服务器网络受限的话也可以直接用官方提供的 Docker 镜像把本地数据目录挂载进去运行。源码编译比较费事要处理 Python 依赖、HMMER、Augustus 路径等除非你要做二次开发否则我不建议。用 conda 装好的 BUSCO 会带上 HMMER、Biopython 等核心依赖省掉很多环境折腾。3.2 谱系数据库怎么选这是整个流程里最需要动脑的一步。数据库选得对不对直接决定结果的真实性和可比性。BUSCO 5.x 默认使用 ODB10 数据库不同谱系对应不同物种范围常见的谱系如下谱系名称适用物种范围举例bacteria_odb10细菌通用fungi_odb10酵母、霉菌、大型真菌eukaryota_odb10广谱真核生物arthropoda_odb10昆虫、甲壳类、蛛形纲vertebrata_odb10鱼类、两爬、鸟类、哺乳动物mammalia_odb10哺乳动物viridiplantae_odb10绿藻与陆生植物选库原则是“跟物种亲缘关系尽量匹配同时兼顾数据库质量”。比如你要评估一个蚂蚁基因组选 arthropoda_odb10 比选 eukaryota_odb10 更合适因为后者覆盖的谱系太宽挑出来的基因往往是极度保守的核心基因数量少、相对容易找全结果会偏乐观对蚁科特有的保守基因缺失反而发现不了。反过来如果手头是个偏门物种没有完全对应的子谱系那就选包含它的大类谱系不要硬套一个只覆盖近缘类群的小库。第一次跑不确定可以先用busco --list-datasets列出可用数据集或者在命令里直接写谱系名让 BUSCO 自动下载。不过我习惯先手动把数据库下到本地 datasets 目录这样第二次跑不会卡在下载和校验环节也能保证离线环境可用。3.3 三种运行模式与参数设置BUSCO 的输入可以是三类基因组序列、蛋白序列、转录本序列。三种模式对应三种评估需求genome 模式评估一个新组装的基因组。最常用需要跑基因结构预测最慢。proteins 模式评估已有注释的蛋白集。如果你已经做了基因注释手头有 proteins.fasta直接拿来测最省时间。transcriptome 模式评估 Trinity 之类的 RNA-seq 从头组装结果通常用于没有基因组时的临时评估。典型命令如下# genome 模式评估昆虫基因组 busco -i insect_genome.fasta \ -l arthropoda_odb10 \ -o insect_busco \ -m genome \ -c 16 # proteins 模式评估注释后的蛋白集 busco -i proteins.fasta \ -l arthropoda_odb10 \ -o insect_prot_busco \ -m proteins \ -c 16 # transcriptome 模式 busco -i transcripts.fasta \ -l arthropoda_odb10 \ -o insect_tx_busco \ -m transcriptome \ -c 16我习惯在需要覆盖已有输出目录时加--force但如果想保留之前的完整结果别加这个参数。另外-o指定的输出名不要跟输入文件名相同否则输出目录和中间文件会混在输入目录里清理时很头疼。关于基因预测器genome 模式下BUSCO 对真核生物会调用基因结构预测工具比较常见的是 Augustus 和 MetaEuk。conda 安装通常会把它们一起带上运行日志里会明确告诉你实际用了哪个预测器。如果你的真核基因组特别大比如超过 3 Gb建议关注一下预测这一步的耗时后续可以用--metaeuk切换到更快的 MetaEuk 流程速度提升明显。3.4 运行耗时与资源预估很多人第一次跑 BUSCO 最担心时间。真实情况差异很大细菌基因组几 Mb默认参数下几分钟到二十分钟。真菌基因组几十 Mb半小时到几小时。节肢动物或脊椎动物基因组几百 Mb 到几 Gb16 线程下通常 2 到 10 小时主要瓶颈在真核基因结构预测。proteins 模式通常分钟级因为省掉了预测步骤。内存一般不会特别恐怖但 genome 模式在比对和排序阶段会吃掉不少内存。我建议 16 线程跑 1 Gb 基因组时至少给 32 Gb 内存配额。线程不是越多越好BUSCO 有部分阶段是单线程的超过物理核心数后并不会线性加速。如果中途因节点故障中断BUSCO 支持断点续跑busco -i insect_genome.fasta -l arthropoda_odb10 -o insect_busco -m genome -c 16 --restart这个参数很实用。我之前跑一个脊椎动物基因组跑了七个多小时在最后阶段断掉用--restart十分钟就续完了。注意--restart和--force不能同时加续跑时输出目录必须保持原样。4. 结果文件不是念一行百分比就完事C、D、F、M的深读4.1 short_summary与full_table各是什么跑完之后最重要的结果是输出目录里的 short_summary 文件文件名类似short_summary.specific.arthropoda_odb10.insect_busco.txt。打开后核心内容如下# BUSCO version is: 5.4.3 # The lineage dataset is: arthropoda_odb10 (Creation date: 2021-02-19, number of BUSCOs: 1013) # Summarized benchmarking in BUSCO notation: # C: 96.2%[S: 91.8%, D: 4.4%], F: 1.5%, M: 2.3%, n: 1013翻译一下n1013arthropoda_odb10 这个库里有 1013 个用于评估的 BUSCO 基因。C96.2%96.2% 的基因在你的输入里是完整的。其中S91.8%是单拷贝完整D4.4%是重复完整D是C的子集所以C S D。F1.5%碎片化只找到部分片段。M2.3%完全缺失。校验规则C F M ≈ 100%。只看这一行综合百分比是不够的。full_table.tsv 里有每个 BUSCO 基因的具体状态、所在序列、起止坐标和得分这才是排查问题的真正入口。missing_busco_list.tsv和fragmented_busco_list.tsv分别列出缺失和破碎的清单。完整基因对应的序列会输出到busco_sequences/目录下按单拷贝、多拷贝、碎片、缺失分目录存放。这些文件在后续分析里非常有用——比如你想确认某个重要基因家族在目标基因组里是不是真的丢了可以直接提取对应的 BUSCO 序列做二次比对验证。4.2 duplicated比例偏高说明什么很多人看到 C 很高就放心了但我会额外盯一眼 D 的比例。D 超过 5% 就要警觉超过 10% 基本说明组装里有冗余。D 高的最常见原因不是物种天然多拷贝而是组装把同一段单倍型拼出了两个拷贝。比如杂合度高的物种如果组装前没有做单倍型分离、组装后没有做 purge_dups两条单倍型的差异区会被当成不同序列分别拼出来。BUSCO 会在两条序列上各找到一份“完整基因”于是报成 Duplicated。判断 D 是真实多拷贝还是组装冗余要看物种背景。对于大多数二倍体动物BUSCO 基因在单倍体视角下应该只有一份D 高基本就是冗余但如果是多倍体植物比如小麦、油菜天然就可能存在多套同源拷贝D 高不一定代表组装有问题需要进一步检查重复拷贝是否分布在预期的同源染色体上。处理思路也不复杂先用 purge_dups 这类工具去除冗余单倍型再重跑 BUSCO 对比。我见过一个项目D 从 16.8% 降到 2.7%同时 C 反而从 78.9% 升到 91.2%整体质量大幅改善。详细过程在第 5 章展开。4.3 missing与fragmented的成因排查M 偏高通常不是单一因素我归纳成四类谱系选错了比如用 vertebrata_odb10 去评估昆虫基因组一大批昆虫特有基因必然报 Missing但这不代表组装不好纯粹是尺子不匹配。组装本身缺了高复杂度区域着丝粒、rDNA 簇、高 AT 区经常组不全如果这些区域富集保守基因就会表现为 M 或 F。输入文件里有太多 N组装时用 N 连接 scaffold 很常见但如果 N 比例过高比如超过 5%会直接干扰下游基因结构预测导致大量 Fragmented。基因预测环节出了问题Augustus 训练不到位、依赖路径没配好都会让候基因质量差进而被误判为 Missing。F 偏高时我建议先把基因组做一遍软屏蔽重复序列再跑看看是否有改善。有些重复元件紧贴着基因上下游会干扰基因预测边界。快速判断可以参照下面表格现象最可能的几个原因建议动作D 高组装冗余、未处理单倍型检查杂合度、做冗余去除M 高谱系不匹配、组装确实缺失换谱系、检查 N 比例、重新评估F 高基因预测被干扰、序列碎片化软屏蔽重复序列、检查 N 比例四项正常组装比较健康可以进入注释等下游分析5. 实战复盘一个节肢动物基因组从C78%到C91%的过程5.1 原始组装评估问题出在冗余与碎片去年有个朋友做一种鞘翅目昆虫的基因组给我发来第一版组装的统计contig N50 有 4.8 Mb组装总长 820 Mb跟流式估算的 780 Mb 也比较接近。他问这版组装能不能继续往下走。我没直接回答先跑了一轮 BUSCO。busco -i beetle_v1.fasta \ -l arthropoda_odb10 \ -o beetle_v1_busco \ -m genome \ -c 24 \ --force结果出来C: 78.9%[S: 62.1%, D: 16.8%], F: 5.6%, M: 15.5%, n: 1013光看 C78.9%对于一个初版组装似乎还行但 D16.8% 太刺眼M15.5% 也不小。我心里基本有数这个组装大概率把不少杂合单倍型拼成了冗余序列D 高说明冗余严重而某些因冗余而“打架”的区段可能根本没组好所以 M 也偏高。这种 N50 和完整度背离的情况在杂合度高的物种里特别常见N50 漂亮不一定代表结构正确。5.2 purge去冗余后的第二次评估我给的建议是先做一遍 purge_dups。它的大致思路是用长读长把 reads 比对回组装基因组根据覆盖度差异识别哪些 contig 是冗余单倍型片段然后从主组装里剔除。这个工具在 HiFi 数据下表现很好对高杂合度物种几乎是标配步骤。朋友手头正好有 HiFi 数据处理很顺利。purge 之后重新组装了一版再跑一轮 BUSCOC: 91.2%[S: 88.5%, D: 2.7%], F: 3.1%, M: 5.7%, n: 1013这组数据的变化很值得分析C 从 78.9% 涨到 91.2%说明去冗余不只是让序列变短反而把之前因冗余冲突而丢失的基因恢复了。D 从 16.8% 跌到 2.7%说明绝大部分重复 BUSCO 命中确实是冗余单倍型造成的假象。M 下降到 5.7%这是可以接受的缺口剩余缺失大概率来自复杂重复区域。F 降到 3.1%碎片化问题也明显缓解。后来这版基因组做基因注释注释完整性明显提升结果跟 BUSCO 趋势对得上。这印证了一件事BUSCO 的 C 和 D 往往联动D 高不能只看成“统计洁癖”它往往掩盖着真实的组装结构问题。5.3 跟注释结果互相印证再分享一个用了蛋白模式的小技巧当基因组版本基本稳定后我会用注释得到的蛋白集再跑一次 proteins 模式作为交叉验证。busco -i gene_annotation.proteins.fasta \ -l arthropoda_odb10 \ -o beetle_v2_proteome_busco \ -m proteins \ -c 24正常情况下蛋白模式的 C 会略低于 genome 模式的 C因为注释流程本身会漏掉一部分基因。如果蛋白模式 C 掉得太多比如相对 genome 模式低了 10 个百分点说明注释软件或参数该调了。这个方法比翻一堆注释统计报告更直观因为 BUSCO 基因本身就是一套固定清单两头一对比注释环节的损失一目了然。想评估自己的注释流程有没有漏基因用这个方法准没错。6. 几个我踩过的坑版本、数据库、可复现性与写作规范6.1 版本不一致导致结果不可比这是最容易掉进去的坑而且一旦掉进去你手里的 C 值就是“薛定谔的 C”。BUSCO 4.x 和 5.x 在底层比对和数据库组织上有不少改动同一份基因组文件用 BUSCO 4.0.5 跑可能是 C: 92%换成 BUSCO 5.4.3 跑变成 89%这不代表组装变差了纯粹是尺子换了。我遇到过一次情况跨实验室合作时对方材料里写 BUSCO C: 91%我们本地复测只有 88%。排查半天才发现对方用的 BUSCO 3.0.2 加 ODB9 数据库我们用的是 5.4.3 加 ODB10。版本和数据库都是变量结果当然对不上。所以我的建议很明确项目记录里固定一个 BUSCO 版本比如全流程统一用 5.4.3。升级 BUSCO 之后不要拿旧数据做直接对比除非你把全部样本重跑一遍。提交结果或写文章时把 BUSCO version、数据库版本、谱系名全部写清楚。6.2 谱系选错导致误判另一个容易犯的错是“一个库打天下”。有些自动化 pipeline 为了省事不管什么输入都用 eukaryota_odb10这非常不严谨。eukaryota_odb10 覆盖所有真核生物里面的 BUSCO 都是极端保守的核心基因数量少缺少谱系特有基因。用它评估节肢动物基因组很可能显示 C 高达 95%看起来很漂亮但根本不能代表基因组完整——很多节肢动物谱系里保守的基因根本不在库内缺失了也查不出来。正确的做法是先确认物种分类位置再选最匹配的谱系。要对齐一个团队不同项目之间的结果也要统一谱系选择策略。如果你不确定用busco --list-datasets看看有哪些选择再根据物种所属的门纲目标去选。这一点值得反复强调选对谱系结果才有意义。6.3 文章里怎么写BUSCO参数才够规范最后说一个发表相关的事。现在基因组学文章里 BUSCO 结果几乎是标配但很多 Methods 部分写得含糊给审稿人留下扣分点。规范的写法应该包含BUSCO 软件版本例如 v5.4.3。数据库版本例如 ODB10。谱系名例如 arthropoda_odb10最好带上库内 BUSCO 数量 n。运行模式genome、proteins 还是 transcriptome。关键参数线程数、是否用了 MetaEuk/Augustus、是否做了 soft-mask。结果一行式例如C: 91.2%[S: 88.5%, D: 2.7%], F: 3.1%, M: 5.7%, n: 1013。我在 Methods 里一般会写这样一段Genome assembly completeness was assessed using BUSCO v5.4.3 with the arthropoda_odb10 (n1013) lineage dataset in genome mode. The final assembly showed C: 91.2%[S: 88.5%, D: 2.7%], F: 3.1%, M: 5.7% (n1013).这一小段信息量完整别人拿着你的原始数据也能完全复现。复现实验是审稿人最爱做的事之一可别在这些细节上丢分。最后再分享一条我在多个项目里反复验证的原则BUSCO 是很好的“完整性仪表盘”但它不是万能牌。它回答的是“保守基因有没有找齐”回答不了“重复序列组装得好不好”“结构变异正确不正确”也替代不了真实基因注释。所以正确用法是把它嵌进整个组装评估链路里和 N50、Hi-C 挂载率、注释统计、reads 回比率一起看。只要记住这一点BUSCO 就是你判断一个基因组能否进入下游分析的最快路径。