667 原核生物比较基因组ANI¶

# FastANI安装
conda install -c bioconda fastani  # conda安装

# 1. 两两比较
fastANI \
  -q query_genome.fna \       # 查询基因组
  -r reference_genome.fna \   # 参考基因组
  -o ani_result.txt            # 输出结果

# 2. 多基因组批量比较
# 准备基因组列表文件
ls genomes/*.fna > genome_list.txt  # 所有基因组路径

# 全 vs 全比较
fastANI \
  --ql genome_list.txt \      # 查询列表
  --rl genome_list.txt \      # 参考列表
  -o all_vs_all_ani.txt \     # 输出
  -t 16 \                     # 线程数
  --matrix                    # 输出矩阵格式

# 输出格式：query  reference  ANI  fragments_mapped  total_fragments
# 只输出ANI≥80%的结果

# skani安装
conda install -c bioconda skani  # conda安装

# 1. 两两比较
skani dist \
  query.fna reference.fna \   # 输入两个基因组
  -o result.txt                # 输出

# 2. 批量比较（sketch + dist 两步）
# 第一步：预计算sketch（只需做一次）
skani sketch \
  genomes/*.fna \              # 所有基因组
  -o sketches/ \               # sketch输出目录
  -t 16                        # 线程数

# 第二步：计算距离
skani triangle \
  genomes/*.fna \              # 输入基因组
  -o triangle_ani.txt \        # 输出三角矩阵
  -t 16                        # 线程数

# 3. 搜索最近参考基因组
skani search \
  query_mag.fna \              # 查询MAG
  -d reference_sketches/ \     # 参考sketch库
  -o search_result.txt \       # 输出
  -n 5                         # 返回前5个最近的

# GGDC (Genome-to-Genome Distance Calculator) 在线使用
# 网址: https://ggdc.dsmz.de/ggdc.php

# 或使用命令行版本
# Formula 2 (推荐): 基于HSP长度/总基因组长度
# dDDH ≥ 70% → 同一物种

# Python调用GGDC API（示例）

# Python计算ANI和可视化
import subprocess  # 调用命令行
import pandas as pd  # 数据处理
import numpy as np  # 数值计算
import seaborn as sns  # 热图
import matplotlib.pyplot as plt  # 绑图

def parse_fastani_output(filepath):
    """解析FastANI输出为矩阵"""
    data = []
    with open(filepath) as f:
        for line in f:
            parts = line.strip().split('\t')
            query = parts[0].split('/')[-1].replace('.fna', '')  # 提取文件名
            ref = parts[1].split('/')[-1].replace('.fna', '')
            ani = float(parts[2])  # ANI值
            data.append({"query": query, "ref": ref, "ANI": ani})

    df = pd.DataFrame(data)
    # 转为方阵
    matrix = df.pivot(index="query", columns="ref", values="ANI")
    matrix = matrix.fillna(100)  # 对角线=100%
    return matrix

# 解析结果
ani_matrix = parse_fastani_output("all_vs_all_ani.txt")

# 热图可视化
fig, ax = plt.subplots(figsize=(12, 10))
sns.heatmap(ani_matrix,
            annot=True, fmt=".1f",  # 显示数值
            cmap="RdYlGn",  # 颜色方案
            vmin=75, vmax=100,  # 色彩范围
            linewidths=0.5,
            ax=ax)
ax.set_title("ANI热图 (≥95%为同种)")
plt.tight_layout()
plt.savefig("ani_heatmap.png", dpi=150)

# 自动物种划分
def delineate_species(ani_matrix, threshold=95):
    """基于ANI阈值自动划分物种"""
    from scipy.cluster.hierarchy import fcluster, linkage
    from scipy.spatial.distance import squareform

    # ANI转距离
    dist = 100 - ani_matrix.values  # 距离 = 100 - ANI
    np.fill_diagonal(dist, 0)  # 对角线为0

    # 层次聚类
    linkage_matrix = linkage(squareform(dist), method='average')  # UPGMA
    clusters = fcluster(linkage_matrix, t=100-threshold, criterion='distance')

    # 输出物种组
    species_groups = {}
    for genome, cluster_id in zip(ani_matrix.index, clusters):
        species_groups.setdefault(cluster_id, []).append(genome)

    print(f"按{threshold}% ANI阈值，划分为 {len(species_groups)} 个物种:")
    for sp_id, members in species_groups.items():
        print(f"  物种{sp_id}: {', '.join(members)}")

    return species_groups

# PyOrthoANI/PyFastANI/Pyskani (2025 NAR Genomics)
# pip install pyorthoani pyfastani pyskani

# PyFastANI使用
import pyfastani  # 导入PyFastANI

# 创建ANI计算器
mapper = pyfastani.Mapper()  # 初始化
mapper.add_draft("ref_genome.fna")  # 添加参考基因组

# 计算ANI
hits = mapper.query_draft("query_genome.fna")  # 查询
for hit in hits:
    print(f"ANI: {hit.identity:.2f}%, 比对片段: {hit.matches}/{hit.fragments}")

# Pyskani使用（最快）
import pyskani

# 批量计算
db = pyskani.Database()
db.sketch("genome1.fna", name="genome1")  # sketch
db.sketch("genome2.fna", name="genome2")
results = db.dist()  # 计算距离
for r in results:
    print(f"{r.query} vs {r.reference}: ANI={r.identity:.2f}%")

# AAI用于属级以上的分类比较
# CompareM工具
comparem aai_wf \
  genome_dir/ \               # 基因组目录
  output_dir/ \               # 输出目录
  --proteins \                # 使用蛋白质序列
  --cpus 16                   # CPU数

# AAI阈值（粗略）
# >90% → 同属
# >60% → 同科
# >45% → 同目

# ANI分析流程
1. 基因组质量评估 (CheckM/CheckM2)
2. 选择ANI工具:
   - 少量高质量基因组 → FastANI
   - 大量/碎片化基因组(MAG) → skani
   - 需要金标准精度 → pyani (ANIb)
   - Python集成 → PyFastANI/Pyskani (2025)
3. 计算全 vs 全 ANI矩阵
4. 基于95%阈值划分物种
5. 热图可视化 + 系统发育树

# 物种划分阈值速查
ANI ≥ 95%: 同一物种
ANI 83-95%: 同属不同种
ANI < 83%: 不同属
dDDH ≥ 70%: 同一物种
AAI > 90%: 同属

# 工具速度排名（快→慢）
skani >> FastANI >> OrthoANI > ANIm > ANIb