基因组注释基础（Genome Annotation）¶

原始测序数据 → 质控 → 组装(Assembly) → 【基因组注释】 → 下游分析
                                            ↓
                                    结构注释（找基因）
                                            ↓
                                    功能注释（查功能）
                                            ↓
                                    富集分析 / 比较基因组学 / 代谢重建

# ========== Prodigal 安装 ==========
conda install -c bioconda prodigal  # conda 安装

# ========== 基本用法 ==========
# -i: 输入基因组 fasta 文件
# -o: 输出 GFF 格式的基因坐标文件
# -a: 输出预测的蛋白质氨基酸序列（用于后续功能注释）
# -d: 输出预测的核酸序列
# -p: 运行模式，single=单个基因组，meta=宏基因组
prodigal \
  -i genome.fasta \        # 输入：组装好的基因组
  -o genes.gff \           # 输出：基因坐标（GFF3格式）
  -a proteins.faa \        # 输出：蛋白质序列（功能注释用）
  -d genes.fna \           # 输出：基因核酸序列
  -p single \              # 模式：single（单基因组）或 meta（宏基因组）
  -f gff                   # 输出格式：gff

# 宏基因组模式（不能自学习，用通用模型）
prodigal -i contigs.fasta -o meta_genes.gff -a meta_proteins.faa -p meta -f gff

# GeneMark 需要先从官网申请 license key
# 原核基因组用法
gmhmmp -m GeneMark_hmm_combined.mod -o genes.lst genome.fasta

# 真核基因组自训练模式（GeneMark-ES）
gmes_petap.pl --ES --sequence genome.fasta

# 安装
conda install -c bioconda augustus

# 查看已有物种模型
augustus --species=help

# 基本用法：指定最接近的物种模型
# --species: 选择最接近的已训练物种模型
# --gff3: 输出 GFF3 格式
augustus --species=human --gff3=on genome.fasta > augustus_genes.gff3

# 整合 RNA-seq 提示信息（hints）可以大幅提高准确率
augustus --species=human --hintsfile=hints.gff --extrinsicCfgFile=extrinsic.cfg genome.fasta > genes.gff3

# 安装（依赖多，推荐 conda）
conda install -c bioconda maker

# MAKER 三步走：
# 1. 生成控制文件
maker -CTL
# 会生成三个文件：
# maker_opts.ctl  → 主要设置（输入文件、物种模型等）
# maker_bopts.ctl → BLAST 参数
# maker_exe.ctl   → 外部程序路径

# 2. 编辑 maker_opts.ctl，设置输入文件
# genome=genome.fasta        # 基因组文件
# est=transcripts.fasta      # 转录组证据（EST/RNA-seq组装）
# protein=proteins.fasta     # 同源蛋白证据
# augustus_species=human      # Augustus 物种模型

# 3. 运行 MAKER
maker maker_opts.ctl maker_bopts.ctl maker_exe.ctl

你要注释什么类型的基因组？
├── 原核生物（细菌/古菌）
│   ├── 单个基因组 → Prodigal (-p single) 或 Prokka（一站式）
│   └── 宏基因组 contigs → Prodigal (-p meta)
├── 真核生物
│   ├── 有 RNA-seq 数据 → BRAKER2 或 MAKER + Augustus + hints
│   ├── 无 RNA-seq → Augustus（选近缘物种模型）+ GeneMark-ES
│   └── 复杂真核（大基因组）→ MAKER pipeline
└── 病毒 → Prodigal-gv 或 PHANOTATE

# ========== 安装 ==========
conda install -c bioconda eggnog-mapper

# ========== 下载数据库（首次运行必须，约50GB）==========
download_eggnog_data.py --data_dir /path/to/eggnog_data

