179_宏基因组宿主预测iPHoP¶

# 安装iPHoP（推荐conda）
conda create -n iphop -c conda-forge -c bioconda iphop
conda activate iphop

# 下载数据库（约150GB解压后）
iphop download --db_dir iphop_db --no_prompt

# 验证安装
iphop --version
iphop predict -h

# 从geNomad输出获取病毒序列
cp genomad_output/contigs_summary/contigs_virus.fna viral_contigs.fa

# 检查序列统计
seqkit stats viral_contigs.fa

# 过滤短序列（建议≥5kb）
seqkit seq -m 5000 viral_contigs.fa > viral_5kb.fa
echo "序列数量: $(grep -c '>' viral_5kb.fa)"

# 运行iPHoP
iphop predict \
    --fa_file viral_5kb.fa \
    --db_dir iphop_db/ \
    --out_dir iphop_output/ \
    --num_threads 16 \
    --min_score 90

# 参数说明
# --fa_file: 病毒序列FASTA
# --db_dir: 数据库路径
# --out_dir: 输出目录
# --num_threads: 线程数
# --min_score: 最低置信度阈值（90=宽松，95=严格）

# 使用自定义宿主数据库（可选，加入自有MAGs）
iphop add_to_db \
    --fna_dir custom_host_genomes/ \
    --gtdb_dir gtdb_taxonomy/ \
    --out_dir custom_db/ \
    --num_threads 16

# 使用自定义数据库预测
iphop predict \
    --fa_file viral_5kb.fa \
    --db_dir custom_db/ \
    --out_dir iphop_custom_output/ \
    --num_threads 16

# 输出目录结构
ls iphop_output/
# Host_prediction_to_genus_m90.csv   # 属级别预测（主要结果）
# Host_prediction_to_genome_m90.csv  # 基因组级别预测
# Detailed_output_by_tool.csv        # 各工具的详细预测
# Summary.txt                         # 汇总统计

# 查看主要结果
head -5 iphop_output/Host_prediction_to_genus_m90.csv

import pandas as pd
import matplotlib.pyplot as plt

# 读取结果
host_pred = pd.read_csv("iphop_output/Host_prediction_to_genus_m90.csv")

# 查看列名
print(host_pred.columns.tolist())
# ['Virus', 'AAI to closest RaFAH reference', 'Host genus',
#  'Confidence score', 'List of methods', 'Host taxonomy']

# 基本统计
total_viruses = len(pd.read_csv("viral_5kb.fa.fai", sep="\t", header=None))
predicted = len(host_pred)
print(f"总病毒数: {total_viruses}")
print(f"成功预测数: {predicted}")
print(f"预测率: {predicted/total_viruses*100:.1f}%")

# 宿主分类分布
# 提取门级别
host_pred['Host_phylum'] = host_pred['Host taxonomy'].apply(
    lambda x: x.split(';')[1] if pd.notna(x) and ';' in str(x) else 'Unknown'
)

# Top 10宿主门
phylum_counts = host_pred['Host_phylum'].value_counts().head(10)
print("\nTop 10 宿主门:")
print(phylum_counts)

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 宿主门分布
phylum_counts.plot(kind='barh', ax=axes[0], color='steelblue')
axes[0].set_xlabel('Number of Viruses')
axes[0].set_title('Host Phylum Distribution')
axes[0].invert_yaxis()

# 置信度分布
axes[1].hist(host_pred['Confidence score'], bins=30, edgecolor='black', alpha=0.7)
axes[1].set_xlabel('Confidence Score')
axes[1].set_ylabel('Count')
axes[1].set_title('Prediction Confidence Distribution')

plt.tight_layout()
plt.savefig('iphop_results.png', dpi=300)

# 整合iPHoP与geNomad结果
virus_info = pd.read_csv("genomad_output/contigs_virus_summary.tsv", sep="\t")
host_info = pd.read_csv("iphop_output/Host_prediction_to_genus_m90.csv")

# 合并
merged = virus_info.merge(host_info, left_on='seq_name', right_on='Virus', how='left')

# 输出整合表
merged[['seq_name', 'length', 'virus_score', 'taxonomy',
        'Host genus', 'Confidence score']].to_csv("virus_host_integrated.tsv",
                                                     sep="\t", index=False)

# 病毒-宿主网络数据
network_data = merged[['taxonomy', 'Host genus']].dropna()
network_data.to_csv("virus_host_network.tsv", sep="\t", index=False)

# 完整流程
iphop predict --fa_file viruses.fa --db_dir iphop_db/ --out_dir output/ -t 16

# 宽松阈值（更多预测，精度略低）
iphop predict --fa_file viruses.fa --db_dir iphop_db/ --out_dir output_m75/ -t 16 --min_score 75

# 严格阈值（更少预测，精度更高）
iphop predict --fa_file viruses.fa --db_dir iphop_db/ --out_dir output_m95/ -t 16 --min_score 95

# 添加自定义宿主基因组
iphop add_to_db --fna_dir my_mags/ --gtdb_dir gtdb/ --out_dir custom_db/ -t 16

# 分割大文件并行处理
seqkit split2 -s 1000 viruses.fa -O split_dir/
for f in split_dir/*.fa; do
    iphop predict --fa_file $f --db_dir iphop_db/ --out_dir ${f%.fa}_out/ -t 4 &
done
wait
cat split_dir/*_out/Host_prediction_to_genus_m90.csv > all_predictions.csv