672 噬菌体-宿主关系预测¶

# iPHoP安装（Integrated Phage Host Prediction）
# 2023年Nature Protocols发表，整合12种方法
conda create -n iphop -c conda-forge -c bioconda iphop  # 创建环境
conda activate iphop  # 激活环境

# 下载数据库（约100GB）
iphop download \
  --db_dir /path/to/iphop_db/ \  # 数据库目录
  --no_prompt                     # 不需要确认

# 运行iPHoP预测
iphop predict \
  --fa_file viral_contigs.fasta \  # 输入噬菌体contigs
  --db_dir /path/to/iphop_db/ \   # 数据库路径
  --out_dir iphop_results/ \       # 输出目录
  --num_threads 16                 # 线程数

# 输出文件说明：
# iphop_results/
# ├── Host_prediction_to_genus_m90.csv    # 属级预测（90%置信度）
# ├── Host_prediction_to_genome_m90.csv   # 基因组级预测
# ├── Detailed_output_by_tool.csv         # 每种方法的详细结果
# └── Host_crispr_links.csv              # CRISPR匹配证据

# iPHoP整合的12种方法包括：
# 1. CRISPR匹配（最可靠）
# 2. WIsH（k-mer似然）
# 3. VirHostMatcher-Net（深度学习）
# 4. PHP（Phage-Host Prediction）
# 5. RaFAH（随机森林）
# 6. 序列同源性比对
# ... 等共12种方法的集成预测

# CHERRY安装（基于图神经网络）
git clone https://github.com/KennthShang/CHERRY  # 克隆仓库
cd CHERRY  # 进入目录
pip install -r requirements.txt  # 安装依赖

# 运行CHERRY预测
python run_Cherry.py \
  --contigs viral_contigs.fasta \   # 输入噬菌体contigs
  --len 3000 \                      # 最小contig长度
  --model pretrain \                # 使用预训练模型
  --topk 1                          # 返回top-k预测

# CHERRY工作原理：
# 1. 构建蛋白质相互作用图（噬菌体蛋白 vs 细菌蛋白）
# 2. 用图神经网络(GNN)学习节点嵌入
# 3. 基于嵌入相似度预测宿主

# CHERRY输出解读
# cherry_prediction.csv包含：
# - contig_id: 噬菌体contig名
# - predicted_host: 预测的宿主（属级/种级）
# - score: 预测置信度分数

# vHULK安装（基于深度学习的噬菌体宿主预测）
conda install -c bioconda vhulk  # conda安装
# 或
pip install vhulk  # pip安装

# 运行vHULK
vhulk predict \
  --input viral_contigs.fasta \    # 输入噬菌体序列
  --output vhulk_results/ \        # 输出目录
  --threads 8                      # 线程数

# vHULK特点：
# - 基于蛋白质家族（Pfam）的特征提取
# - 使用深度神经网络(DNN)分类
# - 支持属级宿主预测
# - 不需要参考基因组数据库（模型内置）

# 基于CRISPR spacer的噬菌体-宿主匹配
import subprocess  # 命令行调用
import pandas as pd  # 数据处理
from Bio import SeqIO  # 序列处理

# 1. 从细菌基因组提取CRISPR spacer
# 使用CCTyper提取（参考671_CRISPR_Cas系统分析.md）
def extract_spacers(genome_dir, output_dir):
    """批量提取CRISPR spacer序列"""
    import os
    import glob
    all_spacers = []  # 存储所有spacer
    for genome in glob.glob(f"{genome_dir}/*.fasta"):
        sample = os.path.basename(genome).replace(".fasta", "")  # 样本名
        # 运行CCTyper
        cmd = f"cctyper {genome} {output_dir}/{sample}/ --prodigal meta"
        subprocess.run(cmd, shell=True)  # 执行命令
        # 读取spacer序列
        spacer_file = f"{output_dir}/{sample}/spacers/all_spacers.fa"
        if os.path.exists(spacer_file):  # 检查文件是否存在
            for record in SeqIO.parse(spacer_file, "fasta"):
                all_spacers.append({
                    "host_genome": sample,  # 宿主基因组
                    "spacer_id": record.id,  # spacer ID
                    "sequence": str(record.seq),  # spacer序列
                    "length": len(record.seq)  # spacer长度
                })
    return pd.DataFrame(all_spacers)

# 2. 将spacer比对到噬菌体基因组
def crispr_host_prediction(spacer_fasta, phage_fasta, output):
    """CRISPR spacer比对预测噬菌体-宿主关系"""
    # 构建噬菌体BLAST数据库
    subprocess.run([
        "makeblastdb",
        "-in", phage_fasta,  # 噬菌体基因组
        "-dbtype", "nucl",   # 核酸类型
        "-out", "phage_db"   # 数据库名
    ])
    # BLASTn比对spacer到噬菌体
    subprocess.run([
        "blastn",
        "-query", spacer_fasta,  # spacer序列
        "-db", "phage_db",       # 噬菌体数据库
        "-out", output,          # 输出文件
        "-outfmt", "6 qseqid sseqid pident length mismatch qlen slen evalue",
        "-evalue", "1e-5",       # E值阈值
        "-word_size", "7",       # 短序列用小word_size
        "-max_target_seqs", "10",  # 每个spacer最多10个hit
        "-num_threads", "8"      # 线程数
    ])

