从NCBI下载到生成进化树：用Prokka+Roary完成细菌泛基因组分析的完整实战记录（附批量脚本）
从NCBI下载到生成进化树：用Prokka+Roary完成细菌泛基因组分析的完整实战记录（附批量脚本）
在微生物基因组学研究中，泛基因组分析已成为揭示物种内基因多样性的重要手段。本文将手把手带您完成从NCBI数据下载到最终生成系统发育树的全流程，特别适合需要在短时间内掌握完整分析链条的研究人员。不同于碎片化的教程，我们聚焦于实际科研场景中的连贯操作，并提供可直接复用的批量处理脚本。
1. 数据获取与预处理
1.1 NCBI基因组数据批量下载
首先登录NCBI Assembly数据库，在搜索框输入目标菌种名称（如"Bacillus subtilis"）。建议设置以下筛选条件：
- Assembly level: Complete Genome或Chromosome
- RefSeq category: Reference/Representative genome
- 排除contig/scaffold级别的组装
勾选所需基因组后，使用"Download Assemblies"功能，选择：
- File type: GenBank (.fna)
- Include: Genomic sequence (FASTA)
实用技巧：对于大规模下载，推荐使用NCBI的datasets命令行工具：
datasets download genome accession --inputfile accessions.txt --filename ncbi_dataset.zip其中accessions.txt包含需要下载的基因组Accession编号列表。
1.2 数据质量控制
解压下载的压缩包后，建议进行基本质量检查：
- 使用
seqkit stats快速统计序列信息：
seqkit stats *.fna > genome_stats.txt重点关注：
- 序列数量（应为1条/基因组）
- 平均序列长度（应与物种预期基因组大小匹配）
典型问题处理：若发现多个contig的基因组，建议：
- 重新筛选更完整的组装版本
- 或使用
--compliant参数运行Prokka
2. 批量基因组注释流程
2.1 Prokka环境配置
推荐使用conda管理分析环境：
conda create -n prokka -c bioconda prokka=1.14.6 conda activate prokka prokka --setupdb # 关键步骤：初始化数据库2.2 自动化注释脚本
创建run_prokka.sh脚本实现批量处理：
#!/bin/bash mkdir -p gff_output for fna in *.fna; do sample=$(basename $fna .fna) prokka $fna \ --prefix $sample \ --outdir ${sample}_annot \ --cpus 8 \ --force # 覆盖已有结果 cp ${sample}_annot/*.gff gff_output/ done关键参数说明：
--cpus：根据服务器资源调整--force：确保重复运行时覆盖旧结果--addgenes：添加gene特征（Roary分析必需）
3. 泛基因组分析核心步骤
3.1 Roary安装与配置
虽然可通过conda安装，但推荐系统级安装确保稳定性：
sudo apt-get install -y roary mafft fasttree3.2 运行泛基因组分析
准备roary_analysis.sh脚本：
#!/bin/bash roary -f roary_results \ -p 16 \ -e \ # 使用MAFFT进行多序列比对 -v \ # 详细日志 -i 90 \ # 蛋白相似度阈值 gff_output/*.gff输出文件解读：
gene_presence_absence.csv：核心/可变基因矩阵summary_statistics.txt：泛基因组大小等统计量core_gene_alignment.aln：核心基因多序列比对
3.3 核心基因系统发育树构建
使用FastTree快速生成进化树：
FastTree -nt -gtr < roary_results/core_gene_alignment.aln > core_genome_tree.nwk对于更精确的分析，可考虑：
iqtree -s core_gene_alignment.aln -m MFP -bb 1000 -nt AUTO4. 结果可视化与解读
4.1 泛基因组特征分析
使用roary_plots.py生成基础图表：
roary_plots.py roary_results/core_gene_alignment.aln roary_results/gene_presence_absence.csv将生成：
- 泛基因组累积曲线（pangenome_matrix.png）
- 核心基因热图（clustered_pangenome_matrix.png）
4.2 进化树美化
推荐使用FigTree进行树形调整，或R语言ggtree包：
library(ggtree) tree <- read.tree("core_genome_tree.nwk") ggtree(tree) + geom_tiplab() + xlim(0, 0.2) + geom_nodepoint(color="red", alpha=0.5)4.3 高级分析拓展
对于深入挖掘，可尝试：
- Panaroo：改进的泛基因组分析工具
- pyseer：基因-表型关联分析
- Scoary：基因与表型特征关联研究
5. 常见问题解决方案
Q1: Prokka运行时出现数据库错误
prokka --listdb # 检查数据库状态 prokka --cleandb # 清理损坏的数据库 prokka --setupdb # 重新初始化Q2: Roary运行内存不足
- 增加
-s参数降低blastp严格度 - 使用
--split_paralogs处理旁系同源基因
Q3: 进化树分支支持率低
- 增加比对基因数量（降低
-i参数） - 使用IQ-TREE代替FastTree
实际项目中，处理30个芽孢杆菌基因组时，使用32核服务器完整流程约需6-8小时。建议将中间结果（特别是GFF文件）妥善保存，便于后续重新分析。