宏基因组功能注释——eggNOG-mapper实战¶

eggNOG v5.0（eggNOG-mapper v2 默认数据库）
├── root (LUCA level)
├── Bacteria (tax_id=2)
│   ├── Proteobacteria (tax_id=1224)
│   │   ├── Gammaproteobacteria
│   │   └── Alphaproteobacteria
│   └── Firmicutes (tax_id=1239)
├── Archaea (tax_id=2157)
└── Eukaryota (tax_id=2759)

# 方法1：conda安装（推荐）
conda create -n emapper python=3.9
conda activate emapper
conda install -c bioconda eggnog-mapper

# 方法2：pip安装
pip install eggnog-mapper

# 方法3：源码安装
git clone https://github.com/eggnogdb/eggnog-mapper.git
cd eggnog-mapper
pip install .

# 下载数据库（必须！约45GB）
download_eggnog_data.py --data_dir /data/eggnog_db/ -y

# 下载DIAMOND数据库（核心）
download_eggnog_data.py -D --data_dir /data/eggnog_db/ -y

# 下载特定tax scope的HMM数据库（可选，用于hmmer模式）
download_eggnog_data.py -H -d 2 --data_dir /data/eggnog_db/ -y  # bacteria
download_eggnog_data.py -H -d 2157 --data_dir /data/eggnog_db/ -y  # archaea

# 验证安装
emapper.py --version
emapper.py --data_dir /data/eggnog_db/ --list_taxa

# 从组装的contigs预测ORF（使用Prodigal）
prodigal -i assembled_contigs.fa \
         -o genes.gff \
         -a proteins.faa \        # 蛋白质序列输出
         -d genes.fna \           # 核苷酸序列输出
         -p meta \                # 宏基因组模式
         -f gff

# 检查蛋白序列数量
grep -c ">" proteins.faa

# 清理序列头（去除Prodigal额外信息）
sed 's/ .*//' proteins.faa > proteins_clean.faa

# 大文件拆分
# 拆成每份10万条序列
awk 'BEGIN{n=0; file="split_00.faa"} /^>/{n++; if(n%100000==1){file=sprintf("split_%02d.faa",n/100000)}} {print > file}' proteins_clean.faa

# 或使用seqkit拆分
seqkit split2 -s 100000 proteins_clean.faa -O split_dir/

# === 基本运行 ===
emapper.py -i proteins_clean.faa \
           --output sample1 \
           --output_dir /results/emapper/ \
           --data_dir /data/eggnog_db/ \
           --cpu 16 \
           --override

# === 宏基因组推荐参数 ===
emapper.py -i proteins_clean.faa \
           -m diamond \
           --itype proteins \
           --tax_scope Bacteria \
           --go_evidence non-electronic \
           --target_orthologs all \
           --seed_ortholog_evalue 1e-5 \
           --seed_ortholog_score 60 \
           --query_cover 20 \
           --subject_cover 20 \
           --output sample1_emapper \
           --output_dir /results/emapper/ \
           --data_dir /data/eggnog_db/ \
           --cpu 32 \
           --override \
           --temp_dir /tmp/emapper/

# === 分步运行（大数据集推荐） ===
# 步骤1：搜索阶段（最耗时）
emapper.py -i proteins_clean.faa \
           -m diamond \
           --itype proteins \
           --data_dir /data/eggnog_db/ \
           --cpu 32 \
           --no_annot \
           --output sample1 \
           --output_dir /results/emapper/ \
           --override

# 步骤2：注释阶段（基于步骤1的结果）
emapper.py --annotate_hits_table /results/emapper/sample1.emapper.seed_orthologs \
           --data_dir /data/eggnog_db/ \
           --cpu 16 \
           --go_evidence non-electronic \
           --tax_scope Bacteria \
           --output sample1 \
           --output_dir /results/emapper/ \
           --override

# === 使用MMseqs2（更快） ===
emapper.py -i proteins_clean.faa \
           -m mmseqs \
           --data_dir /data/eggnog_db/ \
           --cpu 32 \
           --output sample1_mmseqs \
           --output_dir /results/emapper/ \
           --override

# 查看输出文件结构
head -5 sample1_emapper.emapper.annotations

# 统计注释率
total=$(grep -c ">" proteins_clean.faa)
annotated=$(grep -v "^#" sample1_emapper.emapper.annotations | wc -l)
echo "Annotation rate: $annotated / $total = $(echo "scale=2; $annotated*100/$total" | bc)%"

# 提取COG分类统计
grep -v "^#" sample1_emapper.emapper.annotations | \
  cut -f7 | grep -v "^-$" | grep -o . | sort | uniq -c | sort -rn

# 提取KEGG KO列表
grep -v "^#" sample1_emapper.emapper.annotations | \
  cut -f12 | grep -v "^-$" | tr ',' '\n' | sort -u > kegg_ko_list.txt

# 提取GO terms
grep -v "^#" sample1_emapper.emapper.annotations | \
  cut -f10 | grep -v "^-$" | tr ',' '\n' | sort -u > go_terms_list.txt

# R语言处理emapper结果
library(tidyverse)

