泛基因组 | 分享一套“数据下载、质控、组装、矫正、注释到泛基因组统计与绘图“的泛基因组分析组装代码

​
生信服务器活动
网址：

https://dayu.xiyoucloud.net/dayu/api/v1/anonymous/affiliate/BioinfoDu

写在前面

我们前段时间，分享了Linux中部署及使用Codex教程,若是在本地部署，也非常容易。我们直接使用Codex的API密匙即可完成部署。

对于Codex的理解和编程能力，可能未使用的人对此还是有比较大的疑惑。me too。

本期，分享一套基于Codex协助完成的泛基因组分析流程代码。

注：本套代码仅供参考学习。

若要获得本期教程数据和代码，在后台回复：20260624

结果如下

1. 项目结构

pangenome_tutorial/ ├── README.md ├── docs/ │ └── pangenome_workflow.md ├── env/ │ └── pangenome_env.yml ├── script/ │ ├── 00_create_env.sh │ ├── 01_download_public_data.sh │ ├── 02_genome_qc_from_public_assemblies.sh │ ├── 02_qc_trim_assemble.sh │ ├── 03_polish_and_correct.sh │ ├── 04_annotate.sh │ ├── 05_run_pangenome.sh │ ├── run_r_analysis.sh │ └── run_r_analysis_from_existing_matrix.sh ├── R/ │ ├── 01_generate_demo_data.R │ ├── 02_pangenome_analysis.R │ ├── 03_plot_figures.R │ └── 04_import_roary_panaroo.R ├── data/ │ ├── metadata/ │ └── pangenome/ └── results/ ├── figures/ └── tables/

2. 研究对象与公开数据来源

1.1 参考文献

Tettelin 等 2005 年发表的S. agalactiae研究对多个致病菌株进行比较基因组分析，指出单个参考基因组不能代表物种内全部遗传多样性，并提出 pan-genome 的思想。文章摘要中报告，该物种的核心基因组约占单个基因组的 80%，其余为不同菌株间差异明显的附属基因组。

本教程的数据表位于：

data/metadata/tettelin_2005_accessions.tsv data/metadata/publication.tsv

tettelin_2005_accessions.tsv中列出了用于示例重分析的 8 个基因组入口：

1.2 两种可运行入口

由于该经典文章公开沉积的是组装后的 WGS/染色体序列，而不是现代 Illumina FASTQ 原始读段，本教程将流程分为两条入口。

入口 A：复现文献公开基因组数据

cdpangenome_tutorialbashscript/01_download_public_data.shbashscript/02_genome_qc_from_public_assemblies.sh

入口 B：新项目或 SRA FASTQ 原始读段

cdpangenome_tutorialcpdata/metadata/samples_template.tsv data/metadata/samples.tsv# 编辑 data/metadata/samples.tsv，使 fq1/fq2 指向真实 fastq.gz 文件bashscript/02_qc_trim_assemble.shbashscript/03_polish_and_correct.shbashscript/04_annotate.shbashscript/05_run_pangenome.sh

2. 软件环境构建

所有软件安装定义在env/pangenome_env.yml中。推荐使用 mamba，因为生物信息依赖较多。

cdpangenome_tutorialbashscript/00_create_env.sh conda activate pangenome_env

该环境包含：

Bakta 和 CheckM 需要额外数据库。真实项目建议把数据库放在项目目录下，例如：

mkdir-pdata/db bakta_db download--outputdata/db/baktaexportBAKTA_DB="$(pwd)/data/db/bakta/db"mkdir-pdata/db/checkm checkm data setRoot data/db/checkm

3. 数据下载

3.1 文献公开基因组下载

script/01_download_public_data.sh读取data/metadata/tettelin_2005_accessions.tsv，使用 Entrez Direct 从 NCBI nucleotide 数据库下载 FASTA 和 GFF3。

bashscript/01_download_public_data.sh

输出：

