The large-scale blast score ratio (LS-BSR) pipeline: a method to rapidly compare genetic content between bacterial genomes

1Whole genome sequence (WGS) data has changed our view of bacterial relatedness and evolution. Computational analyses available for WGS data include, but are not limited to, single nucleotide polymorphism (SNP) discovery, core genome phylogenetics, and gene based comparative methods. In 2005, a BLAST score ratio (BSR) method was introduced in order to compare peptide identity from a limited number of bacterial genomes. However, the “all vs. all” implementation of this method scales poorly with a larger number of sequenced genomes.

2Here we present the Large Scale BSR method (LS-BSR) that can rapidly compare gene content of a large number of bacterial genomes. Comparable methods have been published in order to group genes into gene families, including OrthoMCL, TribeMCL, and GETHOGs. Although grouping peptides into gene families is not the primary focus of LS-BSR, the output can be parsed to identify the pan-genome structure of a species; scripts are included with LS-BSR that classify coding sequences (CDSs) into pan-genome categories based on user-defined identity thresholds.

3Pipelines have also been established to perform comprehensive pan-genome analyses, including the pan-genome analysis pipeline (PGAP), which requires specific gene annotation from GenBank and complicates the analysis of large numbers of novel genomes. PGAP also doesn’t allow for the screen of specific genes of interest against query genomes in order to identify patterns of distribution. GET_HOMOLOGUES is a recently published tool that can be used for pan-genome analyses, including the generation of dendrograms based on the presence/absence of homologous genes; by only using presence/absence based on gene homology, more distantly related gene relatedness cannot be fully investigated. The integrated toolkit for the exploration of microbial pan-genomes (ITEP) toolkit was recently published and performs similar functions to LS-BSR, including the identification of gene gain/loss at nodes of a phylogeny. ITEP relies on multiple dependencies and workflows, which are available as a pre-packaged virtual machine. The authors of ITEP report that an analysis of 200 diverse genomes would take ∼6 days on a server with 12 processors and scales quadratically with additional genomes.