# 读取注释结果
emapper <- read_tsv("sample1_emapper.emapper.annotations",
                    comment = "##",
                    col_names = TRUE)

# === COG分类统计与可视化 ===
cog_categories <- c(
  "J" = "Translation", "A" = "RNA processing",
  "K" = "Transcription", "L" = "Replication",
  "B" = "Chromatin", "D" = "Cell cycle",
  "Y" = "Nuclear structure", "V" = "Defense mechanisms",
  "T" = "Signal transduction", "M" = "Cell wall",
  "N" = "Cell motility", "Z" = "Cytoskeleton",
  "W" = "Extracellular structures", "U" = "Intracellular trafficking",
  "O" = "PTM/chaperones", "X" = "Mobilome",
  "C" = "Energy production", "G" = "Carbohydrate metabolism",
  "E" = "Amino acid metabolism", "F" = "Nucleotide metabolism",
  "H" = "Coenzyme metabolism", "I" = "Lipid metabolism",
  "P" = "Inorganic ion transport", "Q" = "Secondary metabolites",
  "R" = "General function", "S" = "Function unknown"
)

# 提取并统计COG
cog_data <- emapper %>%
  filter(COG_category != "-") %>%
  mutate(COG_split = strsplit(COG_category, "")) %>%
  unnest(COG_split) %>%
  count(COG_split) %>%
  mutate(description = cog_categories[COG_split])

# 画图
ggplot(cog_data, aes(x = reorder(COG_split, -n), y = n, fill = COG_split)) +
  geom_bar(stat = "identity") +
  labs(x = "COG Category", y = "Gene Count",
       title = "COG Functional Classification") +
  theme_minimal() +
  theme(legend.position = "none")

# === KEGG通路分析 ===
# 提取KO编号
ko_list <- emapper %>%
  filter(KEGG_ko != "-") %>%
  separate_rows(KEGG_ko, sep = ",") %>%
  mutate(KEGG_ko = str_remove(KEGG_ko, "ko:"))

# 统计各KO出现频次
ko_counts <- ko_list %>%
  count(KEGG_ko, sort = TRUE)

# 使用clusterProfiler做KEGG富集（需要转换）
library(clusterProfiler)
# 将KO映射到pathway
# enrichKEGG需要KEGG gene ID或KO ID

# === 多样本比较 ===
# 假设有多个样本的emapper结果
samples <- c("sample1", "sample2", "sample3")
all_cog <- map_dfr(samples, function(s) {
  read_tsv(paste0(s, "_emapper.emapper.annotations"), comment = "##") %>%
    filter(COG_category != "-") %>%
    mutate(sample = s)
})

# Python处理emapper结果
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 读取注释结果
df = pd.read_csv("sample1_emapper.emapper.annotations",
                 sep="\t", comment="#",
                 header=0)

# KEGG模块统计
kegg_modules = df['KEGG_Module'].dropna()
kegg_modules = kegg_modules[kegg_modules != '-']
module_counts = kegg_modules.str.split(',').explode().value_counts().head(20)

# 可视化
fig, ax = plt.subplots(figsize=(10, 6))
module_counts.plot(kind='barh', ax=ax)
ax.set_xlabel('Gene Count')
ax.set_title('Top 20 KEGG Modules')
plt.tight_layout()
plt.savefig('kegg_modules.pdf')

# ===== 一站式宏基因组功能注释流程 =====

# 1. 从组装结果预测蛋白
prodigal -i final_contigs.fa -a proteins.faa -p meta -f gff -o genes.gff

# 2. 清理序列头
sed 's/ .*//' proteins.faa > proteins_clean.faa

# 3. 运行eggNOG-mapper（DIAMOND模式）
emapper.py -i proteins_clean.faa \
  -m diamond \
  --itype proteins \
  --data_dir /data/eggnog_db/ \
  --tax_scope Bacteria \
  --go_evidence non-electronic \
  --seed_ortholog_evalue 1e-5 \
  --seed_ortholog_score 60 \
  --cpu 32 \
  --output metagenome_func \
  --output_dir ./emapper_results/ \
  --temp_dir /tmp/emapper_tmp/ \
  --override

# 4. 统计注释率
echo "Total proteins: $(grep -c '>' proteins_clean.faa)"
echo "Annotated: $(grep -v '^#' emapper_results/metagenome_func.emapper.annotations | wc -l)"

# 5. 提取各功能分类
# COG
grep -v '^#' emapper_results/metagenome_func.emapper.annotations | \
  awk -F'\t' '$7!="-"' | cut -f1,7 > cog_assignments.tsv

# KEGG KO
grep -v '^#' emapper_results/metagenome_func.emapper.annotations | \
  awk -F'\t' '$12!="-"' | cut -f1,12 > kegg_ko_assignments.tsv

# GO
grep -v '^#' emapper_results/metagenome_func.emapper.annotations | \
  awk -F'\t' '$10!="-"' | cut -f1,10 > go_assignments.tsv

# PFAM
grep -v '^#' emapper_results/metagenome_func.emapper.annotations | \
  awk -F'\t' '$21!="-"' | cut -f1,21 > pfam_assignments.tsv

# 6. 可视化（使用R或Python脚本处理上述TSV文件）