739. 宏基因组噬菌体-宿主互作预测¶

信号1: CRISPR间隔区匹配
  细菌的CRISPR系统会记录曾感染它的噬菌体序列
  → 匹配细菌的CRISPR间隔区和噬菌体序列
  → 最可靠但覆盖率低

信号2: 序列同源性
  噬菌体在溶源时会整合到宿主基因组
  → 搜索噬菌体和细菌基因组的同源序列
  → 可能反映历史感染

信号3: 寡核苷酸频率(k-mer)
  噬菌体会适应宿主的密码子偏好
  → 比较噬菌体和细菌的k-mer频率
  → 预测范围广但精度一般

信号4: 共现模式
  噬菌体和宿主在生态上关联
  → 分析丰度的共变关系
  → 需要多样本数据

信号5: 蛋白结构匹配（2025新方法）
  噬菌体尾部蛋白决定宿主范围
  → 用蛋白折叠预测结合特异性
  → 最前沿但计算量大

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

# 下载数据库
virsorter setup -d vs2_db -j 4  # 下载VirSorter2数据库（约20GB）

# ===== 运行VirSorter2 =====
virsorter run \
  -w vs2_output/ \  # 输出目录
  -i contigs.fa \  # 输入：宏基因组组装的contigs
  -d vs2_db/ \  # 数据库路径
  --include-groups dsDNAphage,ssDNA \  # 检测双链DNA噬菌体和单链DNA病毒
  --min-length 5000 \  # 最短序列长度（5kb）
  --min-score 0.5 \  # 最低置信度
  -j 8  # 线程数

# 结果文件
# vs2_output/
# ├── final-viral-combined.fa        # 预测的病毒序列
# ├── final-viral-score.tsv          # 置信度评分
# └── final-viral-boundary.tsv       # 病毒-宿主边界

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

# 运行
VIBRANT_run.py \
  -i contigs.fa \  # 输入：contigs
  -folder vibrant_output/ \  # 输出目录
  -t 8 \  # 线程数
  -l 5000  # 最短序列长度

# 取VirSorter2和VIBRANT的交集作为高置信噬菌体

# ===== iPHoP：整合多种信号的宿主预测工具 =====
# iPHoP = integrated Phage-HOst Prediction

# 安装
conda create -n iphop python=3.8  # 创建环境
conda activate iphop
pip install iphop  # 安装iPHoP

# 下载数据库
iphop download --db_dir iphop_db/  # 下载iPHoP数据库

# 运行宿主预测
iphop predict \
  --fa_file viral_contigs.fa \  # 输入：噬菌体序列
  --db_dir iphop_db/ \  # 数据库目录
  --out_dir iphop_output/ \  # 输出目录
  --num_threads 8  # 线程数

# 结果文件
# iphop_output/
# ├── Host_prediction_to_genus_m90.csv   # 属水平宿主预测
# ├── Host_prediction_to_genome_m90.csv  # 基因组水平宿主预测
# └── Detailed_output_by_tool.csv        # 各方法的详细结果

# ===== 用CRISPR间隔区匹配宿主 =====

# 第一步：从细菌MAG中提取CRISPR间隔区
# 安装CRISPRCasFinder
conda install -c bioconda crisprcasfinder  # 安装

