13 Mar 2025 4479字 15分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

前面使用了eggNOG-mapper快速注释功能基因，我们最后得到了很多结果文件，其中最重要的是两个annotations文件。这里主要讲一下怎么整理结果文件，并且对注释的结果做质量评估。

1. 基因功能注释评估

1.1. 表格信息整理

其实就是整理结果文件中有多少基因注释到了哪些数据库中，以一种直观的方式展现结果。

这里分两种情况，如果自己喜欢以编程的方式处理文件，可以直接处理out.emapper.annotations这个文件。如果前一种方式不是很得心应手的话，那就直接处理out.emapper.annotations.xlsx这个文件，毕竟可以用excel直接打开。

两种处理方式

1. 编程处理：直接解析 out.emapper.annotations 文件，提取有用信息。

2. Excel 处理：使用 out.emapper.annotations.xlsx 文件，利用 Excel 进行筛选和整理。

步骤（以 Excel 处理为例）

1. 删除前两行（运行参数和命令），保留表头及数据部分。

2. 筛选和提取相关信息：

• 基因名（query）：A 列，所有基因名称。

• GO 注释：J 列，去除 - 后，复制 A 列内容（代表这些基因有 GO 注释）。

• KEGG 注释：L 列，去除 - 后，复制 A 列内容（代表这些基因有 KEGG 注释）。

• PFAM 注释：U 列，去除 - 后，复制 A 列内容（代表这些基因有 PFAM 注释）。

可以得到如下形式的表格（另存为annotation.xlsx):

1.2. venny图绘制

可以使用R，用经典的VennDiagram包做韦恩图，或者用UpSetR包做Upset图：

#!/usr/bin/env Rscript
library(VennDiagram)
library(readxl)

# 读取 Excel 文件
data <- read_excel("hap2_annotation.xlsx", sheet = 1)

# 替换 NA 为 ""，然后转换为列表并去重、去空值
data[is.na(data)] <- ""

# 提取非空值的基因 ID 并去重
eggNOG_genes <- unique(na.omit(data$eggNOG[data$eggNOG != ""]))
GO_genes <- unique(na.omit(data$GO[data$GO != ""]))
KEGG_genes <- unique(na.omit(data$KEGG[data$KEGG != ""]))
PFAM_genes <- unique(na.omit(data$PFAM[data$PFAM != ""]))

# 保存 Venn 图到文件
venn.diagram(
  x = list(
    eggNOG = eggNOG_genes,
    GO = GO_genes,
    KEGG = KEGG_genes,
    PFAM = PFAM_genes
  ),
  filename = "hap1_Venny_annotation.png",  # 直接保存文件
  output = TRUE,
  col = "black",
  lwd = 3,
  lty = "solid",
  fill = c("cornflowerblue", "green", "yellow", "darkorchid1"),
  alpha = 0.5,
  label.col = c("orange", "white", "darkorchid4", "white",
                "white", "white", "white", "white", "darkblue", "white",
                "white", "white", "white", "darkgreen", "white"),
  cex = 2.0,
  fontfamily = "serif",
  fontface = "bold",
  cat.col = c("darkblue", "darkgreen", "orange", "darkorchid4"),
  cat.cex = 1.5,
  cat.dist = c(0.05, 0.05, 0.05, 0.05),  # 调整标签位置，防止重叠
  rotation.degree = 0,
  cat.fontfamily = "serif",
  cat.fontface = "bold",
  margin = 0.1
)

# 提示用户图像已保存
message("Venn 图已保存为 hap2_Venny_annotation.png")

1.3. upset图绘制

#!/usr/bin/env Rscript
library(UpSetR)
library(readxl)
data = read_excel("hap1_annotation.xlsx", sheet = 1)
data[is.na(data)] = ""  # 替换NA为空
data = fromList(data)  # 转换成0和1的矩阵（UpSetR包对导入数据有要求，建议看github官网）

upset(
  data,
  nsets = 4,
  matrix.color = "gray23",
  main.bar.color = "gray23",
  mainbar.y.label = "交集基因数",
  mainbar.y.max = 7000,
  sets.bar.color = "gray23",
  sets.x.label = "注释的基因数",
  point.size = 2.2, 
  line.size = 0.7,
  mb.ratio = c(0.7, 0.3),
  show.numbers = "yes",
  set_size.show = FALSE
)

UpSetR参数参考UpSetR: A More Scalable Alternative to Venn and Euler Diagrams for Visualizing Intersecting Sets (r-project.org)

