605 噬菌体宏基因组分析（Phage Viromics）¶

宏基因组组装(MEGAHIT/SPAdes)
    → 病毒序列鉴定(VirSorter2 + VIBRANT + geNomad)
    → 质量控制(CheckV)
    → 共识过滤（多工具交集）
    → 去冗余(CD-HIT/vRhyme)
    → 宿主预测(iPHoP)
    → 功能注释(DRAM-v)
    → 分类学分类(vConTACT2/geNomad)
    → 定量(CoverM)

# ========== 步骤1：组装（获取contigs） ==========
# 病毒鉴定需要较长的contig，建议组装后过滤

# 用MEGAHIT组装
megahit \
    -1 clean_R1.fastq.gz \        # 输入：质控后正向读段
    -2 clean_R2.fastq.gz \        # 输入：质控后反向读段
    -o megahit_output/ \          # 输出目录
    -t 32 \                       # 32个线程
    --min-contig-len 1000 \       # 最短contig 1000bp（病毒基因组一般>3kb）
    --k-min 21 \                  # 最小k-mer值
    --k-max 141 \                 # 最大k-mer值
    --k-step 12                   # k-mer步长

# 过滤短contig（病毒鉴定建议>5kb）
# 用seqkit过滤
seqkit seq \
    megahit_output/final.contigs.fa \  # 输入：组装结果
    -m 5000 \                           # 最短5000bp
    -o contigs_5kb.fa                   # 输出：过滤后的contigs

echo ">=5kb的contig数:" && grep -c ">" contigs_5kb.fa  # 统计数量

# ========== 步骤2：VirSorter2 — 病毒序列鉴定 ==========
# VirSorter2用机器学习模型，支持dsDNA噬菌体、ssDNA病毒、RNA病毒等

# 安装VirSorter2
conda create -n vs2 -c conda-forge -c bioconda virsorter=2  # 创建conda环境
conda activate vs2  # 激活环境

# 下载数据库（首次使用需要）
virsorter setup \
    -d /db/virsorter2/ \  # 数据库安装路径
    -j 8                  # 8个线程下载

# 运行VirSorter2
virsorter run \
    -w vs2_output/ \                    # 输出目录
    -i contigs_5kb.fa \                 # 输入：>=5kb的contigs
    -d /db/virsorter2/ \                # 数据库路径
    --include-groups dsDNAphage,ssDNA \ # 鉴定的病毒类型
    --min-length 5000 \                 # 最短序列5000bp
    --min-score 0.5 \                   # 最低分数阈值（默认0.5，高灵敏度）
    --keep-original-seq \               # 保留原始序列（不切割前噬菌体）
    -j 16 \                             # 16个线程
    --provirus-off                      # 暂时关闭前噬菌体检测（后续用CheckV做）

# 输出解读：
# vs2_output/final-viral-score.tsv  — 每条序列的病毒分数和分组
# vs2_output/final-viral-combined.fa — 鉴定为病毒的序列

# 查看结果分类
# max_score: 病毒分数（0-1，越高越可能是病毒）
# max_score_group: 最可能的病毒类型
head vs2_output/final-viral-score.tsv

# ========== 步骤3：VIBRANT — 第二轮病毒鉴定 ==========
# VIBRANT基于蛋白注释（KEGG/VOG/Pfam）来鉴定病毒
# 同时提供功能注释和基因组质量评估

# 安装VIBRANT
conda install -c bioconda vibrant  # 安装

# 下载数据库
download-db.sh  # VIBRANT自带的数据库下载脚本

# 运行VIBRANT
VIBRANT_run.py \
    -i contigs_5kb.fa \              # 输入：>=5kb的contigs
    -folder vibrant_output/ \        # 输出目录
    -t 16 \                          # 16个线程
    -l 5000 \                        # 最短序列5000bp
    -virome                          # 如果样本是纯virome测序加这个参数
                                     # 如果是普通宏基因组，去掉-virome

# VIBRANT输出结构：
# vibrant_output/VIBRANT_phages_*/    — 鉴定为噬菌体的序列
# vibrant_output/VIBRANT_results_*/   — 详细结果（注释、质量等）

# 查看鉴定结果
echo "VIBRANT鉴定的噬菌体数:"
grep -c ">" vibrant_output/VIBRANT_phages_contigs_5kb/contigs_5kb.phages_combined.fna

# ========== 步骤4：CheckV — 病毒基因组质量评估 ==========
# CheckV做两件核心事情：
# 1. 评估基因组完整性（完整、高质量、中质量、低质量）
# 2. 去除宿主污染区域（trim掉细菌序列）

