Staphylococcus aureus viewed from the perspective of 40,000+ genomes

Staphopia analysis pipeline

The Staphopia analysis pipeline (StAP) processed FASTQ (pFASTQ) files from a single genome through quality control steps and bioinformatic analysis software. StAP (DOI 10.5281/zenodo.1255310) consisted of custom Python3 scripts and open source software organized by the Nextflow (Di Tommaso et al., 2017) (v0.28.2) workflow management platform (Fig. 2). When available we used BioConda (Grüning et al., 2017) to install the open source software. Summary statistics of the original input and subsequent downstream results files were collected at each step of the pipeline. For portability, StAP was wrapped in a Docker container. The version of the pipeline used in this work was Docker Image Tag: 112017(https://hub.docker.com/r/rpetit3/staphopia/).

The input to StAP was either single- or paired-end FASTQ file (or files) generated from Illumina technology. StAP contained an option that allowed FASTQ data to be pulled from the ENA based on the experiment accession number (ena-dl v0.1, https://github.com/rpetit3/ena-dl). A MD5 hash (md5sum) was generated from the input FASTQ data and cross-referenced against a list generated from processed genomes to prevent reanalysis of the same input. BBduk (Bushnell, 2016) (v37.66) was used to filter out adapters associated with Illumina sequencing and trim reads based on quality. Read errors were corrected using SPAdes (Bankevich et al., 2012) (v3.11.1). Based on the corrected reads, low quality reads were filtered out and the total dataset was subsampled to a maximum of 281 Mbases (100× coverage of the N315 reference chromosome (Kuroda et al., 2001)) with Illumina-cleanup (v0.3, https://github.com/rpetit3/illumina-cleanup/). This file (or files, if paired-end) we termed “pFASTQ.”

Processed FASTQ reads were assembled de novo using SPAdes (Bankevich et al., 2012) (v3.11.1). SPAdes also marked assembles as putative plasmids based on evidence such as relative read coverage (Antipov et al., 2016). Summary statistics of the assembly were created using the assembly-summary script (https://github.com/rpetit3/assembly-summary). A BLAST+ nucleotide database was created from the assembled contigs to be used subsequently for sequence query matching. Open reading frames and their putative functions were predicted and annotated using Prokka (Seemann, 2014) (v1.12) and its default database.

The S. aureus strain N315 (Kuroda et al., 2001) chromosome (ST5 methicillin-resistant S. aureus (MRSA); accession NC_002745.2; length 2,814,816 bp) was used as a reference for calling consensus single nucleotide polymorphisms (SNPs) and indels in the pFASTQ reads using the GATK (McKenna et al., 2010) (v3.8.0) pipeline. GATK pipeline also incorporated BWA (Li & Durbin, 2009) (v0.7.17), SAMtools (Li et al., 2009) (v1.6) and Picard Tools (v2.14.1, http://broadinstitute.github.io/picard/) software. Identified variants were annotated using the vcf-annotator script (v0.4, https://github.com/rpetit3/vcf-annotator). Jellyfish (Marçais & Kingsford, 2011) (v2.2.6) was used to count k-mers of length 31 base pairs (31-mers) in the pFASTQ file. We used Ariba (Hunt et al., 2017) (v2.10.2) to make antibiotic resistance and virulence predictions for paired-end reads only. Resistance phenotypes were predicted using the MegaRes reference database (Lakin et al., 2017) and virulence using the Virulence Factor Database (Chen et al., 2016) core dataset.

Multi-locus sequence typing was determined by two or three methods depending on the whether the pFASTQ was paired- or single-end. All methods used the S. aureus MLST allele sequence database downloaded from https://pubmlst.org/saureus/ (November 2017). Alleles for each of the seven loci were aligned against the assembled genome using BLAST+ (v2.7.1+) (Camacho et al., 2009). Alleles and ST were determined based on perfect matches (100% nucleotide identity with no indels). We also used the MentaLiST (Feijao et al., 2018) (v0.1.3) software to call MLST based on k-mer matching of the alleles to the pFASTQ file. Unlike the BLAST+-based MLST method, MentaLiST did not require exact matches to alleles to predict a ST. If the pFASTQ was paired-end, Ariba (Hunt et al., 2017) (v2.10.2) also determined MLST alleles and ST. The default ST call for each genome was determined in the following order: agreement between each method, agreement between MentaLiST and Ariba, agreement between MentaLiST and BLAST+, agreement between Ariba and BLAST+, Ariba alone without a novel or uncertainty call, MentaLiST alone, and finally BLAST+ alone. cgMLST was determined with MentaLiST using the S. aureus cgMLST scheme (Leopold et al., 2014) available at http://www.cgmlst.org.

Evidence for Staphylococcal Cassette Chromosome mec (SCCmec) predictions were based on multiple approaches. The primary approach aligned SCCmec typing primers, downloaded from http://www.staphylococcus.net/, against the assembled genome using BLAST+ (v.2.7.1+) (Camacho et al., 2009). Samples with a perfect match to primer pairs for a given amplicon were assigned an SCCmec type following the Kondo et al. (2007) algorithm. Genes and proteins associated with SCCmec mere aligned against the assembled genome using BLASTN and TBLASTN. We also mapped the pFASTQ to each SCCmec cassette using BWA (Li & Durbin, 2009) (v0.7.17). The overall cassette and mec region coverage statistics were determined as well as the per-base coverage for each cassette using genomeCoverageBed (Quinlan & Hall, 2010) (v2.26.0). The methods described above were based on the 11 SCCmec types currently listed in the http://www.sccmec.org (I–XI) and hence did not include recently described types XII and XIII (Wu et al., 2015; Kaya et al., 2018). We labelled a genome as “MRSA” only if each mecA typing primer (Kondo et al., 2007) had a perfect BLASTN match on the de novo assembly, a predicted mecA gene ortholog had a BLAST score ratio of at least 95%, or Ariba (Hunt et al., 2017), as previously described, predicted reads in the paired-end pFASTQ file matching a mecA target.

Web application, Relational Database and application programming interface

We used Django (v2.0), a Python web framework, to develop a PostgreSQL (v10.1) backed relational database for storing the results from the analysis pipeline (Fig. S1). A Django application was created for each module of the pipeline, automating the creation of database tables for the results. Python scripts building off Django were developed for insertion of results from each StAP module or the StAP as a whole. A web front-end was developed (https://staphopia.emory.edu) using the Bootstrap (v4.0) and jQuery (v3.2.1) web frameworks. We used the Django REST framework to develop an extensive API that allowed users to create queries accessing multiple samples. We also developed an R package, Staphopia-R (DOI 10.5281/zenodo.1255314), to programmatically access the API. The API and its endpoints were documented to allow users to further develop their own packages in a language of their choice. The source code for our web application was made available at DOI 10.5281/zenodo.1255312.

Processing public data

We used the Cancer Genomics Cloud (CGC) Platform, powered by Seven Bridges (http://www.cancergenomicscloud.org/), to process S. aureus genomes through StAP in November 2017. CGC allows users to create custom workflows based on Docker containers, then execute these workflows on the Amazon Web Services (AWS) cloud platform. We obtained a list of publicly available S. aureus sequencing projects from the ENA web API using the following search term:

“tax_tree(1280) AND library_source=GENOMIC AND (library_strategy=OTHER

OR library_strategy=WGS OR library_strategy=WGA) AND

(library_selection=MNase OR library_selection=RANDOM OR

library_selection=unspecified OR library_selection=“size fractionation”).”

We only processed projects which used Illumina sequencing technology. CGC opened AWS r3.xlarge instances (30.5 GB RAM, four processors) that downloaded FASTQ files from the ENA, using ena-dl for each genome, and ran the StAP pipeline. Ouput files were returned to the CGC, then uploaded into the Staphopia database server.

Metadata collection

We used the ENA API to download and store any information linked to the “Experiment,” “Study,” “Run” and “BioSample” accessions into the database for each genome. We also determined each sample’s publication status using three approaches.

The first approach identified existing links between SRA, a mirror of ENA, and PubMed using NCBI’s Entrez Programming Utilities web API (Entrez Programming Utilities Help, 2010). For any links identified, we used the corresponding PubMed ID to extract information corresponding to the publication and stored them in the database.

The second approach searched for accessions within the text of scientific articles. We searched PubMed using the term, “S. aureus,” limited to the years between and including 2010 (the date of the first publicly available Illumina data upload), and 2017. The saved results, stored as XML, were then loaded into Paperpile, a subscription-based reference management tool, and the corresponding main-text PDFs were automatically downloaded. This process did not include supplementary information files, which required a manual operation. For those articles in which a PDF could not be automatically downloaded, attempts to manually acquire the PDF were made. Using the text search program “mdfind,” available on Apple OS X, each accession (BioSample, Experiment, Study and Run) in the Staphopia database was used as a separate query to search all the PDF files. Experiment accessions with a corresponding PubMed ID were then stored in the database. In cases where a Study, BioSample or Run accession was identified in PDF text, each associated Experiment accession was linked to the corresponding PubMed ID.

For the third approach, a collection of PubMed articles with primary descriptions of S. aureus genome sequencing studies was manually curated. For these studies, the PDF and all available supplementary information were downloaded. The process of text-mining the articles and linking Experiment information to PubMed ID was repeated as described in the second approach. A list of these articles is available in BibTeX format from the code repository.

Creating non-redundant S. aureus diversity set

Using available metadata, we selected a non-redundant diversity (NRD) set of genomes that were “Gold” quality (see “Results”), linked to a publication and each had a unique ST. When more than one strain from a ST was available, we randomly selected one individual giving priority to samples with collection date, site of isolation and location of isolation fields filled.

Using predicted variants against N315, we extracted a list of genes that had complete sequence coverage (i.e., “core” genes) but no predicted indels. For each reference gene sequence, we created an alternative gene sequence with SNPs predicted in each sample. If all the 31-mers of the alternative gene sequences were present in a 31-mer database of the pFASTQ file of genome created using the Jellyfish (Marçais & Kingsford, 2011) tool, the reconstructed sequence was considered “validated.” Validated reconstructed gene sequences or all genomes were stored in the database and made available through the API for rapid phylogenetic comparisons.

A total of 380 genomes were selected for phylogenetic analysis as representatives of unique STs with good quality metadata. The set of validated genes were extracted and concatenated into a single sequence for each genome and saved in multi-FASTA and PHYLIP formats. A guide tree was generated with IQ-Tree (Nguyen et al., 2015) (v8.2.11, -fast option) for identification of recombination events with ClonalFrameML (Didelot & Wilson, 2015) (v1.11). A recombination free alignment was created with maskrc-svg (https://github.com/kwongj/maskrc-svg). We used IQ-Tree to generate the final maximum likelihood tree with the GTR model and bootstrap support. Bootstrap support was generated from 1000 UFBoot2 (Hoang et al., 2018) (ultrafast bootstrap) replicates. We annotated the tree using iTOL (Letunic & Bork, 2016).

Design of the Staphopia analysis pipeline and processing 43,000+ genomes

The StAP (Fig. 2) was written to automate processing of individual S. aureus genomes from Illumina shotgun data. The pipeline was designed as a series of modules running individual software packages, organized by the Nextflow (Di Tommaso et al., 2017) workflow language, which made it possible to run the entire pipeline or individual components as needed. The first step of the pipeline was to import single- or paired-end FASTQ files either as local files, or from the ENA database. Following quality-based trimming and down selection of the FASTQ to 281 Mbases (∼100× coverage of the N315 reference chromosome (Kuroda et al., 2001), NC_002745.2), analyses were run on the raw pFASTQ files directly, or on de novo genome assemblies constructed by the SPAdes program (see “Methods” for more details). We decided to downsample the input FASTQ files for two reasons: to manage the computational burden when running thousands of genome projects and also to achieve genome datasets with consistently sized pFASTQ input files. The threshold of ∼100× coverage was chosen after preliminary studies (https://github.com/staphopia/staphopia-paper/blob/master/analysis/07-supplementary.pdf) showed that there was either small or no improvements in outcome for downstream assembly and remapping steps for input files >100× but large increases in processing time and memory requirement. We created a Postgres database to store results from the StAP analysis and a web front end and a web API for mining the data. An R package (Staphopia-R) was written for interacting with the API and was used for most analysis presented in the results.

In November 2017 there were 44,012 publicly-available shotgun sequencing projects with FASTQ files in ENA. Illumina technology was the dominant platform, accounting for 99% of samples (N = 43,972). A total of 81% (N = 35,580) of them had at least 281 Mbases sequence data. We processed all Illumina genomes in parallel through the StAP using cloud servers (see “Methods”). On parallel r3.xlarge instances with 30.5 Gb RAM and four processors, the mean time to process a genome was 52 min with an interquartile range of 47–56 min (Fig. S2).

Sequence and assembly quality trends

We identified samples that were likely mixed-samples or not S. aureus whole genome shotgun projects and/or were of low technical quality and marked them to not be included in subsequent analysis. We removed genomes that failed to match to any known allele of the seven MLST loci (323 genomes), had a total assembly size that differed by more than one Mb from a typical S. aureus chromosome (<1.8Mb or >3.8Mb; 764 genomes), or had a GC content differing more than 5% (<28% or >38%; 467 genomes) of the expected 33% GC content. Failure to complete the StAP pipeline due to poor data quality, and coverages less than 20× were flagged in 101 and 142 genomes, respectively. In total, we removed 1,023 genome projects, leaving 42,949 for further analysis.

We placed genomes into an arbitrary ranking of 1–3 (“Bronze,” “Silver” and “Gold”) based on the pFASTQ coverage and average sequencing quality. Paired-end genomes that had read lengths exceeding 100 bp, a coverage of 100× and an average per base quality score of at least 30 were given a Gold rank. The purpose of the Gold rank was to group together high-quality samples with near-identical coverage. Paired-end genomes with similar read length and quality cutoffs but a lower sequence coverage (between 50× and 100×) were classified as Silver. The remaining samples were given a rank of Bronze. Single-end reads were classified as Bronze no matter the read length, quality or coverage. More than 70% of the samples were of rank Gold (N = 31,014). There were 5,931 Silver and 6,004 Bronze rank samples. Each year since 2012, the number of Gold ranked genomes have exceeded Silver and Bronze (Fig. 3).

Changes in sequence quality and de novo genome assembly metrics over time reflected the development of Illumina technology. Mean per-base quality scores increased from ∼32 in 2010 to >35 in 2012 and have stayed at that level since. The mean sequence read length rose in steps from <50 in 2010 to ∼150 bp in 2017. Assembly metrics such as N50 (Earl et al., 2011), and mean and maximum contig length have gradually increased since 2010. Bronze ranked genome projects had similar (or sometimes even higher) mean per read quality scores than Gold and Silver since 2011. However, Silver and Gold assembly metrics such as N50 and mean contig size were generally quite similar and higher than Bronze.

Genetic diversity measured by MLST

We obtained a view into the genetic diversity of the sequenced S. aureus genomes by in silico MLST using Ariba (Hunt et al., 2017), MentaLiST (Feijao et al., 2018) (both taking pFASTQ as input, but using different algorithms) and BLASTN against assembled contigs. A ST was assigned to 42,337 (98.6%) genomes. Of these, 41,226 (97.3%) calls were in agreement between MentaLiST, BLAST+ and (if paired-end) Ariba methods; 828 had agreement between two methods and a no-call on the other, and 189 were supported by one program with no-calls from the other two. Of the remaining 612 genomes not assigned to a known ST, 306 were predicted to be in a novel ST based on matches to known alleles of each of the seven loci. The remaining 306 genomes had 1–6 known S. aureus MLST alleles.

The 42,337 genomes assigned to existing STs represented only 1,090 STs of 4,466 in the PubMLST S. aureus database (November 2017, https://pubmlst.org/saureus/). The abundance distribution was weighted toward common strains, with the top 10 (STs 22, 8, 5, 239, 398, 30, 45, 15, 36 and 105) representing 70% (N = 29,851) of the genomes (Fig. 4).

The cgMLST set of 1861 loci (November 2017, https://www.cgmlst.org/) were assigned to the genome set using MentaLiST. There were 38,677 distinct patterns, with only 1,850 patterns found in more than one sample, the remaining 36,827 patterns were represented by a single genome.

Antibiotic resistance genes

Treatment of S. aureus infections has been complicated by the evolution of strains resistant to many commonly used antibiotics (Foster, 2017). In particular, MRSA, carrying the mecA gene encoding the PBP2a protein that confers resistance to beta-lactam antibiotics, has become a global problem. We designated a genome as MRSA if each mecA typing primer (Kondo et al., 2007) had a perfect BLASTN match on the de novo assemblies (26,743 strains), a predicted mecA gene ortholog had a BLASTN score ratio of at least 95% (26,430 strains), or the Ariba (Hunt et al., 2017) algorithm predicted reads in the paired-end pFASTQ file matching a mecA target in the MegaRes (Lakin et al., 2017) database (27,120 strains). The number of genomes having at least one of these criteria (27,628) was 64% of the total number. Of these, 95% (26,340) of the samples had agreement between each of the criteria. The top five most common STs had a large portion of MRSA strains (Fig. 4), which reflects the selection bias of the research community in investigating these significant hospital and community pathogen strains over other S. aureus.

The mecA gene is usually horizontally acquired as part of a mobile genetic element called SCCmec (Katayama, Ito & Hiramatsu, 2000). SCCmec elements have been classified into at least 11 classes that vary in composition of mec genes, ccr cassette recombinase genes and spacer regions (http://www.sccmec.org). Knowledge of the SCCmec type can be useful for high-level characterization of MRSA strain types (Kaya et al., 2018). We showed that 10 of the 11 cassettes in the current schema map to at least one genome with highest coverage (an approximate method for assigning SCCmec type) (Table 1). Of the 26,462 (26,185 paired-end) genomes with at least 50% cassette coverage, 96%, 96% and 99% are MRSA based on primer BLASTN, protein BLASTN or MegaRes, respectively. All type XI cassettes were mecA negative by primer BLASTN because these contained the mecC allele (García-Álvarez et al., 2011; Shore et al., 2011), which was sufficiently different to be outside the normal distance for a positive match. Additionally, we found 53 genomes which matched to at least 50% of a SCCmec cassette but were not MRSA and had no reads mapping to the mec region of the cassette.

In addition to mecA, we found numerous other classes of non-core antibiotic resistance genes using the MegaRes (Lakin et al., 2017) class designations (Table 2). We did not consider SNPs/indels in core genes associated with resistance for this analysis. The most common class of resistance genes were beta-lactamases found in 37,758 genomes. Following this, the most common were the genes putatively conferring fosfomycin, macrolide-lincosamide-streptogramin and aminoglycosides resistance (24,205, 22,322, 17,968 genomes respectively). As with MRSA, the other common resistance genes were not distributed evenly among the top ST groups (Fig. 4), reflecting sampling ascertainment bias and also possibly differences in geographic distribution and prevalence of healthcare-isolated strains in the most common genotypes.

Publication, metadata and strain geographic distribution

One challenge to using publicly available datasets through ENA or SRA is determining whether there is a published article describing the sequenced genome. We found through NCBI’s Entrez Tools (eLink) that 6,712 genomes were linked to 48 publications in PubMed (March 2018). We attempted to add to the number by using text-mining methods to find S. aureus accession numbers in PDFs of S. aureus genome publications, ascertaining an additional 5,209 genomes in 30 publications. Therefore, of the 42,949 samples deposited between 2010 and 2017, only 28% (N = 11,921) could be linked to a publication (Fig. 1). Since many genomes have been deposited in the last 1–3 years, this reflected the often significant time lag between depositing sequence data and final publication.

We noted that collection of metadata from public sequencing projects was another challenge. When submitting genome sequences to databases only a limited number of metadata fields are required, leading to the bulk of the information needing to be extracted manually from a publication, if it can be found. Only 40% (N = 17,034) genomes had a collection date, 35% (N = 14,983) had a geographic location and 35% (N = 14,768) had isolate source metadata. Using the available geographic data to geocode the sites of collection, we found that strains were from five continents and at least 40 countries. There was a strong bias toward strains from Europe (N = 7,314) and North America (N = 5,882), reflecting where the funding for most of the early sequencing studies had originated.

A non-redundant S. aureus diversity set

The number of SNPs compared to the N315 reference strain varied from 6 to 141,893 within our collection of 42,949 genomes. The stepped pattern of the distribution (Fig. 5) reflected the organization of S. aureus into CCs. Apart from CC5 strains closely related to N315, the majority of S. aureus had ∼50–50,000 SNPs and ∼500–1500 indels called by the GATK pipeline (McKenna et al., 2010). There were a group of 240 most distant strains with >55,000 SNP (Fig. 5) that were found to be closer to the sister species, S. argenteus (Holt et al., 2011) based on ANI imputed by mash (Ondov et al., 2016), although 230 of these were assigned a S. aureus ST.

Of the 6,904 S. aureus genomes of Gold rank linked to a publication we selected a group of 380, each having a distinct ST as a NRD set of genomes. Of the 2,756 annotated N315 genes (excluding RNAs), 1,113 genes had no indels when reads from each genome in the NRD dataset were mapped. Of these, 878 were “core” genes found in every genome. We reconstructed these genes for each of the NRD genomes starting with the N315 sequence and substituting predicted SNPs. These predicted sequences were then validated by decomposing into 31-mers and cross-checking whether each k-mer was present in pFASTQ files processed by Jellyfish (Marçais & Kingsford, 2011). We concatenated the 878 genes for each member of the NRD set and created a tree based on the 44,377 variant SNP positions (Fig. S3). The structure of the unrooted species tree resembles previous S. aureus phylogenies (Planet et al., 2017). The assembled genomes and concatenated gene sequences from the NRD dataset is available at DOI 10.6084/m9.figshare.6263435.