PubMLST functionality is provided by the BIGSdb web application, which is written in object-oriented Perl and Javascript utilizing a PostgreSQL backend ^{6, 46} . It runs under the Apache web server software on Linux. The API is written as a Perl Dancer2 application interfacing with the same program libraries as the web application, run using the Starman high-performance PSGI/Plack web server ²⁷ .

The PubMLST website can be accessed from any platform using a modern web browser supporting Javascript. The underlying BIGSdb software can be locally installed on a Linux machine with as little as 4GB RAM and a single processor, but at least 4 processor cores and 16 GB RAM are recommended.

Data access. Both manual and automated access to data are available. Most end-users interact with PubMLST using the web interface ( https://pubmlst.org), through which all functionality is available. Querying of sequences can be performed by either pasting in data or uploading files for analysis via web forms. Complex queries can be constructed by combining search elements on a user-modifiable form and results further analysed using integrated plugins. Context-sensitive help is available within the interface linking to specific sections within an online users’ guide ( http://bigsdb.readthedocs.io). All data hosted within PubMLST is also accessible via the API, which enables machine-to-machine direct data exchange without user input ²⁷ . Most methods are discoverable from the root entry point ( http://rest.pubmlst.org/) and are fully documented ( http://bigsdb.readthedocs.io/en/latest/rest.html). The API is frequently used to synchronize allele and profile definitions (nomenclature), but it also supports authenticated data submissions to curator queues, isolate dataset querying and extraction of typing designations from uploaded genome sequences. This latter functionality facilitates incorporating PubMLST allele calling or use of species identification functionality from the rMLST database, directly in to an analysis pipeline ( Supplementary File 1).

Annotation. BIGSdb works exclusively with assembled nucleotide sequences, although these need not be single contiguous sequences and can range from single gene sequences through MLST datasets and draft genome assemblies to complete, finished genomes present as single sequences. Therefore, genomic data present in the sequence read archives must be assembled before deposition into the database. Once these sequences have been deposited, the identification of genes at specific locations (loci) within the deposited sequence data and their precise allelic variant, is central to the BIGSdb approach to annotation. This process is referred to as allele calling and is divided into three stages.

The first stage is the identification and cataloguing of the loci (i.e . genes) present within the group of organisms that are being analysed in a given database. A variety of means of identifying genes within assembled bacterial genomes are available, such as PROKKA ⁴⁷ . Many bacterial species now have an annotated reference genome ⁴⁸ and this can be used to seed the database with genes known to be present in that species. As the majority of bacteria have ‘open genomes’, a complete catalogue of the ‘pan genome’, that is all genes that are available to that group of bacteria ⁴⁹ , requires a continual process of annotation of novel genomes as they are deposited, with an ever-expanding catalogue of genes known to be present in that group of organisms. To this end, BIGSdb has no limit on the number of loci that can be stored in a database other than the storage capacity of the host computer. BIGSdb indexes each gene, which is normally synonymous with one locus in a bacterial genome, with a unique identifier that can be associated with any number of names that have been previously used. This is essential as individual annotations in bacteria from the same or related species often have incompatible names. For example, in the PubMLST Neisseria database a NEIS number is assigned in order of discovery to each unique gene in the database, which is linked to other names that have been used. As novel genes are identified by the on-going process of population annotation, new NEIS numbers are defined. In this way, the pan genome of the organisms in question (in this case the genus Neisseria as the PubMLST Neisseria database is genome wide) is catalogued, and a universal gene and locus nomenclature established and maintained.

The second stage is identifying the presence and location of the known genes within deposited genome data (locus identification). BIGSdb employs different strategies for this, depending on the number of loci being annotated at a given time. The most straightforward method performs a BLAST ⁵⁰ query of the genome sequence against a database of all known alleles for each gene in turn. This is satisfactory if the results for only a few genes are being indexed, but if analysing more loci, for example for a pan-genome or cgMLST analysis where perhaps 1500 or more loci are being catalogued simultaneously, this is a bottleneck that is increasingly time-consuming as more alleles are defined. Therefore, for routine calling, the process is facilitated by reducing the BLAST search space using ‘exemplar alleles’. These are defined such that every known allele for a given gene is within 10% sequence identity of an exemplar allele of the same length. The BLAST query is then performed against a database of exemplar alleles for all loci together, which identifies the location of the loci and allows the individual locus sequences to be extracted and their allelic identities to be determined by a database lookup. This method is much more efficient and allows PubMLST to scan more than 1000 bacterial genomes for a 2000 loci cgMLST scheme within one hour on a single server using 32 cores.

The third stage of population annotation is identifying and defining the alleles present at each locus. Most of this process is automated but it requires curator oversight if a new sequence is too different from an existing allele. Again, two different strategies are employed: (i) using all existing alleles in the search; or (ii) using exemplar or type alleles to constrain the search space and prevent definition of new alleles that are sufficiently divergent to the original type allele to warrant the definition of a new gene. The first strategy is more sensitive for use when developing schemes and for research purposes whereas the second is likely to be employed routinely for schemes used for surveillance or public health purposes. For cgMLST or pan-genome automated allele assignment, new allele sequences must be within 98% identity and 98% total length of a known allele, contain terminal in-frame start and stop codons and have no internal stop codons. With type alleles employed, the allowed percentage identity difference will be increased.

As each database comprises multiple genetic loci, the BIGSdb software has the capacity to organise these loci into any number of ‘schemes’. A scheme is a group of genes (loci) that are grouped together for a particular purpose. Any locus can, in principle, be included in any scheme and there can be any number of schemes of any size. Schemes can be: (i) purely ‘typological’, for example conventional seven-locus MLST schemes ⁵¹ , or (ii) directed to a particular biological or biochemical function, for example the Entner-Doudoroff pathway scheme developed for Campylobacter species ⁵² ; or combine both functional and typing roles, such as the ribosomal MLST (rMLST) scheme ¹⁸ . Schemes enable hierarchical and functional analyses with easy access to defined genetic variation, providing a facile way of extracting specific sequence diversity information from large datasets rapidly and conveniently.

Core genome MLST profile definitions and clustering. The bacterial core genome was originally conceived as those functional genes present in every member of a given group of related microorganisms ⁵³ ; however, for practical purposes it is desirable to use a more relaxed definition. This is because ‘essential’ core genes may be absent from a given WGS data set for several reasons including: (i) the organism from which it is generated is a rare mutant lacking a gene normally present in all members of the species; and (ii) technical issues due to incomplete genome assembly. Both reasons would lead to an ever-decreasing core genome, as isolates without the full complement of ‘core’ genes accumulate. Consequently, bacterial cgMLST schemes, which need to be stable if they are to form the basis of nomenclature, commonly include genes present in 95% or more of WGSs ^{25, 39, 54, 55} . Most isolates analysed with a given cgMLST scheme will therefore have some missing loci. To accommodate this, cgMLST profiles are usually defined with some missing loci, i.e . where a locus position is marked with an ‘N’ instead of an allele number. In pairwise comparisons, this will match to any allele, effectively removing the locus from that comparison, so the number of allowed missing loci needs to be kept at a low level in order to maintain resolution. The cgMLST schemes hosted on PubMLST commonly allow profile definition with up to 50 missing loci. Loci that are frequently absent from many assemblies, possibly due to them containing regions of low complexity longer than the sequence read length, are generally removed from a scheme to minimize missing data.

Once allelic profiles have been defined from WGS data, they can be clustered to identify groups of similar isolates that are likely to share a common ancestor. BIGSdb supports clustering of cgMLST schemes using a single-linkage model with multiple defined thresholds of allelic differences. These thresholds are chosen empirically and depend on the organism and the number of loci used in the scheme, but may range from 200 or more, to identify major lineages, down to 5 or fewer, suitable for identifying point source outbreaks.

Comparative genomics. The BIGSdb platform incorporates a comparative genomics plugin called Genome Comparator. This performs rapid gene-by-gene pairwise comparisons of up to 1000 genomes using either: (i) loci defined in the database, generating MLST profiles consisting of the chosen loci, which can be any defined scheme or user-selected collection of loci; or (ii) using an annotated reference genome, or simply a FASTA file of sequences, as a source of comparator sequences, producing an ad hoc whole genome MLST analysis. A pairwise distance matrix generated from these profiles is then used to generate a NeighborNet analysis ⁵⁶ using the SplitsTree software package ⁵⁷ to visualise the relationships among the isolates being analysed. This is a rapid process which is suitable for identifying clusters of related organisms, for example in a disease outbreak scenario, or to find shared genetic variants when performing functional studies. Alternatively, with the alignment option selected, aligned sequences for each locus are concatenated to produce whole genome coding sequence alignments to be used for analysis in third party applications to investigate phylogenetic relationships among more distantly related isolates ⁵⁸ . As well as using genomes already in the database, it is also possible to include uploaded private genomes for analysis in context with public datasets.

Minimum-spanning trees can be constructed from allelic profiles generated from any set of loci using a plugin that exports to an integrated GrapeTree implementation ⁵⁹ . GrapeTree was specifically designed to visualise tens of thousands of whole genome MLST (wgMLST) profiles and reconstruct relationships despite missing data. Metadata, including provenance or scheme fields (such as ST or clonal complex) can be included in the analysis and used to colour nodes, producing publication quality output.

Links to third-party analysis tools. PubMLST provides rich, structured datasets that can be combined with other data. Other web services provide visualisation tools and have implemented APIs so that data can be programmatically uploaded to them and various BIGSdb plugins are now available that can select datasets, link to appropriate metadata fields, and send them to these sites for analysis ( Figure 5). As an alternative to the integrated GrapeTree implementation, minimum spanning trees can be analysed using a plugin that uploads to PHYLOViZ Online ⁶⁰ . Phylogenetic trees can be generated from concatenated aligned sequences from any selected loci or set of loci defined by a scheme, and visualised with metadata overlays in Interactive Tree of Life ⁶¹ or combined with geographical and temporal data for visualisation in the Microreact software ⁶² .

PubMLST has been used to identify sequence variants for molecular typing since its inception. When most sequencing was still performed on individual genes by Sanger methods, it was necessary to assemble forward and reverse trace files and then compare each locus assembly against known variants. This can still be done on a per-locus basis and if an exact match is not found then the most similar allele is identified along with a list of nucleotide differences, which can be manually checked to see whether the allele is novel. With whole genome data, the rapid identification and classification of microbes is considerably more streamlined. A locally assembled genome sequence can be queried against all loci of interest rapidly ⁸ . This can be performed via the web interface by either pasting assembled contigs in to a search box, or by choosing to upload a FASTA file containing the contigs ( Figure 6A). Any defined scheme, such as MLST, can be selected from a drop-down box, and the analysis run by clicking the submit button ( Figure 6B). The analysis usually takes about a second to perform. In the case of MLST, individual allelic matches will be identified along with the ST if the combination of alleles has been previously defined. Alternatively, this analysis can be performed via the API allowing it to be incorporated directly in to a local analysis pipeline ( Supplementary File 1).

Bacterial populations often comprise distinct lineages of related genotypes that exhibit specific phenotypic characteristics that persist over time. These are shaped by mutation, horizontal genetic transfer and selection from interactions with, for example, the host immune system. Understanding the population structure is therefore important for epidemiology and public health intervention as well as for addressing more fundamental questions concerning evolution, persistence, and adaptation. The population structure of a microorganism can be represented using a cgMLST scheme and the integrated GrapeTree plugin to generate a minimum-spanning tree that can be overlaid with any metadata stored in the database. This provides an intuitive means to investigate and present the clustering of genotypes and how phenotypic or provenance characteristics map to these clusters. For example, performing a query in the Neisseria database for any isolate record where the species is ‘ Neisseria meningitidis’ with a sequence bin containing >2 Mbp data (indicating a complete genome sequence), identifies just over 12,000 records. By selecting the GrapeTree plugin, and analysis by cgMLST, an interactive minimum-spanning tree can be produced with nodes coloured by user-selectable criteria such as country, year of isolation or defined clonal complex designation ( Figure 7). Since the tree can be built using any combination of loci or schemes, and annotated with any set of metadata or field values determined by scheme (combinations of loci), datasets can be readily explored and relationships between genotype and phenotype identified. Analysis of a dataset of this size takes approximately 30 minutes to perform.

Many vaccines target expressed proteins that exhibit natural sequence variation within the bacterial population as a result of interactions with the host immune system. Over time, this may result in vaccines becoming less effective as immune pressure selects against strains with cross-reactive antigens and new lineages or variants take their place. It is therefore important that ongoing surveillance is performed to detect such changes in the bacterial population. The BIGSdb platform is flexible in how loci are defined and it allows small antigenic peptide sequences, such as those expressed on surface-exposed loops of proteins, to be used in addition to the more commonly used nucleotide sequences of complete genes. Schemes can, therefore, be defined that include only the antigen sequences of the protein components of a vaccine formulation, and variants of these defined as for any other sequence. One such scheme is the meningococcal Bexsero Antigen Sequence Type (BAST) scheme ⁴⁰ that is being used to survey structured datasets, including carriage studies, in order to provide early warning of changes in the population of meningococci that might result in reduced vaccine efficacy ⁶³ .

The ability to plot the provenance of an isolate on a geographical map and relate this to its genotypic placement in a phylogenetic tree can be useful when investigating structuring of outbreaks or the global spread of clones. Microreact is a web tool developed to produce these visualisations and the site provides a means for data to be automatically uploaded from other resources, such as PubMLST ⁶² . This can be demonstrated using data from a recent paper describing the development of a MLST scheme for Dichelobacter nodosus, the causative agent of ovine footrot, that indicated that the global bacterial population was geographically structured ³⁶ . All isolates that included a genome record in the D. nodosus PubMLST database (n=171) were selected and the Microreact analysis run using cgMLST loci. This generated a concatenated genome-wide alignment from which a Neighbor-joining tree was reconstructed automatically and uploaded to the Microreact web site ⁶² along with associated metadata. The geographical structuring of clades can be seen in the analysis ( Figure 8).