ClustAGE: a tool for clustering and distribution analysis of bacterial accessory genomic elements

Gene content can vary greatly between closely related strains of bacteria and between other unicellular organisms [1, 2]. Genes within a species can be divided into a conserved core genome and a flexible accessory genome. The core genome of an organism consists of genetic sequence that is conserved among all or nearly all members of the species. Conversely, the accessory genome represents genetic material that is present in some, but not all members of the species. The total complement of genetic material within a species is known as the pangenome [3]. Among bacteria, modification of gene content can arise from one of three major mechanisms: gene loss, gene gain through duplication, and gene gain through horizontal gene transfer (HGT) [4, 5]. Horizontally transferred genetic elements can include such structures as plasmids, integrative and conjugative elements (ICEs), replacement islands, prophages and phage-like elements, transposons, insertion sequences and integrons [6–8]. Collectively, these horizontally transferred elements, as well as any contiguous stretch of genetic material that is not part of the conserved core genome, regardless of source or structure, can be referred to as accessory genomic elements (AGEs).

The accessory genome of bacteria can be an important source of phenotypic diversity [9]. Genes within the accessory genome can drive environmental niche adaptation or pathogenesis within hosts [10, 11]. For instance, in Pseudomonas aeruginosa, genes within the accessory genome have been found to allow the organisms to persist in environments containing heavy metals and toxic organic compounds that would otherwise be unsuitable for P. aeruginosa habitation [12, 13]. In Staphylococcus aureus, the S. aureus pathogenicity islands (SaPIs) are a class of mobile genetic elements that carry genes encoding such virulence factors as TSST1, a toxin important in toxic shock syndrome, or other toxins [14]. Antibiotic resistance genes are frequently found in the accessory genomes of clinically important bacterial pathogens. One example is the carbapenemase-carrying plasmids in Klebsiella pneumoniae and other Gram-negative pathogens that contribute to the spread of this phenotype in healthcare settings [15, 16]. Given that bacterial accessory genomes are known to be enriched in virulence factors [17], directed study of the accessory genome contents and distributions within a population could yield new diagnostics and therapeutics for bacterial infections.

Often AGEs are not discrete structures with well-defined borders and gene compositions, but instead can be mosaic and fragmented with insertions of other elements, structural rearrangements, or partial deletions [18]. Mosaic islands have previously been described in E. coli [19] and Streptococcus pneumoniae [20]. In Pseudomonas aeruginosa, a genomic island containing the type 3 secretion system effector gene exoU was found to have extensive homology and synteny of genes in this island with genes in other P. aeruginosa islands PAPI-1 and pKLC102 [21]. Given the possibly mosaic nature of accessory genomic regions, accessory element characterization is often not as simple as screening genomes for a discrete set of defined genomic islands or other horizontally transferred elements. Therefore, a robust analysis of the pan-accessory genome of a set of bacterial strains must be able to account for potential changes in structure and composition of accessory regions between strains.

With the increase in availability and affordability of whole-genome sequencing, large-scale genomic analyses of populations of isolates have become more feasible. Software packages, such as mga [22] Mauve [23], and Mugsy [24], have been developed to perform segmented alignments of complete genomes for the purposes of aligning shared genomic regions. There are also several bioinformatic tools that exist to characterize the core and pangenome of bacterial species [25–28], but few that specifically examine the accessory genome fraction. To address the accessory genome of bacteria more directly, the previously presented bioinformatic tools Spine and AGEnt [29] were developed to identify the conserved nucleotide core genome sequence in a set of genomic sequences and use this core genome sequence to perform in silico subtractive hybridization to isolate the accessory genome component of each strain. However, software such as Spine and AGEnt or others [30] that characterize the accessory genome of bacterial strains do not focus directly on providing the distribution of accessory elements in a study population. Such distribution analyses are important for answering questions about the contributions of horizontally transferred or subgroup-specific genetic elements that may contribute, for example, to a particular phenotype of interest or to understanding particular niche adaptations.

This report describes ClustAGE, a software package that clusters accessory genomic elements identified by Spine and AGEnt from a set of genomes into discrete AGE units to define the distribution of accessory elements among the analyzed genomes. Several software tools such as DomClust [31], GCQuery [32], PanOCT [33] and OrthoDB [34] have been developed for the purposes of clustering gene sequences into orthologous groups. These programs identify clusters of related genes across bacterial genomes based on gene sequences. The approach to accessory genome characterization taken by ClustAGE differs from these other approaches in that ClustAGE compares the complete nucleotide sequences of the accessory genome rather than just the protein-coding sequences. A nucleotide-sequence-based, gene-agnostic approach offers several advantages in characterizing AGE distributions. First, the identification of shared accessory elements does not depend on annotation techniques, which may differ in technique and results between strains available from public databases or collaborators. Second, intergenic sequence distribution can be studied, allowing distributions of non-protein-coding sequences such as promoter sequences or small RNAs with potential biological relevance in the accessory genome of the population to more easily be analyzed. Third, this approach has the potential to capture variable regions within otherwise conserved genes that may have arisen by homologous recombination or other mechanisms. The data generated by this software allow detailed analysis of the flexible portion of a population’s pangenome.

