Pitfalls of bacterial pan-genome analysis approaches: a case study of Mycobacterium tuberculosis and two less clonal bacterial species

Abstract

Summary: Pan-genome analysis is a fundamental tool for studying bacterial genome evolution; however, the variety in methods used to define and measure the pan-genome poses challenges to the interpretation and reliability of results. Using Mycobacterium tuberculosis, a clonally evolving bacterium with a small accessory genome, as a model system, we systematically evaluated sources of variability in pan-genome estimates. Our analysis revealed that differences in assembly type (short-read versus hybrid), annotation pipeline, and pan-genome software, significantly impact predictions of core and accessory genome size. Extending our analysis to two additional bacterial species, Escherichia coli and Staphylococcus aureus, we observed consistent tool-dependent biases but species-specific patterns in pan-genome variability. Our findings highlight the importance of integrating nucleotide- and protein-level analyses to improve the reliability and reproducibility of pan-genome studies across diverse bacterial populations.

Availability and implementation: Panqc is freely available under an MIT license at https://github.com/maxgmarin/panqc.

Figure 1.

Summary of dataset of 151 complete Mtb genomes. Left: Maximum likelihood phylogeny of all 151 genomes, colored according to their lineage (L1-6, L8). Right: Heatmap of pairwise ANI. Below: Distribution of pairwise ANI values, and corresponding heatmap colorbar.

Figure 2.

Characteristics of Mtb SV pan-genome graph. (A) Left: Circle representing the high-level view of the Mtb SV pan-genome graph. Right: Two bubble regions shown in detail. Bubble Region 20 is representative of regions with a simple insertion/deletion, containing a single SV node (186 bp) in gene pe4 (Rv0160c). Bubble region 309 is representative of a complex bubble region, containing 88 SV nodes (55 759 bp) spanning from gene plcC (Rv2349c) to ppe40 (Rv2356c). (B) Distribution of the number of SV nodes per bubble region. (C) Distribution of SV node length. (D) Hierarchical breakdown of Core and SV nodes in specific categories of interest, showing number of nodes and cumulative length.

Figure 3.

Comparison of Mtb pan-genome predictions across different analysis parameters. (A) Comparison of the number of core and accessory genes estimated for the identical population of 151 Mtb isolates across all tested parameters: Assembly type (hybrid versus short-read), annotation pipeline (Bakta versus PGAP), and pan-genome software (Panaroo, Roary, PPanGGolin, and Pangene). Each data point represents a different set of gene clustering parameters of the specific software. (B) Number of predicted CDS features annotated by Bakta and PGAP across all hybrid Mtb genomes. (C) Number of predicted pseudogene features annotated by Bakta and PGAP across all hybrid Mtb genomes.

Figure 4.

Pan-genome tool comparison across three different bacterial species (A–C) Core and accessory genome estimates for Mtb, Eco, and Sau datasets across all tested parameters: Assembly type (hybrid versus short-read), and pan-genome software (Panaroo, Roary, PPanGGolin, and Pangene). Each data point represents a different set of gene clustering parameters of the specific software. (D) Percentage of gene absences due to CDS annotation discrepancy across Mtb, Eco, Sau. Each data point represents a different set of gene clustering parameters for Panaroo or Roary.

Figure 5.

Overview of the panqc nucleotide correction pipeline and panqc adjustment of Mtb and Eco pan-genome estimates. (A) Diagram of the panqc algorithm: In Step1, all predicted gene absences making up the predicted accessory genome are identified. In Step 2, each absent gene’s nucleotide sequence is aligned against all genomes. In Step 3, alignments are analyzed to identify if the nucleotide sequence is still present despite the previously predicted absence. In Step 4, all genes are clustered based on the similarity of their nucleotide sequences. In Step 5, pan-genome estimates are readjusted accounting for presence/absence of nucleotide sequence. (B) Comparison of Panaroo and Roary pan-genome predictions before and after panqc re-adjustment with default parameters for Mtb and Eco datasets, for both hybrid and short-read assemblies. Each data point represents a different set of gene clustering parameters for Panaroo or Roary before or after panqc adjustment.