原核基因组比较分析¶
一句话概述¶
通过Roary/Panaroo进行泛基因组分析确定核心基因和可变基因，计算ANI值判定物种边界，利用核心基因构建高分辨率系统发育树，是细菌分类学、进化分析和功能基因组学研究的标准流程。
核心知识点表格¶
| 知识点 | 关键内容 | 常用工具 |
|---|---|---|
| 泛基因组概念 | 核心基因组+可变基因组+独有基因 | Roary, Panaroo, PPanGGOLiN |
| ANI计算 | 平均核苷酸一致性（物种划分金标准95%） | FastANI, skani, pyani, OrthoANI |
| 基因预测 | 原核ORF识别与注释 | Prokka, Bakta, PGAP |
| 直系同源聚类 | 基因家族分组 | OrthoFinder, CD-HIT |
| 核心基因系统发育 | 保守基因联合建树 | IQ-TREE, RAxML |
| 基因组共线性 | 大尺度基因排列比较 | Mauve, MUMmer, SyRI |
| 基因岛检测 | 水平转移区域识别 | IslandViewer, Alien_Hunter |
| 泛基因组可视化 | 基因present/absent热图 | Phandango, Anvi'o |
各步骤详解¶
第一步：基因组注释标准化¶
白话解释： 比较基因组分析的第一步是确保所有基因组使用相同的注释流程。如果不同基因组的注释标准不一致，后续的同源基因比较就会有偏差。通常用Prokka或Bakta统一重新注释所有基因组。
技术细节： - Prokka: 快速原核基因组注释，输出GFF3格式 - Bakta: 更新的工具，使用UniProt等数据库，注释质量更高 - 关键：所有样本使用相同版本、相同数据库注释 - 注释质量影响泛基因组分析的准确性
代码示例：
# 安装
conda install -c bioconda prokka bakta
# 方法1: Prokka批量注释
for genome in genomes/*.fasta; do
name=$(basename $genome .fasta)
prokka --outdir prokka_output/$name \
--prefix $name \
--kingdom Bacteria \
--genus Staphylococcus \
--species aureus \
--strain $name \
--cpus 8 \
--compliant \
--force \
$genome
done
# 方法2: Bakta注释（更推荐）
# 先下载数据库
bakta_db download --output /db/bakta_db
for genome in genomes/*.fasta; do
name=$(basename $genome .fasta)
bakta --db /db/bakta_db \
--output bakta_output/$name \
--prefix $name \
--genus Staphylococcus \
--species aureus \
--threads 8 \
$genome
done
# 收集所有GFF文件
mkdir gff_files
cp prokka_output/*/*.gff gff_files/
第二步：ANI计算与物种划分¶
白话解释： ANI（平均核苷酸一致性）就是两个基因组之间共享区域的平均序列相似度。如果ANI≥95%，通常认为是同一个种。这相当于用数字来精确定义"这两个细菌是不是同一个种"。
技术细节： - ANI 95-96% 等价于传统的70% DNA-DNA杂交阈值 - FastANI使用MinHash算法快速估计，适合大规模比较 - pyani使用BLAST或MUMmer精确计算 - 还有AAI（氨基酸一致性）用于更远的进化关系 - dDDH（数字DNA-DNA杂交）是另一种物种划分指标
代码示例：
# 方法1: FastANI (推荐，快速)
conda install -c bioconda fastani
# 准备文件列表
ls genomes/*.fasta > genome_list.txt
# 全部两两比较
fastANI --ql genome_list.txt \
--rl genome_list.txt \
-o fastani_output.txt \
-t 16 \
--matrix
# 输出格式: query reference ANI orthologous_matches total_fragments
# 示例: genome_A.fasta genome_B.fasta 98.7654 800 850
# 方法2: pyani (更精确)
pip install pyani
# 使用ANIb (BLAST)或ANIm (MUMmer)
average_nucleotide_identity.py \
-i genomes/ \
-o pyani_output/ \
-m ANIb \
--workers 8 \
-g --gformat pdf
# 方法3: skani (最快，比FastANI快>20倍，对MAGs更鲁棒)
conda install -c bioconda skani
skani dist -q genomes/*.fasta -r genomes/*.fasta \
-o skani_results.txt -t 16
# skani使用稀疏链式比对（sparse chaining），适合大规模和低质量MAGs
# 绘制ANI热图
python3 << 'EOF'
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
# 读取FastANI结果
df = pd.read_csv("fastani_output.txt", sep="\t", header=None,
names=["query","ref","ANI","ortho","total"])
# 构建矩阵
genomes = sorted(set(df["query"].tolist()))
matrix = pd.DataFrame(100.0, index=genomes, columns=genomes)
for _, row in df.iterrows():
q = row["query"]
r = row["ref"]
matrix.loc[q, r] = row["ANI"]
# 绘制热图
plt.figure(figsize=(12,10))
sns.clustermap(matrix, cmap="YlOrRd", vmin=90, vmax=100,
annot=True, fmt=".1f", figsize=(12,10))
plt.savefig("ANI_heatmap.pdf")
EOF
第三步：泛基因组分析（Roary/Panaroo）¶
白话解释： 泛基因组就是一个物种所有菌株的基因总和。其中"核心基因组"是所有菌株都有的基因（维持基本生命活动），"可变基因组"是部分菌株有的基因（可能与特殊环境适应相关），"独有基因"只在个别菌株中出现。
技术细节： - Roary: 经典工具，使用CD-HIT+MCL聚类 - Panaroo: 更新工具，使用图算法校正注释错误，假阳性更少 - 核心基因定义阈值：通常99%（严格）或95%（软核心） - 泛基因组开放性：Heap's law拟合判断泛基因组是否开放
代码示例：
# === 方法1: Roary ===
conda install -c bioconda roary
# 运行Roary
roary -p 16 \
-f roary_output \
-e -n \ # 生成核心基因比对
-cd 99 \ # 核心基因阈值99%
-i 95 \ # 序列一致性阈值95%
-s \ # 不拆分同源基因
-v \ # 详细输出
gff_files/*.gff
# 输出文件：
# gene_presence_absence.csv - 基因存在/缺失矩阵
# core_gene_alignment.aln - 核心基因比对
# summary_statistics.txt - 统计汇总
# pan_genome_reference.fa - 泛基因组参考序列
# 查看统计
cat roary_output/summary_statistics.txt
# Core genes (99% <= strains <= 100%): 2,145
# Soft core genes (95% <= strains < 99%): 156
# Shell genes (15% <= strains < 95%): 1,823
# Cloud genes (0% < strains < 15%): 4,567
# Total genes: 8,691
# === 方法2: Panaroo (推荐) ===
conda install -c bioconda panaroo
# 运行Panaroo
panaroo -i gff_files/*.gff \
-o panaroo_output \
--clean-mode strict \ # 严格清理错误注释
-t 16 \
--core_threshold 0.99 \
--remove-invalid-genes \
-a core \ # 输出核心基因比对
--aligner mafft
# Panaroo优势：
# 1. 使用图算法检测并修正注释错误
# 2. 处理基因分裂/融合事件
# 3. 三种清理模式：strict/moderate/sensitive
# 查看结果
head panaroo_output/gene_presence_absence_roary.csv
# 与Roary输出格式兼容
第四步：核心基因系统发育¶
白话解释： 用所有菌株都共有的核心基因来建系统发育树，这比用单个基因（如16S rRNA）分辨率高得多。核心基因联合分析提供了基因组水平的进化关系证据。
技术细节： - 核心基因比对连接（concatenation）后建树 - 或者用基因树+物种树方法（coalescent-based） - SNP-based方法更快但信息较少 - Recombination可能影响树的准确性，需要检测和去除
代码示例：
# 从Roary/Panaroo核心基因比对建树
# 方法1: 直接使用核心基因alignment
iqtree2 -s panaroo_output/core_gene_alignment.aln \
-m GTR+G4 \
-B 1000 \
--alrt 1000 \
-T AUTO \
--prefix core_phylo
# 方法2: SNP-based快速建树
# 提取SNP位点
snp-sites -c panaroo_output/core_gene_alignment.aln \
> core_snps.fasta
# 用SNP建树（更快）
iqtree2 -s core_snps.fasta \
-m GTR+ASC+G \ # ASC校正：只有变异位点
-B 1000 \
-T AUTO \
--prefix snp_phylo
# 方法3: 去除重组区域后建树
# 使用Gubbins检测重组
conda install -c bioconda gubbins
run_gubbins.py \
panaroo_output/core_gene_alignment.aln \
--threads 16 \
--prefix gubbins_output \
--tree-builder iqtree
# Gubbins输出去除重组后的alignment
iqtree2 -s gubbins_output.filtered_polymorphic_sites.fasta \
-m GTR+G -B 1000 -T AUTO --prefix recomb_free_phylo
# 方法4: ClonalFrameML（ML框架下的重组检测）
ClonalFrameML core_phylo.treefile \
panaroo_output/core_gene_alignment.aln \
clonalframe_output
第五步：泛基因组可视化与功能分析¶
白话解释： 将泛基因组的基因存在/缺失矩阵与系统发育树结合展示，可以直观看到哪些基因在哪些菌株中存在。同时对可变基因进行功能富集分析，理解细菌环境适应的遗传基础。
代码示例：
# 使用Phandango在线可视化
# 上传: 树文件 + gene_presence_absence.csv
# 网址: https://jameshadfield.github.io/phandango/
# 使用roary_plots.py
python roary_plots.py core_phylo.treefile \
roary_output/gene_presence_absence.csv
# 使用R进行自定义可视化
Rscript << 'EOF'
library(ggtree)
library(pheatmap)
library(tidyverse)
# 读取树
tree <- read.tree("core_phylo.treefile")
# 读取基因存在/缺失矩阵
pa <- read.csv("roary_output/gene_presence_absence.csv",
check.names=FALSE)
# 转为二值矩阵
binary_matrix <- pa %>%
select(Gene, starts_with("genome_")) %>%
column_to_rownames("Gene") %>%
mutate(across(everything(), ~ifelse(is.na(.), 0, 1)))
# 绘制树+热图
p <- ggtree(tree) + geom_tiplab(size=2)
gheatmap(p, t(binary_matrix[1:50,]),
colnames=FALSE, width=2) +
scale_fill_manual(values=c("white","steelblue"))
ggsave("pangenome_tree_heatmap.pdf", width=15, height=10)
EOF
# 泛基因组曲线（Heap's law）
Rscript << 'EOF'
library(micropan)
# 读取基因矩阵
pa_matrix <- read.csv("gene_presence_absence_binary.csv", row.names=1)
# 随机抽样绘制泛基因组增长曲线
heap <- heaps(pa_matrix, n.perm=100)
# alpha < 1: 开放泛基因组（新基因组总能带来新基因）
# alpha > 1: 封闭泛基因组
cat("Heaps law alpha:", heap$alpha, "\n")
EOF
第六步：基因组共线性与结构变异¶
白话解释： 除了基因的有无，基因在染色体上的排列顺序也会发生变化（倒位、转位、重复）。基因组共线性分析可以发现这些大尺度结构变异。
代码示例：
# 使用MUMmer进行两两比较
conda install -c bioconda mummer
nucmer --prefix=ref_vs_query \
reference.fasta query.fasta
# 过滤比对结果
delta-filter -1 ref_vs_query.delta > ref_vs_query.filter
# 生成坐标文件
show-coords -rcl ref_vs_query.filter > ref_vs_query.coords
# 绘制dot plot
mummerplot --png --large \
-p dotplot \
ref_vs_query.filter
# 使用Mauve进行多基因组比对
progressiveMauve --output=mauve_output.xmfa \
genome1.fasta genome2.fasta genome3.fasta
# 使用pyGenomeViz可视化
pip install pygenomeviz
python3 << 'EOF'
from pygenomeviz import GenomeViz
gv = GenomeViz()
gv.set_scale_bar()
# 添加基因组track
track1 = gv.add_feature_track("Genome_A", 3000000)
track2 = gv.add_feature_track("Genome_B", 2950000)
# 添加link（共线性区域）
# 从MUMmer coords文件解析
gv.savefig("synteny_plot.pdf")
EOF
实战命令¶
#!/bin/bash
# === 原核基因组比较分析完整流程 ===
set -e
GENOMES_DIR="genomes"
THREADS=16
OUTPUT="comparative_analysis"
mkdir -p $OUTPUT
# 1. 统一注释
mkdir -p $OUTPUT/gff
for fasta in $GENOMES_DIR/*.fasta; do
name=$(basename $fasta .fasta)
prokka --outdir $OUTPUT/prokka/$name --prefix $name \
--cpus $THREADS --kingdom Bacteria $fasta
cp $OUTPUT/prokka/$name/$name.gff $OUTPUT/gff/
done
# 2. ANI计算
ls $GENOMES_DIR/*.fasta > $OUTPUT/genome_list.txt
fastANI --ql $OUTPUT/genome_list.txt --rl $OUTPUT/genome_list.txt \
-o $OUTPUT/fastani_results.txt -t $THREADS --matrix
# 3. 泛基因组分析
panaroo -i $OUTPUT/gff/*.gff -o $OUTPUT/panaroo \
--clean-mode strict -t $THREADS -a core --aligner mafft
# 4. 去重组+建树
run_gubbins.py $OUTPUT/panaroo/core_gene_alignment.aln \
--threads $THREADS --prefix $OUTPUT/gubbins
iqtree2 -s $OUTPUT/gubbins.filtered_polymorphic_sites.fasta \
-m GTR+G -B 1000 -T AUTO --prefix $OUTPUT/final_tree
echo "Analysis complete!"
面试常问点¶
Q1: 什么是泛基因组？核心基因和可变基因的生物学意义？ A: 泛基因组是一个物种所有基因的总和。核心基因（所有菌株共有）编码基本生命功能（转录、翻译、代谢）；可变基因与环境适应相关（抗生素耐药、毒力因子、代谢特殊底物）；独有基因可能来自水平转移或噬菌体。
Q2: ANI 95%阈值的理论基础是什么？ A: 95% ANI等价于传统的70% DNA-DNA杂交值，是公认的原核物种划分标准。生物学上对应约10^8年的分化时间。但对某些类群（如Bacillus cereus group）这个阈值可能需要调整。
Q3: Roary和Panaroo有什么区别？为什么推荐Panaroo？ A: Roary使用CD-HIT+MCL简单聚类；Panaroo使用图算法可以检测并修正注释错误（假基因、基因分裂/融合）。Panaroo的核心基因组估计更保守准确，假阳性更少，但运行更慢。
Q4: 为什么要去除重组区域再建树？ A: 细菌频繁发生同源重组，重组片段的进化历史与物种树不同。如果包含重组区域，树的拓扑会被扭曲（尤其是支长）。Gubbins/ClonalFrameML可以识别并mask这些区域。
Q5: 如何判断泛基因组是开放还是封闭的？ A: 用Heap's law拟合: n = κ·N^γ。γ>0（α<1）为开放泛基因组（每增加一个基因组都能发现新基因），常见于生态位多样的物种（如E. coli）。γ≈0（α>1）为封闭泛基因组，常见于专性寄生菌。
易错点¶
注释工具不一致：不同工具注释的基因数差异显著。必须用同一工具、同一版本重新注释所有基因组。
draft基因组的contig断裂：基因可能跨越contig边界被截断，导致泛基因组统计偏差。Panaroo对此有一定容错。
ANI计算方向性：ANI(A→B)可能≠ANI(B→A)，因为片段化程度不同。应取双向平均值。
同源基因聚类阈值：过松会合并旁系同源（paralog），过严会分裂真正的直系同源。需要根据物种分化程度调整。
忽略质粒序列：质粒上的基因可能被当作独有基因，实际是可移动元件。需要注意chromosome vs plasmid的分开处理。
样本量不足：泛基因组分析至少需要15-20个基因组才能得到稳定的统计估计。
补充知识¶
泛基因组分析工具选择指南¶
| 工具 | 适用场景 | 优势 | 局限 |
|---|---|---|---|
| Roary | 近缘菌株（同种） | 快速，经典 | 对注释错误敏感 |
| Panaroo | 近缘菌株 | 图算法修正错误 | 较慢 |
| PPanGGOLiN | 大规模分析 | 可处理数千基因组 | 不输出比对 |
| OrthoFinder | 跨属比较 | 处理基因重复 | 不限于原核 |
| GET_HOMOLOGUES | 灵活参数 | 多种聚类算法 | 学习曲线陡 |