data/raw/genomes/*.fna data/raw/genomes/*.gff3 logs/downloaded_genome_stats.tsv

3.2 原始 FASTQ 项目的样本表

如果分析新测序项目或 SRA 下载后的 reads，准备如下样本表：

sample_id fq1 fq2 sample_01 data/raw/reads/sample_01_R1.fastq.gz data/raw/reads/sample_01_R2.fastq.gz sample_02 data/raw/reads/sample_02_R1.fastq.gz data/raw/reads/sample_02_R2.fastq.gz

本项目提供模板：

cpdata/metadata/samples_template.tsv data/metadata/samples.tsv

4. 质控、修剪与组装

原始 reads 质控由script/02_qc_trim_assemble.sh完成。该脚本对每个样本执行：

1. FastQC：检查原始 reads 质量。
2. fastp：去接头、低质量端修剪、过滤过短 reads。
3. FastQC：检查修剪后 reads。
4. SPAdes--careful：进行细菌基因组短读长组装。
5. QUAST：评估 contig 数量、N50、总长度、GC 等指标。
6. MultiQC：整合质控报告。

bashscript/02_qc_trim_assemble.sh

推荐初筛阈值：

5. 基因组矫正与过滤

组装后使用 reads 回贴到 contig，通过 Pilon 矫正小 indel、替换和局部错误。

bashscript/03_polish_and_correct.sh

脚本步骤：

1. BWA index：为 contig 建索引。
2. BWA MEM：将修剪后 reads 比对回组装结果。
3. SAMtools sort/index：整理 BAM。
4. Pilon：根据 read alignment 矫正组装。
5. seqkit：过滤小于 500 bp 的短 contig。

矫正后的文件位于：

results/polished/<sample>/<sample>.polished.min500.fna

质量控制要点：

• 同一项目内所有样本必须采用同一组装、过滤和矫正策略。
• 不建议把完成图基因组、草图基因组、不同测序平台结果不加说明地直接混合。
• 若泛基因组中 singleton 异常偏多，应优先检查污染、碎片化和注释一致性，而不是直接解释为生物学差异。

6. 基因组注释

泛基因组分析依赖基因预测和功能注释。不同注释器会改变 CDS 边界和基因数量，因此同一个项目必须统一注释流程。

bashscript/04_annotate.sh

默认执行 Prokka，若设置BAKTA_DB则同时执行 Bakta。

Prokka 输出：

results/annotation/prokka/<sample>/<sample>.gff results/annotation/prokka/<sample>/<sample>.faa results/annotation/prokka/<sample>/<sample>.ffn

推荐实践：

• 泛基因组聚类使用同一注释器产生的 GFF。
• 发表时报告注释软件版本、数据库版本和关键参数。
• 物种名、菌株名、locustag 应统一且不包含空格。

7. 泛基因组构建

本项目同时给出 Panaroo 和 Roary。Roary 是经典工具，适合快速得到基因 presence/absence 矩阵；Panaroo 更强调错误基因、碎片基因和图结构清理，适合对草图组装进行稳健分析。

bashscript/05_run_pangenome.sh

核心参数：

panaroo\-i*.gff\-oresults/pangenome/panaroo\--clean-mode strict\--remove-invalid-genes\-acore roary\-e--mafft\-p16\-fresults/pangenome/roary\*.gff

常见结果文件：

8. R 语言下游分析与绘图

8.1 演示数据

本教程已生成可直接运行的演示矩阵：

data/pangenome/gene_presence_absence_demo.tsv data/pangenome/gene_family_annotation_demo.tsv data/pangenome/per_genome_gene_counts_demo.tsv

生成、分析和绘图：

bashscript/run_r_analysis.sh

该命令依次运行：

R/01_generate_demo_data.R R/02_pangenome_analysis.R R/03_plot_figures.R

8.2 导入真实 Panaroo/Roary 矩阵

真实项目跑完 Panaroo 后，导入矩阵：

Rscript R/04_import_roary_panaroo.R"$(pwd)"results/pangenome/panaroo/gene_presence_absence.Rtabbashscript/run_r_analysis_from_existing_matrix.sh

如果使用 Roary 输出：

Rscript R/04_import_roary_panaroo.R"$(pwd)"results/pangenome/roary/gene_presence_absence.Rtabbashscript/run_r_analysis_from_existing_matrix.sh

8.3 R 分析内容

R/02_pangenome_analysis.R完成以下统计：

• 计算每个 gene family 在多少个基因组中出现。
• 分类为 Core、Soft-core、Shell、Cloud。
• 统计每个样本的核心与附属基因数量。
• 通过 200 次随机基因组加入顺序估计 pan/core 累积曲线。
• 计算 gene presence/absence 的 Jaccard 距离。
• 用平均连接法构建附属基因组聚类树。
• 汇总功能类别在不同泛基因组组分中的分布。

输出表格：

results/tables/pangenome_class_summary.tsv results/tables/per_genome_gene_family_counts.tsv results/tables/pan_core_accumulation_curves.tsv results/tables/gene_prevalence_distribution.tsv results/tables/jaccard_distance_matrix.tsv results/tables/functional_class_by_pangenome_class.tsv results/tables/accessory_gene_jaccard_tree.nwk

9. 图件索引与图注

所有图件均同时输出 jpg 和 pdf：

示例图注可直接作为报告模板：

若要获得本期教程数据和代码，在后台回复：20260624

10. 参考文献

1. Tettelin H, Masignani V, Cieslewicz MJ, et al. Genome analysis of multiple pathogenic isolates ofStreptococcus agalactiae: implications for the microbial pan-genome. PNAS, 2005. PMID: 16172379. https://pubmed.ncbi.nlm.nih.gov/16172379/
2. PMC 开放全文和数据沉积说明。https://pmc.ncbi.nlm.nih.gov/articles/PMC1216834/
3. Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics, 2014. https://pubmed.ncbi.nlm.nih.gov/24642063/
4. Page AJ, Cummins CA, Hunt M, et al. Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics, 2015. https://pubmed.ncbi.nlm.nih.gov/26198102/
5. Tonkin-Hill G, MacAlasdair N, Ruis C, et al. Producing polished prokaryotic pangenomes with Panaroo. Genome Biology, 2020. https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02090-4
6. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 2018. https://pubmed.ncbi.nlm.nih.gov/30423086/
7. Bankevich A, Nurk S, Antipov D, et al. SPAdes: a new genome assembly algorithm. Journal of Computational Biology, 2012. https://pubmed.ncbi.nlm.nih.gov/22506599/
8. Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics, 2013. https://pubmed.ncbi.nlm.nih.gov/23422339/
9. Walker BJ, Abeel T, Shea T, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE, 2014. https://pubmed.ncbi.nlm.nih.gov/25409509/
10. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Research, 2015. https://pubmed.ncbi.nlm.nih.gov/25977477/
11. Schwengers O, Jelonek L, Dieckmann MA, et al. Bakta: rapid and standardized annotation of bacterial genomes. Microbial Genomics, 2021. https://doi.org/10.1099/mgen.0.000685

若要获得本期教程数据和代码，在后台回复：20260624

若我们的教程对你有所帮助，请点赞+收藏+转发，这是对我们最大的支持。

小杜的生信筆記，注于分享生物信息学相关知识和R语言绘图教程。
发表论文投稿及招聘信息请使用word格式或MarkDown格式发送到

bioinfor2025@163.com

，均为无偿。

2025已离你我而去，2026加油！！

2026年推文汇总 (点击后访问)

2025年推文汇总 (点击后访问)

2024年推文汇总 (点击后访问)

2023年推文汇总 (点击后访问)

2022年推文汇总 (点击后访问)

往期部分文章

1. 最全WGCNA教程（替换数据即可出全部结果与图形）

• WGCNA分析代码六

推荐大家购买最新的教程，若是已经购买以前WGNCA教程的同学，可以在对应教程留言，即可获得最新的教程。（注：此教程也仅基于自己理解，不仅局限于此，难免有不恰当地方，请结合自己需求，进行改动。）

2. 精美图形绘制教程

• 精美图形绘制教程
• 《R语言绘图专栏–50+图形绘制教程》

3. 转录组分析教程

• 转录组上游分析教程[零基础]

• 一个转录组上游分析流程 | Hisat2-Stringtie

• Samll RNA上游分析

4. 转录组下游分析

• 批量做差异分析及图形绘制 | 基于DESeq2差异分析

• GO和KEGG富集分析

• 单基因GSEA富集分析

• 全基因集GSEA富集分析

小杜的生信筆記，注于分享生物信息学相关知识和R语言绘图教程。