# ========== 运行注释 ==========
# -i: 输入蛋白质序列（来自 Prodigal 的 -a 输出）
# --output: 输出文件前缀
# --data_dir: eggNOG 数据库路径
# -m: 搜索模式，diamond（快速）或 mmseqs（更快）
# --cpu: 使用的 CPU 核心数
# --go_evidence: GO 注释的证据级别筛选
# --tax_scope: 分类范围限制（如 Bacteria）
emapper.py \
  -i proteins.faa \                 # 输入：蛋白质序列（来自Prodigal输出的.faa）
  --output t2d_annotation \         # 输出文件前缀
  --data_dir /path/to/eggnog_data \ # 数据库路径
  -m diamond \                      # 搜索引擎：diamond（推荐，速度快）
  --cpu 16 \                        # 使用16个CPU核心
  --go_evidence non-electronic \    # GO注释只保留有实验证据的
  --tax_scope Bacteria \            # 限制搜索范围为细菌（加速+减少假阳性）
  --override                        # 覆盖已有输出文件

# 输出文件解读（.emapper.annotations）：
# 第1列: 查询序列ID
# 第5列: eggNOG 直系同源组（OG）
# 第7列: COG 功能分类字母
# 第10列: GO terms（GO:0006810 等）
# 第12列: KEGG KO 编号（ko:K00001 等）
# 第13列: KEGG Pathway
# 第21列: Pfam 结构域
# 第22列: CAZy 家族

# 安装（较大，约 15GB）
# 从 InterPro 官网下载解压即可
wget https://ftp.ebi.ac.uk/pub/software/unix/iprscan/5/5.77-108.0/interproscan-5.77-108.0-64-bit.tar.gz
tar -xzf interproscan-5.77-108.0-64-bit.tar.gz

# 运行注释
# -i: 输入蛋白质序列
# -f: 输出格式（tsv/xml/json/gff3）
# -goterms: 同时输出 GO 注释
# -pa: 同时输出 Pathway 注释
interproscan.sh \
  -i proteins.faa \    # 输入蛋白质序列
  -f tsv,gff3 \        # 输出格式：tsv表格 + gff3坐标
  -goterms \           # 输出GO注释
  -pa \                # 输出Pathway注释
  -cpu 16 \            # 使用16个CPU
  -dp                  # 关闭预计算查询（离线模式）

# ========== 安装 ==========
conda install -c bioconda prokka

# ========== 注释一个细菌基因组 ==========
# --outdir: 输出目录
# --prefix: 输出文件名前缀
# --kingdom: 生物界（Bacteria/Archaea/Viruses）
# --genus/--species: 物种分类信息
# --cpus: CPU核心数
# --locustag: 基因标签前缀（给每个基因编号用的）
# --compliant: 输出符合NCBI提交标准的格式
prokka \
  --outdir prokka_output \      # 输出目录
  --prefix my_genome \          # 文件名前缀
  --kingdom Bacteria \          # 指定为细菌
  --genus Escherichia \         # 属名（可选但推荐）
  --species coli \              # 种名（可选但推荐）
  --cpus 16 \                   # 使用16个CPU核心
  --locustag ECOLI \            # 基因编号前缀，如 ECOLI_00001
  --compliant \                 # 生成NCBI提交兼容格式
  --addgenes \                  # 在GFF中添加gene特征
  genome.fasta                  # 输入：组装好的基因组fasta

# Prokka 输出文件一览：
# .gff   → GFF3 注释文件（结构+功能信息）
# .gbk   → GenBank 格式（可用Artemis查看）
# .faa   → 蛋白质氨基酸序列
# .ffn   → 基因核酸序列
# .fna   → 基因组核酸序列（带注释头信息）
# .tbl   → NCBI提交格式的特征表
# .tsv   → 注释汇总表（基因ID、产物名、EC号等）
# .txt   → 统计摘要（总基因数、CDS数、tRNA数、rRNA数）
# .log   → 运行日志

# 安装
conda install -c bioconda dram

# 准备数据库（首次运行）
DRAM-setup.py prepare_databases --output_dir /path/to/dram_db

# 注释宏基因组MAGs
# -i: 输入基因组fasta
# -o: 输出目录
DRAM.py annotate \
  -i genome.fasta \          # 输入：MAG或基因组
  -o dram_annotation \       # 输出目录
  --threads 16               # CPU核心数

# 生成代谢通路可视化
DRAM.py distill \
  -i dram_annotation/annotations.tsv \  # 注释结果
  -o dram_distillate                    # 可视化输出目录
# 会生成 metabolism_summary.xlsx 和 product.html（交互式代谢图）

