功能注释数据库（Functional Annotation Databases）¶

Level 1 (大类，共7类)
├── Metabolism（代谢）          ← 宏基因组最关注
├── Genetic Information Processing（遗传信息处理）
├── Environmental Information Processing（环境信息处理）
├── Cellular Processes（细胞过程）
├── Organismal Systems（有机体系统）
├── Human Diseases（人类疾病）
└── Drug Development（药物开发）

Level 2 (亚类，约50类)
├── Carbohydrate metabolism（碳水化合物代谢）
├── Amino acid metabolism（氨基酸代谢）
├── Lipid metabolism（脂质代谢）
└── ...

Level 3 (具体通路，约500条)
├── ko00010 Glycolysis/Gluconeogenesis（糖酵解/糖异生）
├── ko00020 Citrate cycle (TCA cycle)（三羧酸循环）
├── ko00030 Pentose phosphate pathway（磷酸戊糖途径）
└── ...

基因序列 → KO注释 → Pathway映射 → 功能富集分析
                                         ↓
                               比较不同样本/分组间的代谢差异

# ============================================================
# eggNOG-mapper: 综合功能注释（可同时得到COG/GO/KEGG/Pfam注释）
# ============================================================

# ---- Step 1: 下载eggNOG数据库 ----
# download_eggnog_data.py 会下载eggNOG v5.0数据库（约50GB，eggNOG-mapper v2 默认版本）
download_eggnog_data.py \
    --data_dir /path/to/eggnog_db/    # 数据库存放目录

# ---- Step 2: 运行注释（蛋白序列输入） ----
emapper.py \
    -i proteins.faa \                  # 输入：蛋白质序列（FASTA格式）
    -o output_prefix \                 # 输出文件前缀
    --data_dir /path/to/eggnog_db/ \   # eggNOG数据库路径
    -m diamond \                       # 比对方法：diamond（快速）或hmmer（敏感）
    --cpu 16 \                         # 使用16个CPU线程
    --override                         # 覆盖已有输出文件

# ---- Step 2b: 也可以直接注释核酸序列 ----
emapper.py \
    -i genes.fna \                     # 输入：核酸序列
    --itype CDS \                      # 输入类型：CDS（编码序列）
    -o output_prefix \
    --data_dir /path/to/eggnog_db/ \
    -m diamond \
    --cpu 16

# ---- 输出文件说明 ----
# output_prefix.emapper.annotations  ← 主结果文件，包含：
#   - COG类别（COG_category列）
#   - GO注释（GOs列）
#   - KEGG KO号（KEGG_ko列）
#   - KEGG Pathway（KEGG_Pathway列）
#   - Pfam域（PFAMs列）
#   - eggNOG OG（eggNOG_OGs列）
#   - 描述（Description列）

# 输出文件是TSV格式，关键列：
#  列1:  query             基因/蛋白ID
#  列5:  eggNOG_OGs        直系同源组
#  列7:  Description       功能描述
#  列10: GOs               GO注释
#  列12: KEGG_ko           KO号（如 ko:K00001）
#  列13: KEGG_Pathway      通路号（如 ko00010）
#  列14: KEGG_Module       模块号
#  列16: CAZy              CAZy注释
#  列21: COG_category      COG功能类别字母
#  列22: PFAMs             Pfam域

# ============================================================
# HUMAnN3: 宏基因组功能组成定量（MetaPhlAn配套工具）
# 输出通路丰度和基因家族丰度
# ============================================================

# ---- Step 1: 下载数据库 ----
# 下载ChocoPhlAn核酸数据库（微生物pangenome）
humann_databases --download chocophlan full /path/to/db/
# 下载UniRef蛋白数据库
humann_databases --download uniref uniref90_diamond /path/to/db/

