从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 fasttree

3.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 AUTO

4. 结果可视化与解读

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文件）妥善保存，便于后续重新分析。