ClustAGE is a command-line tool built using the Perl programming language for the purpose of analyzing and comparing accessory genomic elements (AGEs) between genomes. The source code is distributed as freeware under the GNU General Public License version 3. The core functionality of ClustAGE requires BLAST+ v2.3.0 [35, 36], of which binaries for OS X or Linux 64-bit are included with the distributions. Optional features require installation of the freeware programs gnuplot v5.0 (http://www.gnuplot.info/) and/or bwa v0.7.13 [37].

ClustAGE takes as input sets of AGE nucleotide sequences previously identified from the genome sequences of at least two separate organisms. These AGE sequence sets can be extracted from complete or draft whole-genome nucleotide sequences using the previously-developed software tools Spine and AGEnt [29]. The ClustAGE algorithm identifies representative contiguous AGEs within the input data set and delineates the distribution of discrete AGEs among the strains’ genomes. It consists of two steps: defining “bins” and defining “subelements” (Fig. 1).

In the first step of this process, clustering of similar AGE sequences into “bins” is performed. First, AGE sequences from all genomes input into ClustAGE are pooled together to create a single nucleotide BLAST database. AGE sequences are then sorted by size. In the initial iteration of the clustering algorithm, the longest contiguous AGE in the dataset is chosen as a bin representative. This bin sequence is then used as the query sequence in a blastn alignment against the database of all input AGE sequences. Alignment results are filtered to remove any hits against AGEs from the same genome as the bin representative, as well as hits below user-defined sequence identity and length cutoffs. BLAST hits against AGE sequences that pass these filters are binned with the representative sequence and removed from the pool of potential bin representative sequences in subsequent iterations. Conversely, all non-aligning AGE sequences remain in the pool of potential bin representative sequences. If only part of an AGE sequence aligns to the bin representative, the non-aligning portion of the AGE sequence is isolated and added to the pool of potential bin representatives. Subsequent iterations of clustering select the next-longest complete or partial AGE sequence that was not previously binned with a larger bin representative sequence and uses it as the query sequence for alignment against the AGE sequence pool. Clustering iterations continue in this fashion until no bin representative sequences above a user-defined length threshold remain in the pool.

Once AGE bins are defined, they are further subdivided into discrete units referred to as “subelements”. Bins are divided into subelements between positions on the reference AGE where the set of input genomes aligning to the reference AGE at the base or bases before the division differs from the set of genomes aligning to the reference AGE base or bases after the position (Fig. 2). In other words, a subelement represents the longest stretch of nucleotide sequence within the bin representative that is contiguous in all strains that contain it. By dividing AGEs into discrete subelements, insertions and deletions contributing to the mosaic nature of genomic islands and other horizontally transferred elements can be identified [38, 39].

Output files from the core function of ClustAGE described above include nucleotide sequences of the bin representative and nucleotide sequences of AGE subelements longer than a user-defined cutoff. A file listing positions within the input sequences from which the bin representative AGEs were derived, as well as a file listing the positions of subelements within each AGE and the distributions of each subelement among the input sequences are also output. Optionally, ClustAGE can produce plots of AGE distributions among the input genomes for each of the bin representative AGEs (Fig. 2). This functionality requires gnuplot (http://www.gnuplot.info) to produce the plots.

ClustAGE allows users to include coordinates and descriptions of protein coding sequences (CDS) within accessory elements as input. If provided, information about which coding sequences are contained within bins and subelements is output for each AGE for which annotations in the bin reference sequence were given. If graphical output was requested, annotated gene positions and directionality will be shown in the images (Fig. 2).

One limitation of working with draft genome sequences generated by de novo assembly of short sequencing reads such as those produced by Illumina sequencing technology is that the assembly process can fail to assemble small portions of the genome even when sufficient reads covering these regions are present in the read data set. This in turn can lead to the false presumption that an AGE is absent from a genome when in fact it simply failed to be properly assembled. To try to account for data missing from de novo generated draft genome sequences assembled from Illumina reads, ClustAGE includes an option to identify missing AGE sequences from raw read sequences. After the set of AGEs is identified from accessory genome sequences as detailed above, whole-genome Illumina sequencing data provided to ClustAGE is aligned to the bin reference AGE sequences using the ‘mem’ function of bwa aligner with default settings [37]. To try to minimize false-positive alignments of core genome read sequences to accessory regions, a core genome nucleotide sequence, such as output by Spine [29], can be provided to ClustAGE. Any reads aligning to both the core genome sequence and an AGE bin sequence will be excluded. Reads aligning to AGEs above a user-defined minimum depth of coverage and producing a contiguous alignment exceeding a minimum user-defined sequence similarity will be added to the binned sequence for that genome. Alignment data from Illumina reads are only added in AGE regions that were not found to have alignments against a bin representative AGE in the original input draft sequence for a genome. To minimize false-positive results, read alignment data are also not added unless the alignment region is either at one or both of the bin representative AGE ends or contiguous with accessory genomic sequence previously aligned by BLAST from assembly data. Subelements are then redefined using the added read alignment data and a separate set of “read-corrected” subelement sequence and coordinate files are output. If optional plotting of AGE distributions was chosen, read-aligned AGE regions are plotted using a different color to distinguish them from AGE alignments derived from accessory genomic sequences (Fig. 2).

The ClustAGE results can be used to visualize and compare relative similarity of total accessory genome content among strains in the population studied. The pipeline script subelements_to_tree.pl is provided with ClustAGE for this purpose. The program quantifies relative amount of shared subelement accessory genomic sequence for each pair of genomes by calculating Bray-Curtis distances [40]. Briefly, the Bray-Curtis distance for a pair of genomes is calculated as d = 1 – (2 S _ij / (S _ii  + S _jj )) where S _ij is the total length, in bases, of subelements identified by ClustAGE in both genomes i and j and S _ii and S _jj are the total accessory genome subelement sizes, in bases, of genomes i and j, respectively. In order to cluster strains by total accessory genome similarity, a matrix of Bray-Curtis distances for each pair of input strains is used to create a neighbor-joining tree using the ‘neighbor’ function of PHYLIP version 3.696 [41]. Optional bootstrap trees from random re-samplings of the data can be generated using PHYLIP’s ‘seqboot’ and ‘neighbor’ functions. Bootstrap support values can then be calculated for each branch of the neighbor-joining tree using the CompareToBootstrap.pl script developed by Morgan N. Price (http://microbesonline.org/fasttree/treecmp.html). In addition to the neighbor-joining tree, a matrix of Bray-Curtis similarity values (1 – d) is output, as well as a file that can be used to add a heatmap of Bray-Curtis similarity values to the neighbor-joining tree in the online tree visualization software Interactive Tree Of Life (http://itol.embl.de) [42].

A utility for visualizing ClustAGE results as a pan-accessory genome figure is also available. ClustAGE Plot (http://vfsmspineagent.fsm.northwestern.edu/cgi-bin/clustage_plot.cgi) uses CGView [43] to produce a representation of ClustAGE results as bins ordered largest to smallest in a circular configuration with concentric rings indicating the distributions of accessory elements for each included strain. Although designed to be flexible, user-friendly, and powerful enough for most users, visualizations with ClustAGE Plot could become less informative with larger (i.e > 100 genomes) and/or high complexity data sets. The xml-formatted file produced by ClustAGE Plot can be downloaded and used to produce higher resolution images on a user’s local version of CGView. Furthermore, the output files generated by ClustAGE provide sufficient data for further processing and can be reformatted to serve as input for other applications capable of visualizations such as R (https://www.r-project.org/), Circos (http://circos.ca/), or other 3rd party programs, depending on the users’ needs and skills.

The ClustAGE scripts and utilities are available for download at https://clustage.sourceforge.io. A web-based implementation of ClustAGE is also available at http://vfsmspineagent.fsm.northwestern.edu/cgi-bin/clustage.cgi. The web version is limited to a maximum of 15 accessory genome sequence sets and does not support read-correction of AGEs.

ClustAGE, in combination with the core and accessory genome identification packages Spine and AGEnt [29], is an easy-to-use and accurate software tool to characterize the distribution of accessory genomic elements (AGEs) within a collection of bacterial whole-genome sequences. It includes utilities for visualizing AGE distributions and comparing and classifying relative accessory genome similarity among strains in the studied population. Taken together, the analysis output provided by ClustAGE can offer researchers a powerful new tool to study the relationships of discrete strain characteristics with flexible genome content in large genomic data sets to gain insight into bacterial evolution and adaptation.

Project name: ClustAGE.

Project home page: https://sourceforge.net/projects/clustage and http://vfsmspineagent.fsm.northwestern.edu/cgi-bin/clustage.cgi.

Operating system(s): Platform independent.

Programming language: Perl.

Other requirements: Perl 5.10 or higher, BLAST+ 2.3.0 or higher. For optional functions, gnuplot 5.0 or higher, bwa 0.7.13 or higher, and/or phylip 3.695 or higher are necessary.

License: GNU GPL v3.

Any restrictions to use by non-academics: None.