# ---- Step 2: 运行HUMAnN3 ----
humann \
    --input sample_R1R2.fastq.gz \       # 输入：质控后的宏基因组reads
    --output humann3_output/ \           # 输出目录
    --nucleotide-database /path/to/db/chocophlan/ \  # ChocoPhlAn数据库
    --protein-database /path/to/db/uniref/ \         # UniRef数据库
    --threads 16 \                       # 线程数
    --metaphlan-options "--db_dir /path/to/metaphlan_db"  # MetaPhlAn 4.2+ 参数变更：--bowtie2db 已更名为 --db_dir

# ---- 输出文件（三个核心结果） ----
# 1. *_genefamilies.tsv     基因家族丰度（UniRef90/KO/EC编号）
# 2. *_pathabundance.tsv    通路丰度（MetaCyc通路）
# 3. *_pathcoverage.tsv     通路覆盖度

# ---- Step 3: 转换为KEGG注释 ----
# 将基因家族丰度从UniRef90转换为KO
humann_regroup_table \
    --input sample_genefamilies.tsv \    # 输入：基因家族丰度
    --output sample_ko.tsv \             # 输出：KO丰度
    --groups uniref90_ko                 # 分组方式：按KO号

# ---- Step 4: 标准化 ----
# CPM标准化（Copies Per Million，每百万拷贝数）
humann_renorm_table \
    --input sample_ko.tsv \              # 输入
    --output sample_ko_cpm.tsv \         # 输出
    --units cpm                          # 标准化单位

# ---- Step 5: 合并多样本 ----
humann_join_tables \
    --input humann3_output/ \            # 输入目录（包含所有样本结果）
    --output merged_ko.tsv               # 合并后的KO丰度表

# ============================================================
# dbCAN3 / run_dbcan v5: CAZy碳水化合物活性酶注释
# 整合三种方法：pyHMMER + DIAMOND + dbCAN_sub（v5用pyHMMER替代HMMER）
# ============================================================

# ---- Step 1: 下载dbCAN数据库 ----
# 下载HMM数据库、DIAMOND数据库等
dbcan_build --cpus 8 --db-dir /path/to/dbcan_db/

# ---- Step 2: 运行CAZy注释 ----
run_dbcan \
    proteins.faa \                       # 输入：蛋白质序列（FASTA）
    protein \                            # 输入类型：protein（蛋白序列）
    --db_dir /path/to/dbcan_db/ \        # dbCAN数据库路径
    --out_dir dbcan_output/ \            # 输出目录
    --dia_cpu 8 \                        # DIAMOND使用的CPU数
    --hmm_cpu 8 \                        # HMMER使用的CPU数
    --tf_cpu 8                           # eCAMI使用的CPU数

# ---- 输出文件 ----
# overview.txt         ← 汇总结果（最重要）
#   包含三种方法的结果，推荐至少2种方法同时注释到才可靠
# hmmer.out            ← HMMER结果
# diamond.out          ← DIAMOND结果
# eCAMI.out            ← 深度学习方法结果

# ============================================================
# ABRicate: 快速注释耐药基因、毒力因子等
# 支持多种数据库：CARD, ResFinder, VFDB, NCBI等
# ============================================================

# ---- 查看可用数据库 ----
abricate --list                          # 列出所有可用数据库

# ---- 耐药基因注释（使用CARD数据库） ----
abricate \
    --db card \                          # 数据库：card（耐药基因）
    --minid 80 \                         # 最低序列一致性：80%
    --mincov 80 \                        # 最低覆盖度：80%
    --threads 8 \                        # 线程数
    contigs.fasta \                      # 输入：组装序列
    > abricate_card_results.tsv          # 输出：TSV格式结果

# ---- 毒力因子注释（使用VFDB数据库） ----
abricate \
    --db vfdb \                          # 数据库：vfdb（毒力因子）
    --minid 70 \                         # 最低一致性：70%
    --mincov 70 \                        # 最低覆盖度：70%
    contigs.fasta \
    > abricate_vfdb_results.tsv

