> For AI agents: the site index is available at https://www.juniortree.com/llms.txt; this page is available as Markdown at https://www.juniortree.com/paper-reading/busco-2015.md.

## Page metadata

- **Read:** 2026-03-29
- **Source URL:** https://www.juniortree.com/paper-reading/busco-2015

# busco-2015

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

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

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

1. 片段连续性指标（N50）：基因组中至少有一半的序列是组装在长度等于或大于 N50 的片段（contigs）上的
2. **错误率与分布评估：** 包括单碱基错误率（per-base error rates）、插入片段大小分布（insert size distributions），以及通过 k-mer 分布来评估基因组偏差等
3. **早期的基因评估尝试：** 在 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）”等指标

![](https://s3.juniortree.com/pic/2026/03/8788acf94880b4826455aa450f79f3d7.webp)

### 针对转录组数据

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

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

### **针对已注释基因集**

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

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

## 源代码赏析时间

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

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

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

```python
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
```

这一步是在组装 HMMER 的命令

```python
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
```

最后组装出来的命令是：

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

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
```

所以真正的流程应该是：

基因组 → 预测蛋白组 → 用每个 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

这里是解析表格的：

```python
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
```

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

```python
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]) # 按照闭区间返回结果
```

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

```python
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,
        )
```

这里会把命中分为三堆：

- `is_complete`
- `is_very_large`
- `is_fragment`

假设得到的是：

```python
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}]
  }
}
```

意思是：

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

第 2 步：`filter()` 之后

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

假设 B3 里 P6 被删掉了：

```python
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}]
  }
}
```

第 3 步：`consolidate_busco_lists()`

它做的事其实很简单：

1. 先看 `is_complete` 和 `is_very_large`

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

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

1. 再看 `is_fragment`

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

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

```python
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}]
  }
}
```

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

```python
# 结果做的结构化
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
```

所以 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 结合了生物学意义，给出了可解释性的理论