173_宏基因组病毒鉴定geNomad¶

# 方法1：conda安装（推荐）
conda create -n genomad -c conda-forge -c bioconda genomad
conda activate genomad

# 方法2：pip安装
pip install genomad

# 下载数据库
genomad download-database ./genomad_db

# 验证安装
genomad --version
genomad --help

# 组装宏基因组（如尚未组装）
megahit -1 sample_R1.fq.gz -2 sample_R2.fq.gz -o megahit_out -t 16

# 过滤短序列（保留≥1000bp）
seqkit seq -m 1000 megahit_out/final.contigs.fa > contigs_1k.fa

# 查看序列统计
seqkit stats contigs_1k.fa

# 一步完成全部分析
genomad end-to-end \
    contigs_1k.fa \
    genomad_output \
    genomad_db/genomad_db \
    --threads 16 \
    --cleanup \
    --splits 8

# 参数说明：
# --threads: 使用线程数
# --cleanup: 完成后清理中间文件
# --splits: 将数据库分成多份降低内存使用

# 如果只想找病毒（不找质粒）
genomad end-to-end \
    contigs_1k.fa \
    genomad_output \
    genomad_db/genomad_db \
    --threads 16 \
    --disable-find-proviruses  # 跳过前噬菌体检测（加速）

# 主要输出文件结构
genomad_output/
├── contigs_1k_annotate/          # 基因注释结果
│   ├── contigs_1k_genes.tsv      # 预测的基因
│   ├── contigs_1k_taxonomy.tsv   # 分类信息
│   └── contigs_1k_mmseqs2/       # MMseqs2比对结果
├── contigs_1k_find_proviruses/   # 前噬菌体检测
│   ├── contigs_1k_provirus.tsv   # 前噬菌体列表
│   └── contigs_1k_provirus.fna   # 前噬菌体序列
├── contigs_1k_summary/           # 汇总结果（最重要）
│   ├── contigs_1k_virus_summary.tsv      # 病毒汇总
│   ├── contigs_1k_plasmid_summary.tsv    # 质粒汇总
│   ├── contigs_1k_virus.fna              # 病毒序列
│   └── contigs_1k_plasmid.fna            # 质粒序列
└── contigs_1k_aggregated_classification/  # 整合分类

# 查看病毒汇总表
head -5 genomad_output/contigs_1k_summary/contigs_1k_virus_summary.tsv
# 关键列：seq_name, length, topology, coordinates, n_genes,
#         genetic_code, virus_score, fdr, taxonomy

# 用Python过滤高置信病毒
python3 << 'EOF'
import pandas as pd

virus = pd.read_csv(
    "genomad_output/contigs_1k_summary/contigs_1k_virus_summary.tsv",
    sep="\t"
)
# 筛选条件：virus_score >= 0.8 且 FDR <= 0.05
high_conf = virus[(virus['virus_score'] >= 0.8) & (virus['fdr'] <= 0.05)]
print(f"高置信病毒序列数: {len(high_conf)}")
high_conf.to_csv("high_confidence_viruses.tsv", sep="\t", index=False)

# 统计分类分布
print(virus['taxonomy'].value_counts().head(10))
EOF

# 用CheckV评估病毒基因组质量
checkv end_to_end \
    genomad_output/contigs_1k_summary/contigs_1k_virus.fna \
    checkv_output \
    -d /path/to/checkv-db \
    -t 16

# 查看CheckV质量评估结果
cut -f1,6,7,8 checkv_output/quality_summary.tsv | head

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 读取结果
virus = pd.read_csv(
    "genomad_output/contigs_1k_summary/contigs_1k_virus_summary.tsv",
    sep="\t"
)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 1. 得分分布
axes[0].hist(virus['virus_score'], bins=30, edgecolor='black', alpha=0.7)
axes[0].set_xlabel('Virus Score')
axes[0].set_ylabel('Count')
axes[0].set_title('Score Distribution')
axes[0].axvline(x=0.7, color='red', linestyle='--', label='Threshold')
axes[0].legend()

# 2. 序列长度分布
axes[1].hist(virus['length'], bins=50, edgecolor='black', alpha=0.7, color='green')
axes[1].set_xlabel('Length (bp)')
axes[1].set_ylabel('Count')
axes[1].set_title('Virus Contig Length')
axes[1].set_xscale('log')

# 3. 拓扑类型
topology_counts = virus['topology'].value_counts()
axes[2].pie(topology_counts, labels=topology_counts.index, autopct='%1.1f%%')
axes[2].set_title('Genome Topology')

plt.tight_layout()
plt.savefig('genomad_summary.png', dpi=300, bbox_inches='tight')
plt.show()

# 完整流程一键运行
genomad end-to-end contigs.fa output genomad_db --threads 16 --cleanup

# 仅注释（不分类）
genomad annotate contigs.fa output genomad_db --threads 16

# 单独运行前噬菌体检测
genomad find-proviruses contigs.fa output genomad_db --threads 16

# 自定义得分阈值
genomad score-calibration contigs.fa output genomad_db --min-score 0.8

# 配合CheckV使用
checkv end_to_end output/contigs_virus.fna checkv_out -d checkv_db -t 16

# 提取特定分类的病毒
grep "Caudoviricetes" output/contigs_virus_summary.tsv > caudovirales.tsv

# 批量处理多个样本
for sample in sample1 sample2 sample3; do
    genomad end-to-end ${sample}_contigs.fa ${sample}_genomad genomad_db -t 16
done