# ---- 多样本结果汇总 ----
abricate --summary \
    sample1_card.tsv \                   # 样本1结果
    sample2_card.tsv \                   # 样本2结果
    sample3_card.tsv \                   # 样本3结果
    > summary_card.tsv                   # 汇总矩阵

# ---- 更新数据库 ----
abricate-get_db --db card --force        # 强制更新CARD数据库

# ============================================================
# RGI: CARD数据库官方耐药基因注释工具
# 比ABRicate更详细，可注释突变型耐药
# ============================================================

# ---- 加载CARD数据库 ----
rgi load \
    --card_json /path/to/card.json \     # CARD数据库JSON文件
    --local                              # 加载到本地

# ---- 运行注释 ----
rgi main \
    --input_sequence contigs.fasta \     # 输入：组装序列或蛋白序列
    --output_file rgi_results \          # 输出文件前缀
    --input_type contig \                # 输入类型：contig / protein / read
    --local \                            # 使用本地数据库
    --clean \                            # 清理临时文件
    --num_threads 8                      # 线程数

# ---- 输出解读 ----
# rgi_results.txt 关键列：
#   Best_Hit_ARO       最佳匹配的ARO编号（抗性基因本体论）
#   Drug Class          耐药的药物类别
#   Resistance Mechanism 耐药机制
#   AMR Gene Family     耐药基因家族

# ============================================================
# KofamScan: 基于KOfam HMM的KEGG KO注释
# KEGG官方推荐，准确性高于BLAST方法
# ============================================================

# ---- 下载KOfam数据库 ----
wget ftp://ftp.genome.jp/pub/db/kofam/ko_list.gz
wget ftp://ftp.genome.jp/pub/db/kofam/profiles.tar.gz
gunzip ko_list.gz
tar xzf profiles.tar.gz

# ---- 运行注释 ----
exec_annotation \
    -f detail-tsv \                      # 输出格式：详细TSV
    -o kofamscan_results.tsv \           # 输出文件
    --cpu 16 \                           # CPU线程数
    --ko-list /path/to/ko_list \         # KO列表文件
    --profile /path/to/profiles/ \       # HMM profile目录
    --tmp-dir /tmp/kofamscan/ \          # 临时目录
    proteins.faa                         # 输入：蛋白质序列

# ---- 输出说明 ----
# 带 * 号的行表示该注释通过了阈值，可信度高
# 列：gene_name | KO | threshold | score | E-value | definition

KO (K00001) → Module (M00001) → Pathway (ko00010)
  基因级别        模块级别          通路级别

                    Gene Ontology
                    /      |      \
                  BP       MF      CC
              生物过程   分子功能  细胞组分
              "做什么"   "怎么做"  "在哪里"

示例（一个转录因子）：
  BP: transcription regulation（转录调控）    → 它参与的过程
  MF: DNA binding（DNA结合活性）             → 它的分子能力
  CC: nucleus（细胞核）                      → 它在细胞的位置

原始reads
    ↓
质控（fastp/Trimmomatic）
    ↓
去宿主（Bowtie2 vs 宿主基因组）
    ↓
组装（MEGAHIT/SPAdes）
    ↓
基因预测（Prodigal/MetaGeneMark）
    ↓ 得到蛋白序列
    ├─→ eggNOG-mapper → COG/GO/KEGG/Pfam 综合注释
    ├─→ KofamScan → KEGG KO精准注释
    ├─→ dbCAN3/run_dbcan v5 → CAZy碳水化合物活性酶注释
    ├─→ ABRicate/RGI → 耐药基因注释
    ├─→ ABRicate → 毒力因子注释
    └─→ InterProScan → Pfam/GO/IPR综合注释
         ↓
    丰度定量（reads mapping回基因 → TPM/RPKM）
         ↓
    差异分析（DESeq2/LEfSe）
         ↓
    富集分析（KEGG通路富集/GO富集）