# 运行CRISPR检测（对每个MAG）
for mag in mags/*.fa; do
    name=$(basename ${mag} .fa)  # 提取文件名
    CRISPRCasFinder.pl \
      -in ${mag} \  # 输入：MAG文件
      -out crispr_output/${name} \  # 输出目录
      -so /path/to/sel392v2.so  # 搜索模型文件
done

# 第二步：收集所有CRISPR间隔区
cat crispr_output/*/Spacers/*.fa > all_spacers.fa  # 合并所有间隔区

# 第三步：将间隔区比对到噬菌体序列
blastn \
  -query all_spacers.fa \  # 查询：CRISPR间隔区
  -subject viral_contigs.fa \  # 目标：噬菌体序列
  -out spacer_matches.txt \  # 输出：匹配结果
  -outfmt "6 qseqid sseqid pident length mismatch" \  # 表格输出格式
  -evalue 1e-5 \  # E值阈值
  -max_target_seqs 5  # 每个间隔区最多5个匹配

# ===== 分析噬菌体-宿主预测结果 =====
import pandas as pd  # 导入pandas
import matplotlib.pyplot as plt  # 导入绑图
import seaborn as sns  # 导入seaborn

# 读取iPHoP结果
predictions = pd.read_csv("iphop_output/Host_prediction_to_genus_m90.csv")
print(f"预测了宿主的噬菌体数: {len(predictions)}")

# 统计宿主分类分布
host_genus = predictions["Host genus"].value_counts().head(20)

plt.figure(figsize=(12, 6))
host_genus.plot(kind="barh", color="#3498db")
plt.xlabel("Number of Phages")
plt.ylabel("Predicted Host Genus")
plt.title("Top 20 Predicted Host Genera")
plt.gca().invert_yaxis()
plt.tight_layout()
plt.savefig("phage_host_distribution.png", dpi=300)
plt.show()

# 读取CRISPR匹配结果
crispr_matches = pd.read_csv("spacer_matches.txt", sep="\t",
    names=["Spacer_ID", "Phage_ID", "Identity", "Length", "Mismatch"])

# 过滤高质量匹配（≥95%一致性，长度≥25bp）
high_quality = crispr_matches[
    (crispr_matches["Identity"] >= 95) &  # 一致性≥95%
    (crispr_matches["Length"] >= 25)  # 匹配长度≥25bp
]

# 从间隔区ID中提取宿主MAG名
high_quality["Host_MAG"] = high_quality["Spacer_ID"].str.extract(r"(bin\.\d+)")

# 构建噬菌体-宿主对
phage_host_pairs = high_quality[["Phage_ID", "Host_MAG"]].drop_duplicates()
print(f"\nCRISPR确认的噬菌体-宿主对: {len(phage_host_pairs)}")
print(phage_host_pairs.head(10))

# ===== 绘制噬菌体-宿主网络 =====
import networkx as nx

# 创建二部图
B = nx.Graph()
for _, row in phage_host_pairs.iterrows():
    B.add_node(row["Phage_ID"], bipartite=0, type="phage")  # 噬菌体节点
    B.add_node(row["Host_MAG"], bipartite=1, type="host")   # 宿主节点
    B.add_edge(row["Phage_ID"], row["Host_MAG"])  # 连边

# 布局
pos = nx.bipartite_layout(B, 
    [n for n, d in B.nodes(data=True) if d.get("type") == "phage"])

plt.figure(figsize=(14, 10))
phage_nodes = [n for n, d in B.nodes(data=True) if d.get("type") == "phage"]
host_nodes = [n for n, d in B.nodes(data=True) if d.get("type") == "host"]

nx.draw_networkx_nodes(B, pos, phage_nodes, node_color="#e74c3c",
                       node_size=200, label="Phage")
nx.draw_networkx_nodes(B, pos, host_nodes, node_color="#3498db",
                       node_size=300, label="Host")
nx.draw_networkx_edges(B, pos, alpha=0.3)
plt.legend(scatterpoints=1, fontsize=12)
plt.title("Phage-Host Interaction Network")
plt.axis("off")
plt.tight_layout()
plt.savefig("phage_host_network.png", dpi=300)
plt.show()

# ===== 噬菌体-宿主预测速查 =====

# 1. 鉴定噬菌体
virsorter run -i contigs.fa -d db/ -w output/ --min-length 5000

# 2. 整合宿主预测
iphop predict --fa_file viral.fa --db_dir db/ --out_dir output/

# 3. CRISPR验证
CRISPRCasFinder.pl -in mag.fa -out output/
blastn -query spacers.fa -subject viral.fa -evalue 1e-5

# 关键判断标准：
# CRISPR匹配: identity≥95%, length≥25bp
# iPHoP: confidence≥90%
# k-mer: 同一聚类的噬菌体和细菌

# 噬菌体检测工具选择：
# VirSorter2 → 最全面
# VIBRANT → 功能注释好
# geNomad → 区分质粒和噬菌体
# DeepVirFinder → 深度学习方法