NASP: an accurate, rapid method for the identification of SNPs in WGS datasets that supports flexible input and output formats

Test datasets.

To demonstrate the speed and functionality of the NASP pipeline, and to compare the output with other pipelines, three datasets were selected. The first includes a set of 21 genome assemblies of members of the genera Escherichia and Shigella used in other comparative studies (Bertels et al., 2014; Touchon et al., 2009) (Table S2). REALPHY was run on self-generated single-ended simulated reads, 100 bp in length. Additional pipelines were run with paired-end reads generated by ART chocolate cherry cake (Huang et al., 2012), using the following parameters: -l 100 -f 20 -p -ss HS25 -m 300 -s 50; this method was not run in conjunction with REALPHY, as the short-read generation is integrated into the method. Unless otherwise noted, the reference genome for SNP comparisons was Escherichia coli K-12 MG1655 (NC_000913) (Blattner et al., 1997). All computations were performed on a single node, 16-core server with 48 Gb of available RAM. For kSNP, the optimum k value was selected by the KChooser script included with the repository.

The second dataset includes a set of 15 Yersinia pestis genomes from North America (Table S3). For those external SNP pipelines that only support raw reads, simulated reads were generated from genome assemblies with ART. A set of SNPs (Table S4) has previously been characterized on these genomes with wet-bench methods (unpublished). This set was chosen to determine how many verified SNPs could be identified by different SNP pipelines. All computations were performed on a single node, 16-core server with 48Gb of available RAM.