# 安装CheckV
conda install -c bioconda checkv  # 安装

# 下载数据库
checkv download_database /db/checkv/  # 下载数据库

# 先合并VirSorter2和VIBRANT的结果
cat vs2_output/final-viral-combined.fa \
    vibrant_output/VIBRANT_phages_*/contigs_5kb.phages_combined.fna \
    > all_viral_candidates.fa  # 合并所有候选病毒序列

# 去重复序列名
seqkit rmdup all_viral_candidates.fa \    # 去除重复序列
    -s \                                   # 按序列去重
    -o viral_candidates_dedup.fa           # 输出

# 运行CheckV
checkv end_to_end \
    viral_candidates_dedup.fa \    # 输入：候选病毒序列
    checkv_output/ \               # 输出目录
    -d /db/checkv/ \               # 数据库路径
    -t 16                          # 16个线程

# CheckV输出解读：
# quality_summary.tsv — 核心结果，每条序列的质量等级
# proviruses.fna      — 前噬菌体序列（去除了宿主区域）
# viruses.fna         — 纯病毒序列

# 查看质量分布
cut -f8 checkv_output/quality_summary.tsv | sort | uniq -c | sort -rn
# 质量等级：Complete > High-quality > Medium-quality > Low-quality > Not-determined

# ========== 步骤5：共识过滤 — 只保留多工具都认可的病毒 ==========
import pandas as pd  # 导入pandas

# 读取各工具的鉴定结果
vs2_ids = set()  # VirSorter2鉴定的序列ID集合
with open("vs2_output/final-viral-combined.fa") as f:
    for line in f:
        if line.startswith(">"):
            vs2_ids.add(line.strip().split()[0][1:])  # 提取序列ID

vibrant_ids = set()  # VIBRANT鉴定的序列ID集合
with open("vibrant_output/VIBRANT_phages_contigs_5kb/contigs_5kb.phages_combined.fna") as f:
    for line in f:
        if line.startswith(">"):
            vibrant_ids.add(line.strip().split()[0][1:])

# 读取CheckV质量评估结果
checkv = pd.read_csv("checkv_output/quality_summary.tsv", sep="\t")  # 读取CheckV结果

# 策略1：严格共识（VirSorter2 AND VIBRANT 都鉴定）
consensus_strict = vs2_ids & vibrant_ids  # 取交集
print(f"严格共识: {len(consensus_strict)} 条病毒序列")

# 策略2：基于CheckV质量过滤
# 保留中质量以上 + 至少有1个病毒基因的序列
quality_pass = checkv[
    (checkv["checkv_quality"].isin(["Complete", "High-quality", "Medium-quality"])) |  # 中等质量以上
    (checkv["viral_genes"] >= 1)  # 或至少有1个病毒基因
]["contig_id"].tolist()

print(f"CheckV质量过滤: {len(quality_pass)} 条病毒序列")

# 最终保留：共识 + 质量过滤
final_viruses = set(quality_pass) & (vs2_ids | vibrant_ids)  # 至少一个工具鉴定 + 质量通过
print(f"最终保留: {len(final_viruses)} 条高可信度病毒序列")

# 提取最终的病毒序列
with open("final_viral_sequences.txt", "w") as f:  # 写入最终序列ID列表
    for vid in final_viruses:
        f.write(vid + "\n")

# 用seqkit提取最终病毒序列
seqkit grep \
    -f final_viral_sequences.txt \        # 按ID列表筛选
    checkv_output/combined.fna \          # 从CheckV的输出中提取
    -o final_viruses.fa                   # 最终的病毒序列文件

echo "最终病毒序列数:" && grep -c ">" final_viruses.fa

# ========== 步骤6：去冗余（Dereplication） ==========
# 如果有多个样本，需要将所有样本的病毒序列去冗余

# 方法1：CD-HIT去冗余（简单直接）
cd-hit-est \
    -i final_viruses.fa \       # 输入：所有样本合并的病毒序列
    -o viral_clusters.fa \      # 输出：去冗余后的代表序列
    -c 0.95 \                   # 95%序列一致性阈值（同一种病毒）
    -aS 0.85 \                  # 较短序列的85%被覆盖
    -n 10 \                     # 词长
    -M 16000 \                  # 内存16G
    -T 16                       # 16线程

# 方法2：vRhyme分箱（将碎片化的病毒基因组拼回完整基因组）
vRhyme \
    -i final_viruses.fa \           # 输入：病毒序列
    -o vrhyme_output/ \             # 输出目录
    -b mapped_reads.bam \           # BAM比对文件（用于覆盖度信息）
    -t 16                           # 16线程

