注释对比：Prokka/Bakta/DRAM¶

# === 安装 ===
conda install -c bioconda prokka           # Prokka
conda install -c bioconda bakta            # Bakta
conda install -c bioconda dram             # DRAM（需要数据库配置）

# === 数据库 ===
# Bakta数据库下载（~30GB完整版，~2GB精简版）
bakta_db download --output bakta_db --type full     # 完整版
bakta_db download --output bakta_db --type light    # 精简版（推荐开始用）

# DRAM数据库配置（需要KEGG许可或使用免费替代）
DRAM-setup.py prepare_databases \
  --output_dir dram_db \
  --skip_uniref                            # 跳过UniRef（太大）

# === Prokka 注释 ===
prokka \
  genome.fa \                              # 输入基因组fasta
  --outdir prokka_output \                 # 输出目录
  --prefix my_genome \                     # 输出文件前缀
  --kingdom Bacteria \                     # 界: Bacteria/Archaea/Viruses
  --genus Escherichia \                    # 属名（可选，提高准确度）
  --species coli \                         # 种名（可选）
  --cpus 8 \                              # 线程数
  --mincontiglen 200 \                     # 最短contig
  --evalue 1e-06 \                         # E值阈值
  --rfam                                   # 使用Rfam注释ncRNA

# 输出文件:
# my_genome.gff    - GFF3格式注释（★最重要）
# my_genome.gbk    - GenBank格式
# my_genome.faa    - 蛋白质序列
# my_genome.ffn    - 核酸序列
# my_genome.fna    - 基因组序列（输入副本）
# my_genome.txt    - 注释统计摘要
# my_genome.tsv    - 基因信息表格

# === Bakta 注释（新项目首选） ===
bakta \
  genome.fa \                              # 输入基因组
  --db bakta_db/db \                      # 数据库路径
  --output bakta_output \                  # 输出目录
  --prefix my_genome \                     # 文件前缀
  --genus Escherichia \                    # 属名（可选）
  --species coli \                         # 种名（可选）
  --threads 8 \                            # 线程数
  --min-contig-length 200 \               # 最短contig
  --locus-tag ECOL                         # 基因位点标签前缀

# Bakta额外功能:
# --skip-trna/--skip-rrna                  # 跳过tRNA/rRNA预测
# --gram +/-/?                             # 革兰氏阳性/阴性/未知
# --keep-contig-headers                    # 保留原始contig名

# 输出文件（与Prokka类似）:
# my_genome.gff3   - GFF3格式注释
# my_genome.gbff   - GenBank格式
# my_genome.faa    - 蛋白质序列
# my_genome.json   - JSON格式结果
# my_genome.tsv    - 表格格式

# === Bakta vs Prokka 主要优势 ===
# 1. 数据库持续更新（Prokka依赖的数据库已老旧）
# 2. 符合INSDC提交标准（直接用于GenBank/ENA提交）
# 3. 包含AMR（抗性基因）和毒力因子注释
# 4. 蛋白质命名更规范（避免"hypothetical protein"过多）

# === DRAM 注释（宏基因组MAG功能分析） ===
# 第1步：注释
DRAM.py annotate \
  -i mags/*.fa \                           # 输入多个MAG基因组
  -o dram_output \                         # 输出目录
  --threads 8 \                            # 线程数
  --min_contig_size 2500                   # 最短contig

# 第2步：蒸馏（★DRAM特色功能）
DRAM.py distill \
  -i dram_output/annotations.tsv \        # 注释结果
  -o dram_distillate \                     # 蒸馏输出
  --trna_path dram_output/trnas.tsv \     # tRNA结果
  --rrna_path dram_output/rrnas.tsv       # rRNA结果

# 蒸馏输出（★面试重点）:
# product.html          - 代谢概览图（交互式网页）
# genome_stats.tsv      - 基因组统计
# metabolism_summary.xlsx - 代谢通路汇总表
# 功能模块覆盖度: 看每个MAG的代谢能力

# === 在R中分析注释结果 ===
# 读取Prokka/Bakta的TSV输出
genes <- read.table("my_genome.tsv", header=TRUE, sep="\t")

# 统计各类特征
table(genes$ftype)                          # CDS/tRNA/rRNA/等
cat("CDS总数:", sum(genes$ftype=="CDS"), "\n")
cat("假定蛋白:", sum(grepl("hypothetical", genes$product)), "\n")

# 计算注释率
annotated <- sum(!grepl("hypothetical", genes$product))
total_cds <- sum(genes$ftype=="CDS")
cat("注释率:", round(annotated/total_cds*100, 1), "%\n")
# 好的注释率: >60% | 一般: 40-60% | 差: <40%

# === 读取DRAM蒸馏结果 ===
library(readxl)
metabolism <- read_excel("metabolism_summary.xlsx")
# 看碳代谢、氮代谢、硫代谢等模块的覆盖度

# === 注释速查 ===
# Prokka（快速经典）
prokka genome.fa --outdir out --prefix name --kingdom Bacteria --cpus 8

# Bakta（★推荐新项目）
bakta genome.fa --db bakta_db/db --output out --prefix name --threads 8

# DRAM（宏基因组功能）
DRAM.py annotate -i mags/*.fa -o out --threads 8
DRAM.py distill -i out/annotations.tsv -o distillate

# 选择: 快速→Prokka | 新项目→Bakta | 功能/代谢→DRAM
# 评估: 注释率>60%好 | 假定蛋白<30%好 | tRNA≈20正常