The last dataset includes simulated data from Y. pestis. Reads and assemblies from 133 Y. pestis genomes (Cui et al., 2013) were downloaded from public databases and processed with NASP using the Colorado 92 (CO92) genome as the reference to produce a reference phylogeny for WGS data simulation. Assemblies and reads were simulated from this reference phylogeny and a reference genome (CO92 chromosome) using TreeToReads (https://github.com/snacktavish/TreeToReads), introducing 3501 mutations; in this process, mutations are introduced into genomes to reproduce the phylogeny, although the mutations are completely manufactured. A phylogeny was inferred from the concatenated SNP alignment (3501 simulated SNPs produced by TreeToReads) with RAxML v8 (Stamatakis, 2014) to provide a ‘true’ phylogeny for the simulated data. Simulated reads (250 bp) and assemblies were both processed with pipelines to identify how many of these introduced SNPs could be identified.

To test the scalability of NASP on genome assemblies, a set of 3520 E. coli genomes was selected (Table S5). Genomes were randomly selected with a python script (https://gist.github.com/jasonsahl/990d2c56c23bb5c2909d) at various levels (100–1000) and processed with NASP. In this case, NASP was run on multiple nodes across a 31-node high performance computing (HPC) cluster at Northern Arizona University. The elapsed time was reported only for the step where aligned files are compiled into the resulting matrix. Time required for the other processes is dependent on the input file type and the amount of available resources on a HPC cluster.

External SNP pipelines.

Multiple SNP pipelines, both reference-dependent and reference-independent, were compared with NASP, including kSNP v3.9.1 (Gardner et al., 2015), ISG v0.16.10–3 (Sahl et al., 2015a), Parsnp v1.2 (Treangen et al., 2014), REALPHY v112 (Bertels et al., 2014), SPANDx v3.1 (Sarovich & Price, 2014), Mugsy v1r2.2 (Angiuoli & Salzberg, 2011), lyve-SET v1.1.6 (Katz et al., 2013), and the CFSAN SNP pipeline (https://github.com/CFSAN-Biostatistics/snp-pipeline). Exact commands used to run each method are shown in Supplemental Data File 1. An overview of all tested methods is shown in Table 1. Most of the methods output FASTA or nexus files, which were used to infer phylogenies. For Mugsy, the MAF file was converted to FASTA with methods described previously (Sahl et al., 2011).

Phylogenetics.

Phylogenies were inferred using a maximum likelihood algorithm implemented in RAxML v8.1.7 (Stamatakis, 2014), except where noted. The exact commands used to infer the phylogenies are shown in Supplemental Data File 1. Tree topologies were also compared on the same input data using compare2trees (Nye et al., 2006). Commands to infer these phylogenies using FastTree2 v2.1.7 SSE3 (Price et al., 2010), ExaBayes v1.4.1 (Aberer et al., 2014), and Parsimonator v1.0.2 (github.com/stamatak/Parsimonator-1.0.2) are shown in Supplemental Data File 1.

Dendrogram of multiple methods.

To visually represent the performance of different methods, a dendrogram was generated. Each phylogeny was compared against a maximum likelihood phylogeny inferred from the reference test set with compare2trees and a congruence score was calculated. A unweighted pair group method with arithmetic mean (UPGMA) dendrogram was then calculated with Phylip v3.6 (Felsenstein, 2005) on the resulting similarity matrix.

Pipeline functionality and post-matrix scripts

NASP is a reference-dependent pipeline that can incorporate both raw reads and assemblies in the SNP discovery process; NASP was not developed for the identification and annotation of short insertions/deletions (indels). NASP can use multiple aligners and SNP callers to identify SNPs and the consensus calls can be calculated across all methods. A complete workflow of the NASP method is shown in Fig. 1. Several post-matrix scripts are included with NASP in order to convert between file formats, including generating input files for downstream pipelines [e.g. Plink (Renteria et al., 2013)]. An additional script can annotate a NASP SNP matrix using SnpEff (Cingolani et al., 2012) to provide functional information for each SNP.

NASP run time scalability

To visualize how NASP scales on processing genome assemblies, a set of 3520 E. coli genomes was sampled at 100-genome intervals and processed with NASP with 10 replicates. The results demonstrate that the matrix building step in NASP scales linearly up to three and a half minutes with the processing of additional genomes (Fig. 2a). The memory footprint of this step also scales linearly (Fig. 2b) and doesn’t exceed 4Gb on a large set of genomes (n = 1000). If raw reads are used, additional time is required for the alignment and SNP calling methods, and the overall wall time would scale with the number of reads that needed to be processed.

Pipeline comparisons on E. coli genomes data set

To test differences between multiple pipelines, a set of 21 genomes of members of the genera Escherichia and Shigella used in other comparative genomics studies (Bertels et al., 2014; Touchon et al., 2009) were downloaded and processed with Parsnp, SPANDx, kSNPv3, ISG, REALPHY, CFSAN, lyve-SET, Mugsy, and NASP. For methods that do not support genome assemblies, paired-end reads were simulated with ART, while single-end reads were used by REALPHY, as this method is integrated into the pipeline.

To identify how well the simulated paired end reads represent the finished genomes, a NASP run was conducted on a combination of completed genome assemblies as well as simulated raw reads. The phylogeny demonstrates that assemblies and raw reads fall into identical locations (Fig. S1, available in the online Supplementry Material), suggesting that the paired-end reads are representative of the finished genome assemblies.

The developers of REALPHY assert that their analysis of this dataset demonstrates the utility of using their approach to avoid biases in the use of a single reference genome by using multiple references (Bertels et al., 2014). To test differences between methods, SNPs were identified with multiple reference-dependent and -independent methods, and maximum likelihood (ML) phylogenies were compared with compare2trees. The results demonstrate that all methods, with the exception of kSNPv3 and lyve-SET, returned a phylogeny with the same topology as the published phylogeny (Bertels et al., 2014) (compare2trees topological score =100 %) (Table 2). The run wall time demonstrates that most other methods were significantly faster than REALPHY (Table 2), even when REALPHY was invoked using a single reference. Wall time comparisons between methods are somewhat problematic, as some pipelines infer phylogenies and others, including NASP, do not. Additionally, using raw reads is generally expected to be slower than using a draft or finished genome assembly. Finally, some methods are optimized for job management systems, whereas others were designed to run on a single node.

One of the other assertions of the REALPHY developers is that phylogenies reconstructed using an alignment of concatenated SNPs are unreliable (Bertels et al., 2014; Touchon et al., 2009), especially with regards to branch length biases (Leache et al., 2015). However, the phylogeny inferred from a NASP alignment of monomorphic and polymorphic sites was in complete agreement with the topology of the phylogeny inferred from a concatenation of SNPs (compare2trees topological score =100 %); tree lengths were indeed variable with use of these two different input types using the same substitution model (Fig. S2). We also employed an ascertainment bias correction (Lewis correction) (Leache et al., 2015) implemented in RaxML, in order to correct for the use of only polymorphic sites, and found no difference between tree topologies using substitution models that did not employ this correction (data not shown). For this dataset of genomic assemblies, there appears to be no effect of using a concatenation of polymorphic sites on the resulting tree topology, although branch lengths were affected compared with an alignment including monomorphic sites.

To understand how the choice of the reference affects the analysis, NASP was also run using E. coli genome assemblies and simulated reads against the outgroup, E. fergusonii, as the reference. The results demonstrate that the same tree topology was obtained by using a different, and much more divergent, reference (compare2trees topology score =100 %). However, in both cases, fewer SNPs were identified by using a divergent reference (Table 2).

Some researchers suggest that reference-independent approaches are less biased and more reliable than reference dependent-approaches (Gardner & Hall, 2013). For the case of this E. coli dataset, the phylogeny inferred by Mugsy, a reference-independent approach, was in topological agreement with other reference-dependent approaches (Table 2). In fact, kSNPv3 was one of the only methods that returned a topology that was inconsistent with all other methods (Table 2); an inconsistent kSNP phylogeny has also been reported in the analysis of other datasets (Pettengill et al., 2014). To analyze this further, we identified SNPs (n = 826) from the NASP run using simulated paired-end reads that were uniquely shared on a branch of the phylogeny that defines a monophyletic lineage (Fig. S3). We then calculated how many of these SNPs were identified by all methods and found widely variable results (Table 2). Using kSNP with only core genome SNPs identified only five of these SNPs, which explains the differences in tree topologies.

In many cases, the same tree topology was returned even though the number of identified SNPs differed dramatically (Table 2). This result could be due to multiple factors, including if and how duplicates are filtered from the reference genome or other genome assemblies. With regards to NASP, erroneous SNPs called in genome assemblies are likely to be artifacts from the whole-genome alignments using NUCmer. The default value for aligning through poorly scoring regions before breaking an alignment in NUCmer is 200, potentially introducing spurious SNPs into the alignment, especially in misassembled regions in draft genome assemblies. By changing this value to 20, the same tree topology was obtained, although many fewer SNPs (n = approximately 100 000) were identified (Table 2). This value is easily altered in NASP and should be appropriately tuned based on the inherent expected diversity in the chosen dataset. Additional investigation is required to verify that SNPs in divergent regions are not being lost by changing this parameter. Another option is to use simulated reads from the genome assemblies in the SNP identification process.

Phylogeny differences for the same dataset

Previously, it has been demonstrated that different phylogenies can be obtained for the same dataset using either RAxML or FastTree2 (Pettengill et al., 2014). To test this result across multiple phylogenetic inference methods, the NASP E. coli read dataset was used. Phylogenies were inferred using a maximum likelihood method in RAxML, a maximum parsimony method implemented in Parsimonator, a minimum evolution method in FastTree2, and a Bayesian method implemented in Exabayes (Aberer et al., 2014). The results demonstrate variability in the placement of one genome (UMN026) depending on the method. FastTree2 and Exabayes agreed on the topology, including 100 % congruence of the replicate trees. The maximum-likelihood and maximum-parsimony phylogenies were slightly different (Fig. S3) and included low bootstrap support values at the variable node. The correct placement of UMN026 is unknown and is likely to be confounded by the extensive recombination observed in E. coli (Dykhuizen & Green, 1991).

Pipeline comparisons on a well characterized dataset

To test the functionality of different SNP calling pipelines, a set of 15 finished Y. pestis genomes were processed with NASP. This set of genomes was selected because 26 SNPs in the dataset have been verified by wet-bench methods (Table S4). Additionally, 13 known errors in the reference genome, Y. pestis CO92 (Parkhill et al., 2001), have been identified (Table S4) and should consistently be identified in SNP discovery methods. The small number of SNPs in the dataset requires accurate SNP identification to resolve the phylogenetic relationships of these genomes.

The results demonstrate differences in the total number of SNPs called between methods (Table 3). Most of the methods identified all 13 known sequencing errors in CO92, although Parsnp, REALPHY and kSNPv3 failed to do so. The number of verified SNPs identified also varied between methods, from 21 in kSNPv3 to all 26 in multiple methods (Table 3). An analysis of wet-bench-validated SNPs (n = 9) that are identified in more than one genome demonstrated that some methods failed to identify all of these SNPs, which could lead to a very different phylogeny and incorrect resolution of important phylogenetic relationships. In fact, such SNPs could represent critical markers, resulting in the inappropriate linkage or separation of strains in an outbreak event.

Pipeline comparisons on a simulated set of assemblies and reads

Simulated data for Y. pestis were used to compare SNP identification between pipelines. In this method, 3501 mutations (Supplemental Data File 2) were inserted into genomes based on a published phylogeny and FASTA file. Raw reads were also simulated from these artificially mutated assemblies with ART to generate paired end sequences. Reads and assemblies were run across all pipelines, where applicable.

The results demonstrate that NASP identified all of the inserted SNPs using raw reads, although 67 SNPs failed the proportion filter (0.90) and 232 SNPs fell in duplicated regions (Table 4); some of the duplicated SNPs would also fail the proportion filter. Of all other methods, only ISG identified all inserted mutations. Parsnp identified the majority of the mutations, although duplicate regions appear to have also been aligned.

To understand how the SNPs called would affect the overall tree topology, a phylogeny was inferred for each set of SNPs with RAxML. A similarity matrix was made for each method based on the topological score compared with the ML phylogeny inferred from the known mutations. The UPGMA dendrogram demonstrates that the NASP results generally return a phylogeny that is more representative of the ‘true’ phylogeny than other methods (Fig. 3). Without removing SNPs found in duplicated regions, the NASP phylogeny was identical to the phylogeny inferred from the known SNPs.

Comparisons between short-read aligners and SNP callers

One of the benefits of NASP is that it implements multiple SNP callers and aligners. To identify potential differences between methods, short reads simulated from the E. coli dataset were aligned against K-12 MG1655. Simulated reads were aligned against the reference genome with BWA-MEM, bowtie2 and Novoalign, and SNPs were called with the UnifiedGenotyper method in GATK. The results demonstrate that clear differences were observed between SNP sets between aligners (Fig. S4a). This demonstrates that aligners using different algorithms and alignment stringencies can differ dramatically in the set of SNPs called. To determine the effect of including a divergent outlier, the analysis was repeated excluding E. fergusonii. The results demonstrate that many more consensus SNPs were identified (79 %) between all three methods by excluding E. fergusonii compared with including the outlier (45 %) (Fig. S4a). These results indicate that some aligners are better at aligning divergent sequences using default parameters.

To test the effect of different SNP callers on the set of SNPs identified, BWA-MEM was used to align reads and SNPs were called with the UnifiedGenotyper method in GATK v2.7.2, Samtools v0.1.19 and VarScan v2.3.6. The results of this analysis demonstrated variability in the set of SNPs identified, although the variation was much less than observed between aligners (Fig. S4b). In this dataset, these small differences between SNP callers are unlikely to affect the tree topology.

To test the impact of different read aligners in a smaller dataset, reads from the Y. pestis North America dataset were aligned against CO92 and the number of called SNPs were identified. In this case, six SNPs were called that were only identified by bowtie2 and three SNPs were not identified by Novoalign. In outbreak situations where all SNPs may be needed to understand relationships, the union of SNPs called by all methods could be used. In situations where a large number of SNPs define the population structure, the intersection of all aligners and SNP callers provides high confidence, consistent calls.