BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs • JuniorTree

BUSCO（Benchmarking Universal Single-Copy Orthologs, 基准通用单拷贝同源基因）通过基于进化上对基因内容的预期，来定量评估基因组组装和注释的完整性。

BUSCO 前基因组组装质量的检验#

在 BUSCO 普及之前，基因组组装质量的检验主要依赖于纯技术或统计学指标，这些指标主要反映测序或组装方法的特征

片段连续性指标（N50）：基因组中至少有一半的序列是组装在长度等于或大于 N50 的片段（contigs）上的
错误率与分布评估： 包括单碱基错误率（per-base error rates）、插入片段大小分布（insert size distributions），以及通过 k-mer 分布来评估基因组偏差等
早期的基因评估尝试： 在 BUSCO 之前，曾有一种名为 CEGMA（Core Eukaryotic Gene Mapping Approach）的工具，它通过比对 248 个保守的核心真核基因来尝试评估基因组完整性

为什么需要 BUSCO#

传统的以 N50 为代表的统计学指标，无法从“基因内容（gene content）”的维度来评估基因组组装的完整性，作者在论文中提到 N50 和 BUSCO 的相关性很差，这说明即使 N50 很高（序列很长），但是也有可能丢掉关键基因

因此我们需要 BUSCO 不仅仅是需要基因组能拼得长，更需要生物学上的基因完整度

BUSCO 的生物学原理#

BUSCO 的核心在于利用那些在进化上高度保守、且在绝大多数物种中只以单拷贝形式存在的同源基因

同源基因 (Orthologs)： 指的是在不同物种中，由同一个共同祖先基因经过物种形成（speciation）而演化来的基因
进化上高度保守： 意味着在漫长的进化过程中，这些基因在某个庞大的系统发育分支（如脊椎动物、节肢动物、真菌等）的各个物种里都被完整地保留了下来，它们通常负责非常基础且核心的生命活动
在绝大多数物种中只以单拷贝形式存在： 研究表明，这类基因在进化上受到“单拷贝控制（single-copy control）”，通过对数百个物种的基因组进行取样分析，研究人员发现这些基因在超过 90% 的被测物种中都只存在唯一的一个拷贝，真正的基因重复（duplication）在这些基因上极为罕见

单拷贝 ≠ 序列保守#

单拷贝（single-copy）描述的是基因组中的拷贝数，而保守性（conservation）描述的是序列在进化中的变化程度，这是两个完全独立的维度

一个基因可以是：

	单拷贝	多拷贝
高度保守	组蛋白H3（拷贝少但序列几乎不变）	rRNA基因（多拷贝且高度保守）
快速演化	很多单拷贝免疫基因	嗅觉受体基因家族

单拷贝是基因组结构的描述；保守性是进化动力学的描述，前者可以提示后者（单拷贝基因承受不起功能丧失，可能受强负选择），但绝不等同于后者

BUSCO 的算法实现#

针对基因组组装序列#

初步候选定位 (约占总运行时间的 15%)： 算法首先调用 tBLASTn 工具，利用预先构建好的 BUSCO 群体共有序列（consensus sequences），在待评估的基因组测序数据中进行全基因组搜索，从而定位出可能包含目标单拷贝基因的候选基因座（genomic loci）
**基因预测与提取 (约占总运行时间的 80%)：**定位到候选区域后，算法使用 Augustus 软件执行基因预测，在这一步中，算法会在特定的氨基酸 BUSCO 模块谱（block-profiles）的引导下，对这些候选基因座进行基因结构的注释和序列提取
**核心同源性评估与分类 (约占总运行时间的 5%)：**最后，算法将预测出的基因序列交由 HMMER 3 软件处理，HMMER 3 利用从氨基酸多重比对中建立的隐马尔可夫模型（HMM）谱进行打分，来评估这些匹配到的基因是否为真正的同源基因（即得分是否达到特定 BUSCO 群体的阈值/bitscore cut-offs），并将确认为同源基因的结果分类为“完整（Complete）”或“碎片化（Fragmented）”等指标

针对转录组数据#

如果评估的是转录组数据，算法会省略 tBLASTn 和 Augustus 的复杂步骤：

算法采用寻找最长开放阅读框 (longest open reading frame, ORF) 的直接方法来进行初始的基因预测
随后，直接将预测出的 ORF 序列送入 HMMER 3 流程进行核心同源性评估

针对已注释基因集#

如果输入的数据已经是提取完毕的注释基因集：

