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

📅 2026/6/25 14:14:06
泛基因组 | 分享一套“数据下载、质控、组装、矫正、注释到泛基因组统计与绘图“的泛基因组分析组装代码
​生信服务器活动网址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.tsvtettelin_2005_accessions.tsv中列出了用于示例重分析的 8 个基因组入口菌株accession类型说明18RS21AAJO01000000WGS masterTettelin 2005 新沉积515AAJP01000000WGS masterTettelin 2005 新沉积CJB111AAJQ01000000WGS masterTettelin 2005 新沉积COH1AAJR01000000WGS masterTettelin 2005 新沉积H36BAAJS01000000WGS masterTettelin 2005 新沉积A909CP000114complete chromosomeTettelin 2005 新沉积NEM316AL732656complete chromosome已公开数据库基因组2603V/RAE009948complete chromosome已公开数据库基因组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.sh2. 软件环境构建所有软件安装定义在env/pangenome_env.yml中。推荐使用 mamba因为生物信息依赖较多。cdpangenome_tutorialbashscript/00_create_env.sh conda activate pangenome_env该环境包含类型软件数据下载entrez-direct, ncbi-datasets-cli, sra-tools原始 reads 质控FastQC, fastp, MultiQC组装与评估SPAdes, QUAST, CheckM矫正BWA, SAMtools, Pilon注释Prokka, Bakta泛基因组Panaroo, Roary, MAFFT, FastTreeR 分析和绘图R, ggplot2, tidyverse, ape, vegan, pheatmapBakta 和 CheckM 需要额外数据库。真实项目建议把数据库放在项目目录下例如mkdir-pdata/db bakta_db download--outputdata/db/baktaexportBAKTA_DB$(pwd)/data/db/bakta/dbmkdir-pdata/db/checkm checkm data setRoot data/db/checkm3. 数据下载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.tsv3.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.tsv4. 质控、修剪与组装原始 reads 质控由script/02_qc_trim_assemble.sh完成。该脚本对每个样本执行FastQC检查原始 reads 质量。fastp去接头、低质量端修剪、过滤过短 reads。FastQC检查修剪后 reads。SPAdes--careful进行细菌基因组短读长组装。QUAST评估 contig 数量、N50、总长度、GC 等指标。MultiQC整合质控报告。bashscript/02_qc_trim_assemble.sh推荐初筛阈值指标建议阈值解释Q30 比例大于 80%低于该值说明测序质量偏差需要谨慎adapter 含量小于 5%高接头比例会影响组装reads 数每个基因组 50 倍以上覆盖度更稳妥细菌组装通常需要足够覆盖contig N50越高越好与数据质量和重复序列有关contig 数越少越好过多 contig 可能导致泛基因组假阳性总长度接近物种预期明显偏大可能污染偏小可能缺失CheckM completeness大于 95%代表基因组较完整CheckM contamination小于 5%高污染会夸大 pan-genome5. 基因组矫正与过滤组装后使用 reads 回贴到 contig通过 Pilon 矫正小 indel、替换和局部错误。bashscript/03_polish_and_correct.sh脚本步骤BWA index为 contig 建索引。BWA MEM将修剪后 reads 比对回组装结果。SAMtools sort/index整理 BAM。Pilon根据 read alignment 矫正组装。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常见结果文件文件来源用途gene_presence_absence.RtabPanaroo/RoaryR 绘图和统计矩阵gene_presence_absence.csvPanaroo/Roary含注释的宽表core_gene_alignment.alnPanaroo/Roary核心基因系统发育summary_statistics.txtPanaroo/Roarypan/core/shell/cloud 摘要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.R8.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.sh8.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.nwk9. 图件索引与图注所有图件均同时输出 jpg 和 pdf图号文件源数据绘图脚本图意Figure 1results/figures/Figure1_pangenome_composition.jpg/.pdfpangenome_class_summary.tsvR/03_plot_figures.R核心、软核心、壳层和云基因数量Figure 2results/figures/Figure2_pan_core_accumulation.jpg/.pdfpan_core_accumulation_curves.tsvR/03_plot_figures.Rpan-genome 和 core-genome 随基因组数量变化Figure 3results/figures/Figure3_gene_prevalence_distribution.jpg/.pdfgene_prevalence_distribution.tsvR/03_plot_figures.Rgene family 在样本中的流行度分布Figure 4results/figures/Figure4_per_genome_gene_content.jpg/.pdfper_genome_gene_family_counts.tsvR/03_plot_figures.R单个基因组的核心和附属基因数量Figure 5results/figures/Figure5_accessory_jaccard_clustering.jpg/.pdfjaccard_distance_matrix.tsvR/03_plot_figures.R附属基因组 Jaccard 距离聚类Figure 6results/figures/Figure6_functional_class_distribution.jpg/.pdffunctional_class_by_pangenome_class.tsvR/03_plot_figures.R功能类别在不同泛基因组组分中的分布示例图注可直接作为报告模板若要获得本期教程数据和代码在后台回复2026062410. 参考文献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/PMC 开放全文和数据沉积说明。https://pmc.ncbi.nlm.nih.gov/articles/PMC1216834/Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics, 2014. https://pubmed.ncbi.nlm.nih.gov/24642063/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/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-4Chen 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/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/Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics, 2013. https://pubmed.ncbi.nlm.nih.gov/23422339/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/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/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格式发送到bioinfor2025163.com均为无偿。2025已离你我而去2026加油2026年推文汇总 (点击后访问)2025年推文汇总 (点击后访问)2024年推文汇总 (点击后访问)2023年推文汇总 (点击后访问)2022年推文汇总 (点击后访问)往期部分文章1. 最全WGCNA教程替换数据即可出全部结果与图形WGCNA分析代码六推荐大家购买最新的教程若是已经购买以前WGNCA教程的同学可以在对应教程留言即可获得最新的教程。注此教程也仅基于自己理解不仅局限于此难免有不恰当地方请结合自己需求进行改动。2. 精美图形绘制教程精美图形绘制教程《R语言绘图专栏–50图形绘制教程》3. 转录组分析教程转录组上游分析教程[零基础]一个转录组上游分析流程 | Hisat2-StringtieSamll RNA上游分析4. 转录组下游分析批量做差异分析及图形绘制 | 基于DESeq2差异分析GO和KEGG富集分析单基因GSEA富集分析全基因集GSEA富集分析小杜的生信筆記注于分享生物信息学相关知识和R语言绘图教程。