# ========== 步骤7：宿主预测（这个噬菌体感染谁？） ==========
# iPHoP整合多种方法预测病毒宿主

iphop predict \
    --fa_file viral_clusters.fa \      # 输入：病毒序列
    --out_folder iphop_output/ \       # 输出目录
    --db_dir /db/iphop/ \              # iPHoP数据库
    --num_threads 16                   # 16线程

# iPHoP输出中最重要的文件：
# Host_prediction_to_genus_m90.csv — 属水平的宿主预测

# ========== 步骤8：功能注释（DRAM-v） ==========
# DRAM-v专门对病毒基因组做注释，特别关注辅助代谢基因(AMG)

DRAM-v.py annotate \
    -i viral_clusters.fa \                   # 输入：病毒序列
    -v vs2_output/final-viral-score.tsv \    # VirSorter2的分数文件
    -o dramv_output/ \                       # 输出目录
    --threads 16                             # 16线程

# 生成DRAM-v汇总图
DRAM-v.py distill \
    -i dramv_output/annotations.tsv \    # 输入：注释结果
    -o dramv_distill/                    # 输出目录

# DRAM-v特别关注的AMG（辅助代谢基因）：
# 这些基因被噬菌体"偷"来帮助宿主提高代谢效率
# 常见AMG：光合作用、碳固定、硫代谢、核苷酸合成相关基因

# ========== 步骤9：定量（用CoverM计算覆盖度） ==========

# 将reads比对到病毒序列
bowtie2-build viral_clusters.fa viral_index  # 建立索引

bowtie2 \
    -x viral_index \                  # 索引
    -1 clean_R1.fastq.gz \           # 正向读段
    -2 clean_R2.fastq.gz \           # 反向读段
    --threads 16 \                   # 16线程
    --no-unal \                      # 不输出未比对的reads
    | samtools sort -@ 8 -o mapped.bam  # 排序并输出BAM

samtools index mapped.bam  # 建立BAM索引

# 用CoverM计算RPKM
coverm genome \
    -b mapped.bam \                    # 输入：BAM文件
    --genome-fasta-files viral_clusters.fa \  # 病毒基因组
    -m rpkm \                          # 计算RPKM
    -o viral_abundance.tsv \           # 输出：丰度表
    --min-covered-fraction 0           # 最小覆盖比例

# ========== geNomad：一站式病毒+质粒鉴定 ==========
# 2024年Nature Biotechnology发表，越来越多流程用它替代VirSorter2

# 安装
conda install -c bioconda genomad  # 安装geNomad

# 下载数据库
genomad download-database /db/genomad/  # 下载

# 运行geNomad
genomad end-to-end \
    contigs_5kb.fa \              # 输入：contigs
    genomad_output/ \             # 输出目录
    /db/genomad/genomad_db \      # 数据库路径
    --threads 16 \                # 16线程
    --cleanup                     # 清理临时文件

# geNomad优势：
# 1. 同时鉴定病毒和质粒
# 2. 自带分类学分类
# 3. 集成CRF模型，前噬菌体边界检测更准
# 4. 速度比VirSorter2快

========== 噬菌体宏基因组标准流程 ==========
组装(MEGAHIT) → 鉴定(VirSorter2+VIBRANT+geNomad)
    → 质控(CheckV) → 共识过滤 → 去冗余(CD-HIT)
    → 宿主预测(iPHoP) → 功能注释(DRAM-v)
    → 分类(vConTACT2/geNomad) → 定量(CoverM)

========== 核心工具矩阵 ==========
病毒鉴定 ：VirSorter2 / VIBRANT / geNomad / DeepVirFinder
质量评估 ：CheckV（必用）
宿主预测 ：iPHoP / CHERRY / WIsH
功能注释 ：DRAM-v / Pharokka
分类学   ：vConTACT2 / geNomad / PhaGCN
去冗余   ：CD-HIT-EST / dRep
分箱     ：vRhyme
定量     ：CoverM / Bowtie2+samtools
整合流程 ：MVP / ViWrap / ViroProfiler

========== 病毒基因组质量等级（CheckV） ==========
Complete       ：完整基因组（100%完整度）
High-quality   ：>90%完整度
Medium-quality ：50-90%完整度
Low-quality    ：<50%完整度
Not-determined ：无法判断

========== 推荐共识策略 ==========
高可信度：VirSorter2 ∩ VIBRANT ∩ CheckV(medium+)
中可信度：(VirSorter2 ∪ VIBRANT) ∩ CheckV(至少1个viral gene)
探索性  ：任一工具鉴定 + CheckV不含host gene