组装后的基因组 (contigs.fasta)
    │
    ├─ 方案A：一站式 → Prokka → 全套输出
    │
    └─ 方案B：分步式
         ├── Step 1: Prodigal 预测基因 → genes.gff + proteins.faa
         ├── Step 2: eggNOG-mapper 功能注释 → GO/KEGG/COG
         ├── Step 3: InterProScan 结构域注释 → Pfam/TIGRFAM
         └── Step 4: 专用数据库注释（可选）
              ├── dbCAN3/run_dbcan v5 → CAZy碳水化合物酶
              ├── ABRicate → 耐药基因(CARD)、毒力因子(VFDB)
              └── SignalP/TMHMM → 信号肽/跨膜区

组装后的基因组 (genome.fasta)
    │
    ├── Step 0: RepeatMasker → 屏蔽重复序列（先把"噪音"去掉）
    │
    ├── Step 1: RNA-seq 比对 → HISAT2/STAR → 转录组证据
    │
    ├── Step 2: 基因预测（多种方法综合）
    │    ├── Ab initio: Augustus, GeneMark-ES
    │    ├── 同源证据: GeMoMa, Exonerate
    │    └── 转录组证据: StringTie, PASA
    │
    ├── Step 3: 证据整合 → MAKER/EvidenceModeler → 共识基因模型
    │
    ├── Step 4: 功能注释
    │    ├── BLAST vs UniProt/Nr → 同源功能推断
    │    ├── InterProScan → 结构域和GO注释
    │    └── eggNOG-mapper → KEGG/COG注释
    │
    └── Step 5: 手动审查 → Apollo/JBrowse → 修正错误预测

蛋白质序列 → eggNOG-mapper → KO 编号 → Pathway 映射
                                ↓
                        K00844 (HK, hexokinase)
                                ↓
                        ko00010 (Glycolysis/Gluconeogenesis)
                                ↓
                        用于差异代谢通路分析

# 使用 HMMER 直接搜索 Pfam
# -E: E-value 阈值（越小越严格）
# --domtblout: 输出每个结构域的详细比对结果
hmmscan \
  --cpu 8 \                    # CPU 核心数
  -E 1e-5 \                    # E-value 阈值
  --domtblout pfam_results.txt \  # 结构域结果输出
  /path/to/Pfam-A.hmm \       # Pfam HMM 数据库
  proteins.faa                 # 输入蛋白质序列

超几何分布公式：

P(X = k) = C(K,k) * C(N-K, n-k) / C(N,n)

其中：
N = 总基因数（全班100人）
K = 某GO term包含的基因数（全班30个男生）
n = 差异基因数（20个好学生）
k = 差异基因中属于该GO term的数量（15个男生好学生）
C(a,b) = 组合数，从a中选b个的方式数

白话：从100人里随机抽20人，恰好抽到>=15个男生的概率是多少？
如果概率很低（p < 0.05），说明"好学生里男生多"不是偶然的 → 富集了！

1. 确定背景基因集（N）：通常是所有被注释到的基因
2. 确定目标基因集（n）：差异表达基因或特定条件的基因
3. 对每个GO/KEGG term：
   - 计算在背景中有多少基因属于该term（K）
   - 计算在目标集中有多少基因属于该term（k）
   - 用超几何检验计算p值
4. 多重检验校正（BH/FDR）：因为同时测试了成百上千个term
5. 筛选显著富集的term（通常 FDR < 0.05）

# ========== clusterProfiler GO 富集分析 ==========
# 安装（如果没有的话）
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("clusterProfiler", "org.Hs.eg.db", "enrichplot"))

# 加载包
library(clusterProfiler)  # 富集分析核心包
library(enrichplot)        # 可视化包
library(org.Hs.eg.db)     # 人类基因注释数据库（肠道菌群宿主侧分析用）
# 注意：如果分析细菌基因，不用 org.Hs.eg.db，而是用自己的注释表

# ---------- 方法1：使用自定义注释（宏基因组常用）----------
# 读取 eggNOG-mapper 输出的注释文件
anno <- read.delim("t2d_annotation.emapper.annotations", 
                    comment.char = "#",    # 跳过#开头的注释行
                    header = FALSE)        # eggNOG输出无表头

