617_基因组注释Prokka_Bakta¶

conda create -n annotation python=3.8  # 创建环境
conda activate annotation
conda install -c bioconda prokka  # 安装Prokka
prokka --version  # 查看版本

# === 最简单用法 ===
prokka genome.fasta  # 默认参数注释

# === 标准用法 ===
prokka \
    --kingdom Bacteria \       # 原核生物（Bacteria/Archaea）
    --genus Escherichia \      # 属名
    --species coli \           # 种名
    --strain K12 \             # 菌株名
    --outdir prokka_output \   # 输出目录
    --prefix ecoli_K12 \       # 输出文件前缀
    --locustag ECOLI \         # 基因座位标签
    --cpus 8 \                 # CPU数
    --rfam \                   # 启用Rfam搜索ncRNA（较慢但更全）
    --addgenes \               # 在GFF中添加基因特征
    --compliant \              # 生成GenBank兼容的输出
    genome.fasta               # 输入文件

# === 使用自定义蛋白数据库 ===
prokka \
    --proteins custom_db.fasta \  # 优先用这个数据库注释
    --outdir prokka_custom \
    --cpus 8 \
    genome.fasta

ls prokka_output/
# .gff   — GFF3格式注释文件（最常用，Roary等需要）
# .gbk   — GenBank格式（含序列+注释）
# .fna   — 核苷酸序列（与输入相同，但重命名了）
# .faa   — 蛋白序列
# .ffn   — 基因核苷酸序列
# .fsa   — GenBank提交格式的序列
# .tbl   — GenBank提交格式的特征表
# .sqn   — GenBank提交用的ASN.1格式
# .log   — 日志文件
# .txt   — 统计摘要
# .tsv   — 基因信息表格

# 查看统计摘要
cat prokka_output/ecoli_K12.txt
# 输出类似：
# organism: Escherichia coli K12
# contigs: 1
# bases: 4641652
# CDS: 4319
# rRNA: 22
# tRNA: 86
# tmRNA: 1

# === conda安装 ===
conda install -c bioconda bakta  # 安装Bakta

# === 下载数据库（必须！约53GB）===
bakta_db download --output /path/to/bakta_db  # 下载完整数据库

# 或下载轻量版数据库（约1.6GB，功能注释较少）
bakta_db download --output /path/to/bakta_db --type light

# 验证
bakta --version  # 查看版本
bakta_db info --db /path/to/bakta_db  # 查看数据库信息

# === 标准注释 ===
bakta \
    --db /path/to/bakta_db \   # 数据库路径
    --output bakta_output \    # 输出目录
    --prefix ecoli \           # 输出文件前缀
    --genus Escherichia \      # 属名
    --species coli \           # 种名
    --strain K12 \             # 菌株名
    --threads 8 \              # 线程数
    --locus-tag ECOLI \        # 基因座位标签
    --compliant \              # GenBank兼容输出
    genome.fasta               # 输入文件

# === 注释MAG（宏基因组组装基因组）===
bakta \
    --db /path/to/bakta_db \
    --output bakta_mag \
    --prefix mag_bin1 \
    --threads 8 \
    --meta \                   # MAG/宏基因组模式
    bin1.fasta

# === 使用Bakta Web（在线版，无需安装数据库）===
# 访问 https://bakta.computational.bio
# 上传FASTA文件即可

ls bakta_output/
# .gff3     — GFF3注释文件
# .gbff     — GenBank格式
# .fna      — 核苷酸序列
# .faa      — 蛋白序列
# .ffn      — 基因核苷酸序列
# .tsv      — 基因信息表格
# .json     — 机器可读的JSON格式（Bakta独有）
# .embl     — EMBL格式
# .txt      — 统计摘要
# .hypotheticals.tsv — 假定蛋白信息
# .hypotheticals.faa — 假定蛋白序列

#!/bin/bash
# 批量注释多个细菌基因组

GENOMES_DIR="genomes"    # 基因组文件目录
OUTPUT_DIR="annotations" # 输出目录
DB="/path/to/bakta_db"   # Bakta数据库路径
THREADS=8                # 线程数

mkdir -p ${OUTPUT_DIR}

# === 用Bakta批量注释 ===
for genome in ${GENOMES_DIR}/*.fasta; do
    name=$(basename "$genome" .fasta)  # 提取文件名
    echo "正在注释: ${name}"

    bakta \
        --db ${DB} \
        --output ${OUTPUT_DIR}/${name} \
        --prefix ${name} \
        --threads ${THREADS} \
        --locus-tag ${name} \
        ${genome}

    # 打印统计
    echo "--- ${name} 统计 ---"
    grep -E "CDS|tRNA|rRNA|sORF" ${OUTPUT_DIR}/${name}/${name}.txt
    echo ""
done

# === 汇总统计 ===
echo "========== 注释统计汇总 =========="
for dir in ${OUTPUT_DIR}/*/; do
    name=$(basename "$dir")
    cds=$(grep "CDS" ${dir}/${name}.txt | awk '{print $2}')
    hypo=$(grep -c "hypothetical" ${dir}/${name}.tsv 2>/dev/null || echo 0)
    echo "${name}: CDS=${cds}, 假定蛋白=${hypo}"
done

# === Python解析Bakta的JSON输出 ===
import json  # 导入json模块

# 读取JSON结果
with open("bakta_output/ecoli.json", 'r') as f:
    data = json.load(f)

# 统计基因类型
feature_types = {}
for feature in data.get("features", []):
    ftype = feature.get("type", "unknown")
    feature_types[ftype] = feature_types.get(ftype, 0) + 1

print("基因类型统计:")
for ftype, count in sorted(feature_types.items(), key=lambda x: -x[1]):
    print(f"  {ftype}: {count}")

# 提取所有有功能注释的基因
annotated = [f for f in data["features"] 
             if f.get("type") == "cds" and "hypothetical" not in f.get("product", "").lower()]
print(f"\n有功能注释的CDS: {len(annotated)}/{len([f for f in data['features'] if f.get('type')=='cds'])}")