不过，两个图展示的信息其实是一样的（如果不理解upset图，可以对照韦恩图，通过柱状图的数字以及底部的点和线就能明白）。作图的目的不是为了炫耀技巧，而是为了让信息更加直观。upset图相对更清晰一些。当然，也可以直接整理成表格呈现。

2. BUSCO评估

2.1. busco下载

这边使用conda下载：

conda create -n busco_env
conda activate busco_env
conda install -c conda-forge -c bioconda busco

安装之后通过busco -h查看是否安装成功，如果提示缺什么软件就用conda补上。

通过busco --list-datasets可以查看当前有哪些物种的数据库，我的植物是香椿，这边选择陆生植物数据库（embryophyta_odb12），适用于所有高等植物，包括被子植物、裸子植物。

通过选定好植物数据库后，可以使用busco --download embryophyta_odb12下载该植物数据库。

busco的详细参数可以看官网的user guide User guide BUSCO v5.4.4 (ezlab.org)

简单讲一讲格式和能用到的参数：

busco -i [SEQUENCE_FILE] -l [LINEAGE] -o [OUTPUT_NAME] -m [MODE] [OTHER OPTIONS]
'''
主要参数：
 -i		序列文件位置
 -l		下载的同源物种保守基因数据库位置
 -o		输出文件名
 -m		模式，分为genome,proteins,transcriptome三种
 其他参数：
 --cpu		设置cpu数量
 --download		在线下载数据库，根据分类有"all"、"prokaryota"、"eukaryota"和"virus" （不推荐，速度慢）
 --offline		离线模式，不会更新数据库
'''

2.2. 运行busco评估注释完整度

接下来就可以使用一下参考脚本运行busco评估了:

#!/bin/bash

# 设置变量
DATABASE_PATH="/data/wangxingbin/dataBase/Busco/busco_downloads/lineages/embryophyta_odb12"
SEQUENCE_PATH="/data/wangxingbin/00.Project/01.Toona_Annotation/hap1/03.functionalGene/03.TEsorter/hap1_rmTE.aa"
OUTPUT_DIR="hap1"
CPU=8

# 确保 BUSCO 处于 conda 环境
CONDA_ENV="busco_env"

# 日志文件
LOG_FILE="${OUTPUT_DIR}_busco.log"

# 颜色输出
GREEN="\e[32m"
RED="\e[31m"
NC="\e[0m"  # 还原默认颜色

echo -e "${GREEN}[INFO] 激活 Conda 环境：${CONDA_ENV}${NC}"
source /data/wangxingbin/anaconda3/bin/activate ${CONDA_ENV}

# 检查数据库路径是否存在
if [[ ! -d "$DATABASE_PATH" ]]; then
    echo -e "${RED}[ERROR] BUSCO 数据库路径不存在：$DATABASE_PATH${NC}"
    exit 1
fi

# 检查输入序列文件是否存在
if [[ ! -f "$SEQUENCE_PATH" ]]; then
    echo -e "${RED}[ERROR] 输入序列文件不存在：$SEQUENCE_PATH${NC}"
    exit 1
fi

# 运行 BUSCO
echo -e "${GREEN}[INFO] 开始运行 BUSCO${NC}"
busco -i ${SEQUENCE_PATH} -l ${DATABASE_PATH} -o ${OUTPUT_DIR} -m proteins --cpu ${CPU} --offline | tee ${LOG_FILE}

# 检查 BUSCO 是否成功运行
if [[ $? -eq 0 ]]; then
    echo -e "${GREEN}[INFO] BUSCO 运行完成，结果存放在：${OUTPUT_DIR}${NC}"
else
    echo -e "${RED}[ERROR] BUSCO 运行失败，请检查日志文件：${LOG_FILE}${NC}"
    exit 1
fi

运行大概十来分钟，得到结果文件```short_summary.specific.embryophyta_odb12.hap1.txt``：

embryophyta_odb12高等植物数据库有2026个BUSCO groups，2014个（99.4%）BUSCO groups能够完整比对上（包括905个单拷贝和1109个多拷贝），说明这个基因组注释的完整性还是不错的。

结果解读：

• 高完整性：99.4% 的总 BUSCOs 完整性非常高，代表基因组注释相对完整。

• 单拷贝与重复拷贝：44.7% 是单拷贝的，这意味着在大部分情况下，每个基因组只包含一个副本，但有 54.7% 的基因组是重复的，可能与基因组的扩增或基因复制相关。

• 低片段化和缺失率：只有 1 个片段化的 BUSCO 和 11 个缺失的 BUSCO，基因组注释质量还是较好，仅有较少重要的基因丢失。

基因组注释从0开始（7）：功能基因注释质量评估