Figure

Figure. Steps in the bioinformatics pipeline implemented at Utah Public Health Laboratory.

The bioinformatic pipeline developed and implemented at UPHL consists of 8 steps (Figure): 1) read quality control, 2) reference strain determination, 3) read mapping to the reference strain, 4) single-nucleotide polymorphism (SNP) and small insertion or deletion (indel) detection, 5) de novo genome assembly, 6) genome annotation, 7) phylogenetic tree construction, and 8) phylogenetic analysis. Although such processes are standard, several software solutions are available for the respective steps.

The first step in almost all WGS bioinformatics analyses is quality control of the raw sequencing data. It is important to remove poor-quality sequence data and technical sequences (i.e., adapter sequences). Highly accurate sequences are required for SNP detection, enabling the detection of actual SNPs and distinguishing from sequencing artifacts. Quality control in our pipeline is performed by using Trimmomatic (6), a multithreaded command line tool that removes adapter sequences, trims low-quality sequence from the beginning or end of a sequence, removes reads that fall below a user-defined threshold for length, and validates paired-end sequence reads.

The second step in the pipeline is reference sequence determination. To determine SNPs, a reference sequence is needed against which to compare sequencing reads. The choice of reference sequence might have a substantial effect on the number and type of SNPs that are detected, making this step important. We use Mash for reference sequence determination (7). Mash enables us to quickly compare the large set of sequencing reads generated against the reference set of 54,118 National Center for Biotechnology Information RefSeq genomes (https://www.ncbi.nlm.nih.gov/refseq) to determine nucleotide distance and relatedness (8).

Once a reference sequence is determined, the next step in the analysis pipeline is mapping the quality-controlled sequencing reads to the reference genome. We perform read mapping by using the Burrows-Wheeler Aligner (BWA) software package with the bwa-mem option (9). BWA uses a Burrow-Wheeler Transform to efficiently align sequencing reads to reference genomes allowing for gaps and mismatches. The output of BWA is the standard sequence alignment map format known as SAM, which facilitates the next step in the pipeline.

The fourth step in the pipeline uses mapping of the sequencing reads to the reference sequence to identify SNPs and indels. We perform SNP and indel determination by using SAMtools and VarScan2, which also calculate SNP frequency in the sequence data (10,11). The output of VarScan2 can be easily viewed in the Integrative Genomics Viewer, which enables the interactive viewing of large genomic datasets (12). The output file of VarScan2 can also be used in more complex downstream analyses (i.e., to build SNP matrixes and phylogenetic trees).

The quality-controlled sequencing reads are then used for de novo genome assembly in the sixth step of the pipeline. We perform de novo genome assembly on individual isolates by using the St. Petersburg genome assembler, also known as SPAdes (13,14). The SPAdes assembler has 3 modules: sequencing read error correction; SPAdes assembly; and a mismatch corrector module. The first module error corrects the quality-controlled sequencing reads by using advanced algorithms based on Hamming graphs and Bayesian subclustering. Sequencing error correction in this manner has shown to dramatically improve genome assemblies of NGS data (15). The SPAdes assembly module uses the error-corrected reads and performs the actual assembly in an iterative manner making use of de Bruijn graphs. The resulting genome assembly is then used as input for the third module, which greatly reduces the number of mismatches and small indels by using BWA and results in highly accurate contigs (contiguous sequence data made up of overlapping sequencing reads) and scaffolds (ordered and oriented contigs based on paired-end read data).

We then annotate the resulting genome assembly to identify protein-coding genes, tRNAs, and rRNAs. We use Prokka for annotation of protein-coding genes, tRNA, and rRNA on the contigs and scaffolds generated by SPAdes (16). Prokka can fully annotate a bacterial genome in approximately 10 minutes on a high-end quad-core desktop computer by making use of a suite of existing software, tools, and sequence databases, such as UniProt (17) and NCBI RefSeq (8).

We then use shared orthologous genes to construct phylogenetic trees that provide insight into the relatedness of isolates. Once multiple genomes have been annotated, we calculate the pan genome of the combined genomes by using Roary (18). The pan genome consists of the union of genes shared by genomes of interest, and Roary can compute the pan genome of 1,000 bacterial genomes on a single CPU computer in 4.5 hours (19). In addition to determining the pan genome of the genomes of interest, Roary also generates a concatenated nucleotide alignment of the pan genome, which can be used to build a phylogenetic tree of these sequences. This pan genome alignment is used as the input to RAxML for phylogenetic tree construction (20). RAxML is a program that has been designed and optimized for conducting phylogenetic analyses on large datasets by using maximum-likelihood techniques to estimate evolutionary trees from nucleic acid sequence data (21).

The last step in the pipeline is phylogenetic analyses. These analyses can detect a signature of selection on individual genes and provide knowledge about the evolutionary forces acting on the genes of the sequenced isolates. The pan genome alignment can also be used to detect signatures of selection by calculating the ratio of the number of nonsynonymous substitutions per nonsynonymous site to the number of synonymous substitutions per synonymous site. The value of this ratio is used to infer the direction and magnitude of natural selection, with values >1 implying positive selection (i.e., driving change), values <1 implying purifying selection (i.e., acting against change), and values of exactly 1 indicating neutral selection (i.e., no selection). To determine the ratios for detecting signatures of selection, we use the YN00 model (22) implemented in the PAML software package (23). The PAML results are a plain text file that can be viewed in any word processor or imported into statistical analysis software, such as R, for further analysis or plotting.

The bioinformatic pipeline we describe can be partitioned as a function of computer resources (i.e., the number of CPUs, the amount of RAM, and the amount of storage space). Typical laptop or desktop computers might only have enough power to perform the first steps in the pipeline, whereas a high-end workstation would have enough power to perform all the steps for hundreds of samples at once. In many cases, the limiting factor is how much RAM a computer has. Many of the more complex steps in the pipeline require large amounts of RAM, often more than what many laptops and desktops can hold. All the software described can easily be installed and run on a typical desktop or laptop computer (Figure). At UPHL, we performed steps 1–4 of the described analyses on bacterial isolates by using an Apple MacBook Pro laptop (Apple, Inc., Cupertino, CA, USA) with a single 3.2-GHz Intel Core i5 processor, 16 gigabytes (GB) of RAM, and 500 GB of storage space (Technical Appendix Table 2). Many PHLs might already have the computational resources needed to perform these bioinformatic analyses on a small number of samples in a reasonable amount of time. However, some basic command-line instructions would be needed to execute software. Numerous online resources, many of them free, will help novices learn the basics of the command-line interface. One such resource is the Biostar Handbook (https://www.biostarhandbook.com). This online document and e-book is an excellent resource that introduces bioinformatics and covers all of the major areas of focus in bioinformatics, including a crash course in the command-line interface.

Computers with an increased number of processing cores, more RAM, and more storage space than the typical desktop or laptop computer will allow PHLs to perform all the analyses described here as well as more advanced and computationally intensive analyses (Figure). High-end desktops are relatively inexpensive to purchase, and it might be possible to upgrade desktops a PHL already has. All the analyses we describe here were performed at UPHL on an Apple iMac equipped with a single 3.2-GHz Intel Core i5 processor, 32 GB of RAM, and 2 TB of storage space (Technical Appendix Table 2). For 10 isolates, the analyses took ≈5 days to complete. Theoretically, the number of isolates that could be analyzed can be increased to up to hundreds of isolates on a similar high-end desktop computer; however, the amount of time to perform these analyses would also increase substantially.

With a high-end Linux-based workstation (http://www.linux.org) and a network-attached storage array, several hundred genomes can be analyzed in a reasonable timeframe (Figure). At UPHL, we invested in a high-end Hewlett-Packard workstation (HP, Inc., Palo Alto, CA, USA) with four 3.0-GHz Intel Xeon processors (Intel Corp., Santa Clara, CA, USA), each 3.0 GHz with 12 processing cores; 256 GB of RAM; and a Synology network-attached storage array (Synology, Inc., Taipei, Taiwan) with 24 TB of storage (Technical Appendix Table 2). With such a system and bioinformatics personnel in place, hundreds of genomes can be generated and analyzed in 2–3 days, providing near real-time results for disease outbreak surveillance and monitoring. In addition to high-end computer hardware, experienced personnel are needed to deploy, maintain, curate, and automate bioinformatics pipelines (i.e., bioinformaticians). To take full advantage of computational resources, programs should be automated and linked together so that as data are generated by the sequencer, they are automatically added to the bioinformatics pipelines.

With NGS becoming more and more important for public health laboratories, the need for bioinformatic analyses in greatly increasing. Unfortunately, the pace of WGS implementation is far outpacing the number of bioinformaticians being hired to work in PHLs and, understandably, not all PHLs will have the need, desire, or financial capacity to hire a full-time bioinformatician. The objective of this perspective is to show that bioinformatic analyses can be performed on everything from a simple laptop to a high-end Linux workstation and the user can have little to no experience in bioinformatics or can be a full-fledged bioinformatician. As the volume of sequencing data increases, the ability to connect phenotype to genotype becomes a reality. Knowing a priori that a microorganism is likely to be resistant to antimicrobial drugs or could be a highly virulent strain would greatly improve patient outcomes, improve outbreak surveillance, and help prioritize resources to combat outbreaks. By using molecular evolutionary analyses, PHLs can investigate the evolution of antimicrobial-resistance genes to track in near real-time mutations that are linked to newly acquired resistance genes or novel mutations that result in resistance.

NGS has the potential to revolutionize public health. NGS is not only replacing PFGE, but has the potential to replace traditional culture-based testing as well. Culture-independent diagnostic testing though metagenomic sequencing and analysis has the ability to quickly identify pathogens without applying any type of selection.