算法会完全跳过初步定位和基因预测步骤，直接把这些基因输入给 HMMER 3 执行核心的同源性评估和分类

源代码赏析时间#

这里的代码是v6的，托管在GitLab上：https://gitlab.com/ezlab/busco ↗

我比较好奇 HMMER 比对是怎么做的，所以就只看了这一段的代码

HMMER 这一步比对的不是原始基因组序列，而是先从基因组里预测出来的蛋白序列：

def run_analysis(self):
    """
    This function calls all needed steps for running the analysis.
    """
    super().run_analysis()
    self._run_prodigal() # 先跑 prodigal
    self.gene_details = self.prodigal_runner.gene_details # 预测蛋白对应的基因坐标信息
    self.run_hmmer(self.prodigal_runner.output_faa) # 将预测出来的蛋白文件再传给 hmmer
    self.hmmer_runner.write_buscos_to_file()
    return

python

这一步是在组装 HMMER 的命令

def configure_job(self, busco_id, seq_filename, output_filename):

        hmmer_job = self.create_job()
        hmmer_job.add_parameter("--domtblout")
        hmmer_job.add_parameter(os.path.join(self.results_dir, output_filename))
        hmmer_job.add_parameter("--cpu")
        hmmer_job.add_parameter("1")
        hmmer_job.add_parameter(
            os.path.join(self.lineage_dataset, "hmms", "{}.hmm".format(busco_id))
        )
        hmmer_job.add_parameter(seq_filename)
        return hmmer_job

python

最后组装出来的命令是：

hmmsearch \
    --domtblout <results_dir/output_filename> \
    --cpu 1 \
    <lineage_dataset/hmms/<busco_id>.hmm> \
    <seq_filename>

bash

BUSCO 给的参考的数据集中会有很多的 hmm 模型，比如我测试使用的是 lepidoptera_odb12 ，解压会给一个 dataset.cfg ，其中的 number_of_BUSCOs 就是 hmm 模型的数量，也就是5760个保守基因家族

name=lepidoptera_odb12
domain=eukaryota
creation_date=2025-07-01
number_of_BUSCOs=5760
number_of_species=79
ncbi_taxid=7088
max_intron=130000
max_seq_len=160000
species=fly
OrthoDB_version=12.0
dataset_version=01

plaintext

所以真正的流程应该是：

基因组 → 预测蛋白组 → 用每个 BUSCO 的 HMM 去搜这些蛋白 → 判断这个 BUSCO 是完整、重复、碎片化还是缺失

输入基因组：contig1，contig2，contig3

基因预测后得到蛋白：P1, P2, P3, P4, P5, P6

BUSCO 数据集里有 4 个 marker：B1, B2, B3, B4

然后分别搜索：

B1.hmm 去搜 P1..P6，发现 P2 是完整命中 -> B1 = Single-copy
B2.hmm 去搜 P1..P6，发现 P3 和 P5 都是完整命中 -> B2 = Duplicated
B3.hmm 去搜 P1..P6，只发现 P4 命中了半截 -> B3 = Fragmented
B4.hmm 去搜 P1..P6，一个都没找到 -> B4 = Missing

这里是解析表格的：