4The LS-BSR method can either use a defined set of genes, or can use Prodigal to predict CDSs from a set of query genomes. When using Prodigal, all CDSs are concatenated and then de-replicated using USEARCH at a pairwise identity of 0.9 (identity threshold can be modified by the user). Each unique CDS is then translated with BioPython (www.biopython.org) and aligned against its nucleotide sequence with TBLASTN to calculate the reference bit score; if BLASTN or BLAT is invoked, the nucleotide sequences are aligned. Each query sequence is then aligned against each genome with BLAT, BLASTN, or TBLASTN and the query bit score is tabulated. The BSR value is calculated by dividing the query bit score by the reference bit score, resulting in a BSR value between 0.0 and 1.0 (values slightly higher than 1.0 have been observed due to variable bit score values obtained by TBLASTN). The results of the LS-BSR pipeline include a matrix that contains each unique CDS name and the BSR value in each genome surveyed. CDSs that have more than one significant BSR value in at least one genome are also identified in the output. A separate file is generated for CDSs where one duplicate is significantly different than the other in at least one genome; these regions could represent paralogs and may require further detailed investigation. Once the LS-BSR matrix is generated, the results can easily be visualized as a heatmap or cluster with the Multiple Experiment Viewer (MeV) or R; the heatmap represents a visual depiction of the relatedness of all peptides in the pan-genome across all genomes. The Interactive Tree Of Life project can also be used to generate heatmaps from LS-BSR output and correlate heatmap data with a provided phylogeny. A script is included with LS-BSR (compare_BSR.py) to rapidly compare CDSs between user-defined sub-groups, using a range of BSR thresholds set for CDS presence/absence. Annotation of identified CDSs can then be applied using tools including RAST and prokka (http://www.vicbioinformatics.com/software.prokka.shtml). LS-BSR source code, unit tests, and test data can be freely obtained at https://github.com/jasonsahl/LS-BSR under a GNU GPL v3 license.

5To determine the scalability of the LS-BSR method, 1,000 Escherichia coli and Shigella genomes were downloaded from Genbank; E. coli was used as a test case due to the large number of genomes deposited in Genbank. Genomes were sub-sampled at different depths (100 through 1000, sampling every 100) with a python script (https://gist.github.com/jasonsahl/115d22bfa35ac932d452) and processed with LS-BSR using 16 processors. A plot of wall time and the number of genomes processed demonstrates the scalability of the method (Fig. 1A) using three different alignment methods. To demonstrate the parallel nature of the algorithm, 100 E. coli genomes were processed with different numbers of processors. The results demonstrate decreased runtime of LS-BSR with an increase in the number of processors used (Fig. 1B).

6The LS-BSR method is an improvement on a previous BSR implementation (http://bsr.igs.umaryland.edu/) in terms of speed and ease of use. The former BSR algorithm requires peptide sequences and genomic coordinates of CDSs to run. LS-BSR only requires genome assemblies in FASTA format, which is the standard output of most genome assemblers. To test the speed differences between methods, 10 E. coli genomes (Table S1) were processed with both methods. Using the same number of processors (n = 2) on the same server, the original BSR method took ∼14 h (wall time) to complete, while the LS-BSR method, using TBLASTN, took ∼25 min to complete (wall time). Because the original BSR method is an “all vs. all” comparison and the LS-BSR method is a “one vs. all” comparison, this difference is expected to be more pronounced as the number of genomes analyzed increases.

7To demonstrate the utility of the LS-BSR pipeline, a set of 96 E. coli and Shigella genomes were processed (Table S1); these genomes are in various stages of assembly completeness and have been generated with various sequencing technologies from Sanger to Illumina. The BSR matrix was generated with TBLASTN in 2 h 34 min from a set of ∼20,000 unique CDSs using 16 processors. In addition to the LS-BSR analysis, a core genome single nucleotide polymorphism (SNP) phylogeny was inferred on 96 genomes using methods published previously; the SNP phylogeny with labels is shown in Fig. S1. Briefly, all genomes were aligned with Mugsy and the core genome was extracted from the whole genome alignment; the alignment file was then converted into a multiple sequence alignment in FASTA format. Gaps in the alignment were removed with Mothur and a phylogeny was inferred on the reduced alignment with FastTree2.

8The compare_ BSR.py script included with LS-BSR was used to identity CDS markers that are unique to specific phylogenetic clades (Fig. 2). Identified CDSs had a BSR value ⩾0.8 in targeted genomes and a BSR value <0.4 in non-targeted genomes; the gene annotation of all marker CDSs is detailed in Table S2. The conservation and distribution of all clade-specific markers was visualized by correlating the phylogeny with a heatmap of BSR values (Fig. 2). This presentation provides an easy way for the user to highlight features conserved in one or more phylogenomic clades.

9E. coli and Shigella pathogenic variants (pathovars) are delineated by the presence of genetic markers primarily present on mobile genetic elements. The conservation of these markers was used as a validation of the LS-BSR method. A representative sequence from each pathovar-specific marker (Table S2) was screened against the 96-genome test set and the BSR values (Table S3) were visualized as a heatmap (Fig. 2). The BSR matrix demonstrates that pathovar-specific genes were accurately identified in each targeted genome (Table S3, Fig. 2). For example, the ipaH3 marker was positively identified in all Shigella genomes and the Shiga toxin gene (stx2a) was conserved in the clade including O157:H7 E. coli (Fig. 2). A sub-set of these 96 E. coli genomes is included with LS-BSR as test data to characterize the conservation and distribution of pathovar specific genes.

10Finally, the BSR values were used to cluster all 96 genomes with an average linkage algorithm implemented in MeV and the structure of the resulting dendrogram was compared to the core SNP phylogeny. The BSR based clustering method incorporates both the core and accessory genome, while the SNP phylogeny relies on core genomic regions alone. A comparison of the tree structures demonstrates that while Shigella genomes share a diverse evolutionary history (Fig. 3A), they all cluster together based on gene presence and conservation (Fig. 3B). This result was also observed using a k-mer frequency method, which uses all possible k-mer values to infer a phylogeny and validates the findings of the LS-BSR pipeline. The dendrogram also differed from the core SNP phylogeny in other genomes, which could represent either assembly problems, or more likely the acquisition of accessory genomic regions that are not a product of direct descent.

11The functionality of LS-BSR was compared to recently released pan-genome software packages including GET_HOMOLOGUES, ITEP, and PGAP. A set of 11 Streptococcus pyogenes genomes was chosen for the comparative analysis, as it was also used as a test set in the PGAP publication; the comparative analysis and results are shown in Table 1. Overall, the size of the core genome was comparable between methods, with LS-BSR (BLASTN) and GET_HOMOLOGUES calculating differing core genome numbers compared to the other methods. However, small differences were expected due to differing thresholds and clustering algorithms. Based on these results, LS-BSR represent a significant improvement in terms of speed and ease of use compared to comparable methods, while having comparable utility.

12One application in comparative genomics is the analysis of the pan-genome, or the combined genome, of isolates within a species. Post matrix-building scripts are available to visualize the pan-genome of a given dataset. One script (BSR_ to_ PANGP.py) creates a matrix compatible with PanGP, for visualization of pan-genome statistics. The pan_ genome_ stats.py script provides data that can be used to visualize the conservation of CDSs at different genome depths (Fig. 4A). An additional script randomly subsamples the CDS distribution at all depths and produces data that can be plotted to visualize core genome convergence (Fig. 4B), accumulation of CDSs (Fig. 4C), and the number of unique CDSs for each genome analyzed (Fig. 4D). All analyses were conducted on a set of 100 E. coli genomes, with 100 iterations.

13The LS-BSR method can rapidly compare the gene content of a relatively large number of bacterial genomes in either draft or complete form, though with more fragmented assemblies LS-BSR is likely to perform sub-optimally. As sequence read lengths improve, assembly fragmentation should become less problematic due to more contiguous assemblies. LS-BSR can also be used to rapidly screen a collection of genomes for the conservation of known virulence factors or genetic features. By using a range of peptide relatedness, instead of a defined threshold, homologs and paralogs can also be identified for further characterization.

14LS-BSR is written in Python with many steps conducted in parallel. This allows the script to scale well from hundreds to thousands of genomes. The LS-BSR method is a major improvement on a previous BSR implementation in terms of speed, ease of use, and utility. As more WGS data from bacterial genomes become available, methods will be required to quickly compare their genetic content and perform pan-genome analyses. LS-BSR is an open-source software package to rapidly perform these comparative genomic workflows.