# 3. 解析匹配结果
def parse_crispr_matches(blast_output, min_identity=90, min_coverage=0.8):
    """解析CRISPR-噬菌体匹配结果"""
    columns = ["spacer_id", "phage_id", "identity", "length",
               "mismatch", "spacer_len", "phage_len", "evalue"]
    results = pd.read_csv(blast_output, sep='\t', names=columns)

    # 过滤高质量匹配
    good_hits = results[
        (results["identity"] >= min_identity) &  # 一致性≥90%
        (results["length"] / results["spacer_len"] >= min_coverage)  # 覆盖≥80%
    ].copy()

    # 提取宿主信息（从spacer_id中解析）
    good_hits["host"] = good_hits["spacer_id"].apply(
        lambda x: x.rsplit("_", 2)[0]  # 提取宿主名
    )

    # 统计噬菌体-宿主关系
    interactions = good_hits.groupby(["phage_id", "host"]).agg(
        n_spacers=("spacer_id", "nunique"),  # 匹配的spacer数
        avg_identity=("identity", "mean"),    # 平均一致性
        min_evalue=("evalue", "min")          # 最小E值
    ).reset_index()

    print(f"鉴定到 {len(interactions)} 对噬菌体-宿主关系")
    print(f"涉及 {interactions['phage_id'].nunique()} 个噬菌体")
    print(f"涉及 {interactions['host'].nunique()} 个宿主")
    return interactions

# 4. 可视化噬菌体-宿主网络
def plot_phage_host_network(interactions, output_file="phage_host_network.png"):
    """绘制噬菌体-宿主互作网络"""
    import matplotlib.pyplot as plt  # 绑图
    import networkx as nx  # 网络图

    G = nx.Graph()  # 创建无向图
    for _, row in interactions.iterrows():
        G.add_edge(
            row["phage_id"], row["host"],  # 添加边
            weight=row["n_spacers"]  # 权重=spacer数
        )

    # 设置节点颜色
    colors = []
    for node in G.nodes():
        if node in interactions["phage_id"].values:
            colors.append("#FF6B6B")  # 噬菌体=红色
        else:
            colors.append("#4ECDC4")  # 宿主=青色

    fig, ax = plt.subplots(figsize=(14, 10))
    pos = nx.spring_layout(G, seed=42)  # 弹簧布局
    nx.draw(G, pos, node_color=colors, with_labels=True,
            node_size=800, font_size=8, ax=ax)  # 绘制网络
    ax.set_title("噬菌体-宿主互作网络（基于CRISPR证据）")
    plt.tight_layout()
    plt.savefig(output_file, dpi=150)
    print(f"网络图已保存: {output_file}")

# 整合多种预测方法的结果
import pandas as pd  # 数据处理

def integrate_predictions(iphop_file, cherry_file, crispr_file):
    """整合iPHoP、CHERRY和CRISPR三种方法的预测结果"""
    # 读取各方法结果
    iphop = pd.read_csv(iphop_file)  # iPHoP结果
    cherry = pd.read_csv(cherry_file)  # CHERRY结果
    crispr = pd.read_csv(crispr_file)  # CRISPR结果

    # 标准化列名
    iphop = iphop.rename(columns={"Virus": "phage_id",
                                   "Host genus": "host_genus"})
    cherry = cherry.rename(columns={"contig_id": "phage_id",
                                     "predicted_host": "host_genus"})
    crispr = crispr.rename(columns={"phage_id": "phage_id",
                                     "host": "host_genus"})

    # 合并结果
    all_phages = set(iphop["phage_id"]) | set(cherry["phage_id"]) | set(crispr["phage_id"])
    consensus = []  # 共识结果
    for phage in all_phages:
        predictions = {}  # 各方法预测
        # iPHoP预测
        ip = iphop[iphop["phage_id"] == phage]
        if len(ip) > 0:
            predictions["iPHoP"] = ip.iloc[0]["host_genus"]
        # CHERRY预测
        ch = cherry[cherry["phage_id"] == phage]
        if len(ch) > 0:
            predictions["CHERRY"] = ch.iloc[0]["host_genus"]
        # CRISPR预测
        cr = crispr[crispr["phage_id"] == phage]
        if len(cr) > 0:
            predictions["CRISPR"] = cr.iloc[0]["host_genus"]

        # 计算共识
        if len(predictions) >= 2:  # 至少2种方法预测
            from collections import Counter
            host_counts = Counter(predictions.values())  # 统计最多的宿主
            best_host = host_counts.most_common(1)[0]  # 最常见预测
            consensus.append({
                "phage_id": phage,
                "consensus_host": best_host[0],  # 共识宿主
                "agreement": best_host[1],  # 一致方法数
                "n_methods": len(predictions),  # 总方法数
                "confidence": "high" if best_host[1] == len(predictions)
                              else "medium"  # 置信度
            })

    consensus_df = pd.DataFrame(consensus)
    print(f"共识预测: {len(consensus_df)} 个噬菌体")
    print(f"高置信度: {(consensus_df['confidence']=='high').sum()}")
    return consensus_df

# 噬菌体-宿主预测流程
噬菌体contigs → 多方法预测
  方法1: iPHoP（整合12种方法，推荐首选）
  方法2: CHERRY（图神经网络）
  方法3: CRISPR spacer匹配（最可靠证据）
  方法4: vHULK/RaFAH（深度学习/随机森林）
  → 共识预测（≥2种方法一致=高置信度）
  → 噬菌体-宿主网络可视化

# 工具选择
综合预测: iPHoP（最推荐，准确率最高）
图神经网络: CHERRY（快速，支持新型噬菌体）
直接证据: CRISPR匹配（最可靠但覆盖率低）
快速筛查: RaFAH（随机森林，速度快）

# 证据可靠性排序
CRISPR匹配 > 序列同源性 > 受体蛋白匹配 > k-mer/ML预测

# iPHoP置信度阈值
≥90%: 高置信度（推荐使用）
75-90%: 中置信度（需其他证据支持）
<75%: 低置信度（仅供参考）