def parse_hmmer_output(self, filename, busco_query):

    processed_genes = []
    matched_genes = []
    records = defaultdict(list)

    data = []
    with open(filename, "r") as f:
        lines = f.readlines()

    for line in lines:
        if line.startswith("#"): # 跳过注释行
            continue
        else:
            values = line.split() # 把每一行拆成字段
            description = "".join(values[22:])
            row_data = tuple(values[:22]) + (description,)
            data.append(row_data)

    arr = np.array(
        data,
        dtype=[
            ("target_name", "U500"),
            ("target_accession", "U500"),
            ("tlen", "i4"),
            ("query_name", "U500"),
            ("query_accession", "U500"),
            ("qlen", "i4"),
            ("eval", "f8"),
            ("score", "f8"),
            ("bias", "f8"),
            ("dom_num", "i4"),
            ("ndom", "i4"),
            ("c-eval", "f8"),
            ("i-eval", "f8"),
            ("dom_score", "f8"),
            ("dom_bias", "f8"),
            ("hmm_coord_from", "i4"),
            ("hmm_coord_to", "i4"),
            ("ali_coord_from", "i4"),
            ("ali_coord_to", "i4"),
            ("env_coord_from", "i4"),
            ("env_coord_to", "i4"),
            ("acc", "f8"),
            ("description", "U500"),
        ],
    )

    target_hits_sorted, hmm_match_lengths = self.sort_hmmer_results(arr) # 排好序的基因ID列表以及每个基因的HMM匹配长度

    for target in target_hits_sorted:
        all_domain_hits = arr[arr["target_name"] == target]
        gene_id = self.get_gene_id(target)
        score = all_domain_hits[0]["score"]

        if not self.mode == "proteins" and not target in processed_genes:
            if self._check_overlap(matched_genes, gene_id):
                continue

        # Extract frame information (present in transcriptome mode)
        frame = (
            all_domain_hits[0]["description"]
            if "frame" in all_domain_hits[0]["description"]
            else None
        )
        # Store bitscore matches for each gene match. If match below cutoff, discard.
        if score < float(self.cutoff_dict[busco_query]["score"]):
            continue

        hmm_coords = [
            (hit["hmm_coord_from"], hit["hmm_coord_to"]) for hit in all_domain_hits
        ]
        env_coords = [
            (hit["env_coord_from"], hit["env_coord_to"]) for hit in all_domain_hits
        ]

        records[gene_id].append(
            {
                "hmm_matched_length": hmm_match_lengths[target],
                "hmm_profile_length": all_domain_hits[0]["qlen"],
                "hmm_coords": hmm_coords,
                "env_coords": env_coords,
                "score": score,
                "frame": frame,
                "orig gene ID": target,
                "ref gene ID": target,
            }
        )

        if gene_id not in matched_genes:
            matched_genes.append(gene_id)
            processed_genes.append(target)
    return records

python

计算蛋白的 HMM 覆盖总长度：

def sum_hmm_len(self, hmm_coords):
    hmm_regions_sorted = sorted(hmm_coords, key=lambda x: x[0]) # 按照坐标对的第一个元素（起始位置）进行升序排序
    hmm_regions_used = []
    for region in hmm_regions_sorted:
        for used in hmm_regions_used:
            if (
                used[0] <= region[0] <= used[1]
            ):  # if any new coord contained in previous region, expand if necessary
                if region[1] > used[1]:
                    used[1] = region[1]
                break
        else:
            hmm_regions_used.append(list(region))

    return sum([region[1] - region[0] + 1 for region in hmm_regions_used]) # 按照闭区间返回结果

python

然后是对 match 到的进行排序：

def _sort_matches(self, matched_record, busco_query):
	# 元组的初始化...
	if self.config.get("busco_run", "datasets_version") == "odb10":
		# odb10数据集（使用统计学z-score方法）
	else:
	 # odb12及以后数据集（使用百分比方法）
	 if size >= 0.8 * profile_length:
      complete
   else:
      fragment
	return (
            busco_complete,
            busco_vlarge,
            busco_fragment,
            matched_genes_complete,
            matched_genes_vlarge,
            matched_genes_fragment,
        )

python

这里会把命中分为三堆：

is_complete
is_very_large
is_fragment

假设得到的是：

is_complete = {
  "B1": {
    "P1": [{"bitscore": 200, "length": 461}]
  },
  "B2": {
    "P3": [{"bitscore": 180, "length": 430}],
    "P4": [{"bitscore": 170, "length": 421}]
  }
}
is_very_large = {
  "B5": {
    "P9": [{"bitscore": 210, "length": 530}]
  }
}
is_fragment = {
  "B3": {
    "P5": [{"bitscore": 130, "length": 181}],
    "P6": [{"bitscore": 108, "length": 210}]
  }
}

python

意思是：

B1 有 1 个完整命中 P1
B2 有 2 个完整命中 P3/P4
B5 有 1 个 very_large 命中 P9
B3 有 2 个 fragment 命中 P5/P6

第 2 步：filter() 之后

这里会去重，并去掉低分命中。

假设 B3 里 P6 被删掉了：

is_complete = {
  "B1": {
    "P1": [{"bitscore": 200, "length": 461}]
  },
  "B2": {
    "P3": [{"bitscore": 180, "length": 430}],
    "P4": [{"bitscore": 170, "length": 421}]
  }
}
is_very_large = {
  "B5": {
    "P9": [{"bitscore": 210, "length": 530}]
  }
}
is_fragment = {
  "B3": {
    "P5": [{"bitscore": 130, "length": 181}]
  }
}

python

第 3 步：consolidate_busco_lists()

它做的事其实很简单：

先看 is_complete 和 is_very_large

\text{一个 BUSCO 只剩 1 个命中} \to \text{single\_copy\_buscos}