# 构建 gene2go 映射表（term2gene格式）
# 假设第1列是基因ID，第10列是GO terms
gene2go <- data.frame(
  gene = anno$V1,       # 基因ID
  go = anno$V10         # GO注释（可能有多个，用逗号分隔）
)

# 展开多值GO注释（一个基因可能有多个GO term）
library(tidyr)
gene2go_expanded <- gene2go %>%
  separate_rows(go, sep = ",") %>%    # 把 "GO:0001,GO:0002" 拆成两行
  filter(go != "-") %>%               # 去掉没有注释的("-")
  filter(go != "")                    # 去掉空值

# 转换为 enricher 需要的 TERM2GENE 格式
term2gene <- gene2go_expanded[, c("go", "gene")]  # 注意顺序：先term再gene

# 定义差异基因列表（从差异分析结果中获得）
diff_genes <- c("gene_001", "gene_002", "gene_003")  # 替换为实际差异基因ID

# 运行 GO 富集分析（自定义注释版）
ego <- enricher(
  gene = diff_genes,        # 差异基因列表
  TERM2GENE = term2gene,    # 自定义的GO注释映射
  pvalueCutoff = 0.05,      # p值阈值
  pAdjustMethod = "BH",     # 多重检验校正方法（Benjamini-Hochberg）
  qvalueCutoff = 0.2        # q值（FDR）阈值
)

# ---------- 方法2：使用内置数据库（模式生物用）----------
# 如果分析的是人类基因（比如T2D宿主转录组）
ego2 <- enrichGO(
  gene = diff_gene_ids,       # 差异基因的Entrez ID
  OrgDb = org.Hs.eg.db,       # 人类注释数据库
  ont = "BP",                  # 选择GO分支：BP/MF/CC/ALL
  pAdjustMethod = "BH",       # BH法校正p值
  pvalueCutoff = 0.05,        # p值阈值
  qvalueCutoff = 0.2          # q值阈值
)

# ---------- 可视化 ----------
# 气泡图（最常用）
dotplot(ego, showCategory = 20) +   # 显示top20个GO term
  ggtitle("GO Enrichment Analysis - T2D Gut Microbiome")

# 柱状图
barplot(ego, showCategory = 15)

# 网络图（展示GO term之间的关系）
emapplot(ego, showCategory = 20)

# KEGG 富集分析（类似逻辑）
# 构建 term2gene 用 KEGG pathway
kegg_term2gene <- data.frame(
  pathway = anno$V13,    # KEGG pathway列
  gene = anno$V1         # 基因ID
)
kegg_term2gene <- kegg_term2gene %>%
  separate_rows(pathway, sep = ",") %>%
  filter(pathway != "-" & pathway != "")

ekegg <- enricher(
  gene = diff_genes,
  TERM2GENE = kegg_term2gene,
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH"
)

dotplot(ekegg, showCategory = 15) +
  ggtitle("KEGG Pathway Enrichment - T2D vs Control")

宏基因组测序数据
    ↓
质控 → 组装（MEGAHIT）→ contigs
    ↓
Prodigal -p meta → 预测基因（蛋白质序列）     ← 结构注释
    ↓
eggNOG-mapper → GO/KEGG/COG 功能注释          ← 功能注释
    ↓
差异分析（T2D组 vs 健康对照组）
    ↓
clusterProfiler 富集分析                       ← 功能富集
    ↓
发现 T2D 相关的差异代谢通路
    例如：短链脂肪酸(SCFA)合成通路 ↓
         丁酸合成基因(bcd, buk)在T2D中显著降低
         LPS合成通路在T2D中富集（促炎）

# 原核一站式注释
prokka --kingdom Bacteria --cpus 16 genome.fasta --outdir anno_out

# 宏基因组基因预测
prodigal -i contigs.fasta -o genes.gff -a proteins.faa -p meta -f gff

# eggNOG-mapper 功能注释
emapper.py -i proteins.faa --output result -m diamond --cpu 16

# InterProScan 结构域注释
interproscan.sh -i proteins.faa -f tsv,gff3 -goterms -pa -cpu 16

# DRAM 代谢重建
DRAM.py annotate -i genome.fasta -o dram_out --threads 16
DRAM.py distill -i dram_out/annotations.tsv -o dram_distill