\text{一个 BUSCO 剩多个命中} \to \text{multi\_copy\_buscos}

再看 is_fragment

如果一个 BUSCO 有多个 fragment，只保留得分最高的那个
然后放进 fragmented_buscos

所以，上面的数据会变成：

single_copy_buscos = {
  "B1": {
    "P1": [{"bitscore": 200, "length": 461}]
  },
  "B5": {
    "P9": [{"bitscore": 210, "length": 530}]
  }
}
multi_copy_buscos = {
  "B2": {
    "P3": [{"bitscore": 180, "length": 430}],
    "P4": [{"bitscore": 170, "length": 421}]
  }
}
fragmented_buscos = {
  "B3": {
    "P5": [{"bitscore": 130, "length": 181}]
  }
}

python

最终指标的计算，把每个 BUSCO 家族最终分到的类别做计数和百分比：

# 结果做的结构化
def record_results(self):
    self._get_busco_percentages()
    extra = ",E:{}%".format(self.e_percent) if self.e_percent else ""
    self.one_line_summary_raw = (
        "C:{}%[S:{}%,D:{}%],F:{}%,M:{}%,n:{}{}\t{}\n".format(
            self.complete_percent,
            self.s_percent,
            self.d_percent,
            self.f_percent,
            self.missing_percent,
            self.total_buscos,
            extra,
            "   ",
        )
    )
    self.one_line_summary = "Results:\t{}".format(self.one_line_summary_raw)

python

def _get_busco_percentages(self):
    self.single_copy = len(self.single_copy_buscos)  # 单拷贝 BUSCO 家族数
    self.multi_copy = len(self.multi_copy_buscos)  # 重复 BUSCO 家族数
    self.only_fragments = len(self.fragmented_buscos)  # 碎片化 BUSCO 家族数
    self.total_buscos = len(self.cutoff_dict) # 总 BUSCO 家族数
    self.num_missing = (
        self.total_buscos - self.single_copy - self.multi_copy - self.only_fragments
    ) # 缺失 BUSCO 家族数

    # Get percentage of each kind of BUSCO match
    self.s_percent = abs(round((self.single_copy / self.total_buscos) * 100, 1))
    self.d_percent = abs(round((self.multi_copy / self.total_buscos) * 100, 1))
    self.f_percent = abs(round((self.only_fragments / self.total_buscos) * 100, 1))
    self.complete_percent = abs(
        round(((self.single_copy + self.multi_copy) / self.total_buscos) * 100, 1)
    )
    self.missing_percent = abs(
        round((self.num_missing / self.total_buscos) * 100, 1)
    )

    if self.miniprot_pipeline:
        self._get_stop_codon_percent_and_avg_identity()
    else:
        self.e_percent = 0
        self.avg_identity = None
        self.complete_stop_codon_count = 0

    return

python

所以 BUSCO 本质上回答的是：“这份基因组里，预期应该存在的保守基因家族，找回了多少？”

总结#

看了论文和源码，感觉豁然开朗，单从底层计算工具的角度来看，BUSCO 确实没有发明全新的序列比对或基因预测算法，它本质上是一个整合了现有优秀软件的自动化分析流水线（Pipeline）

BUSCO 的伟大之处和核心专有贡献并不在于“写了一个新搜索代码”，而在于它提供了“标尺”和“规则”，可以说，没有 BUSCO 团队前期的海量生物学工作，这些工具根本不知道该去寻找什么

研究团队并没有让工具去瞎找，而是基于 OrthoDB 数据库，对数百个物种的基因组进行了极其庞大的比较基因组学分析，他们大海捞针般地筛选出了在超过 90% 的物种中都保持单拷贝的同源基因，并经过严格的唯一性和保守性过滤，最终为不同的进化分支“量身定制”了专属的测试集（例如脊椎动物的 3023 个基因、节肢动物的 2675 个基因、真菌的 1438 个基因等）

然 BUSCO 调用了 HMMER 3 软件，但 HMMER 3 只是一个计算器。BUSCO 团队为这成千上万个单拷贝基因构建了基于氨基酸多重比对的隐马尔可夫模型谱，并且为每一个 BUSCO 基因群设定了专属的特异性得分阈值（bitscore cut-offs），只有当 HMMER 3 算出的得分满足了 BUSCO 设定的这个严苛阈值，才会被判定为真正的同源基因

对于结果 C[D], F, M, n ，在统计学上给出了明确的定义

另外比起 N50 这种统计学算法或者理论，BUSCO 结合了生物学意义，给出了可解释性的理论