snp-search: simple processing, manipulation and searching of SNPs from high-throughput sequencing

Background A typical bacterial pathogen genome mapping project can identify thousands of single nucleotide polymorphisms (SNP). Interpreting SNP data is complex and it is difficult to conceptualise the data contained within the large flat files that are the typical output from most SNP calling algorithms. One solution to this problem is to construct a database that can be queried using simple commands so that SNP interrogation and output is both easy and comprehensible. Results Here we present snp-search, a tool that manages SNP data and allows for manipulation and searching of SNP data. After creation of a SNP database from a VCF file, snp-search can be used to convert the selected SNP data into FASTA sequences, construct phylogenies, look for unique SNPs, and output contextual information about each SNP. The FASTA output from snp-search is particularly useful for the generation of robust phylogenetic trees that are based on SNP differences across the conserved positions in whole genomes. Queries can be designed to answer critical genomic questions such as the association of SNPs with particular phenotypes. Conclusions snp-search is a tool that manages SNP data and outputs useful information which can be used to test important biological hypotheses.


Background
Next Generation Sequencing (NGS) technologies have facilitated large scale resequencing of genomes for the purpose of Single Nucleotide polymorphism (SNP) discovery. The thousands of discovered SNPs during such projects can cause some difficulty for biological interpretation and computational analyses.
To address this challenge, various databases and tools have been created that annotate SNPs with information derived from a reference genome and store that information within a database [1][2][3][4][5][6][7].
While these tools have excelled at the management of the SNP data, less attention has been directed towards simple processing and searching of the stored data. For example, the tool developed by Wang et al for SNP management includes limited search functions and does not allow for advanced user database querying. Some databases like that described by Dereeper et al have implemented a web tool for SNP discovery and analysis. However, it does not accept Variant Call Format (VCF) files as input and is only applicable to diploid species. Another tool is SNPpy by Mitha et al which provides a database management system that allows for SQL querying and data manipulation. The SNPpy system is dependent on a relational database (PostgreSQL) that requires considerable expertise to administer and does not accept VCF. VCF files were an initiative from the 1000 Genomes Project [8] and are an output produced by many SNP calling programs such as dbSNPs [9], Samtools mpileup [10] and GATK [11].
There is a need therefore for a tool that extracts SNPs from VCF files, stores them into a simple database and provides multiple output options for analysis. To meet such a need, we have developed snp-search, an easy to use tool for management of SNPs generated from haploid next generation sequencing data. Given a VCF file, snp-search stores the SNPs generated by the variant calling algorithm into a sqlite database. Sqlite is a relational database contained in a single file and does not require complex management software. snp-search can then be used to extract useful information from the database. For example, by running snp-search using the command line syntax, the user can extract unique SNPs for a specified set of strains; generate a SNP phylogeny and provide detailed information about each individual SNP such as if it is a synonymous or non-synonymous SNP, gene information, etc (see Table 1). The snpsearch database schema has been purposely designed to be simple so that it is easy to interpret findings from large and complex biological queries. Therefore a simple SQL join is able to extract data that will provide the power to answer questions about variation within the core genome quickly and easily. snp-search also has the ability to filter SNPs depending on the function of genes within which the SNPs are found. Furthermore there is no file size limit and no dependence on any external tool. The sole requirements are a Unix environment, SQLite, Ruby 1.8.7 or greater and optionally FastTree 2 [12] for phylogenetic tree construction. As far as we are aware, there are no available tools with similar features as snp-search in the literature.

Implementation
snp-search is written in Ruby, a popular scripting language, and uses the Ruby ActiveRecord library that maps database tables to Ruby objects. It will run on most modern Unix-based architectures and can be installed with the command gem install snp-search.
snp-search has two fundamental features: 1. Creation of a local SQLite database and schema To create the database, data is sourced from two input files: reference genome (in genbank or EMBL format) and a VCF file. snp-search requires the reference genome used in the mapping process to import the features (genes, etc) and annotations into the database and requires the VCF file to populate the database tables with the SNP information. The database schema created by snp-search is designed to facilitate construction of simple queries that will address complex biological questions. Annotations from the genbank file are populated in the annotations Thus queries that relate strains to snps or annotations may be easily composed. These relationships are schematised in Figure 1.

Output requested objects from database
Once the database has been populated, snp-search provides several filtering and output options ( Figure 2).

Output formats Concatenated FASTA file
SNPs for each sample are collected from the database and concatenated and converted to FASTA format.

Tab-delimited tabular file
Information for individual SNPs are provided in tabdelimited tabular format.

Newick-file format
snp-search uses FastTree 2 [12] to generate a Newick file for SNP phylogeny.

Filtering options Filter SNPs according to the SNP Quality
A Phred-scaled score that the SNP exists at the given site (default 90).

Filter SNPs according to the Genotype Quality
A Phred-scaled score that the genotype is true (also known as Genotype quality score, default 30).

Filter SNPs according to the AD ratio
The ratio of the unfiltered count of all reads that carried that specific allele compared to other REF and ALT alleles in that site (default 0.9).

Ignore SNPs from feature
Ignoring particular SNP calls in features, such as phages, transposases, and insertion sequences to remove horizontally transferred DNA as a result of a recombination event.

Ignore SNPs in specified range
SNPs in a defined specified range will be ignored.

Ignore samples
Certain samples can be excluded from the SNP output.

Ignore non-informative SNPs
A SNP that is found in all samples is not included in the output.

Analysis
Unique SNPs for a set of strains Given a set of defined strains as input, snp-search queries the database for unique SNPs shared by the defined strains and outputs the results in concatenated FASTA or tabular format.

SNP phylogeny
Concatenated FASTA files are used to generate robust phylogenetic trees that are based on SNP differences across whole core genomes.

SNP annotation
snp-search processes the data in the database and outputs information on each SNP. The following are the description of this output: -Synonymous or non-synonymous SNP. This is calculated by translating the coding sequence and report if the SNP has whether caused an amino acid change or not.  -Original amino acid. Provide the amino acid for the SNP region in the non-mutated sequence. -Has amino acid has changed due to the mutation?
An answer of 'yes' will be given if there are at least one amino acid changes between the original nonmutated sequence and the mutated sequence.
-Is there a change in hydrophobicity of the amino acid? The answer to this question is based on the predefined set of hydrophobic and non-hydrophobic amino acids taken from Livingstone et al., [13]:  -Is there a change in the polarity of the amino acid?
The answer to this question is based on the predefined set of polarised and non-polarised amino acids taken from Hausman et al., [14]:

Database query
SQLite database can be interrogated by user defined queries (requires some knowledge of SQL). To view, administer and query the database, one may download a SQL GUI tool or simply use the Unix version of sqlite3.

Application
The snp-search tool was used to generate a SNP database of 200 Group A Streptococcus pyogenes strains (GAS) sequenced using 50 bp reads on the Illumina HiSeq platform. The reads were mapped onto the MGAS315 core genome using the bowtie tool. Samtools mpileup [10] was used for SNP discovery. The resulting VCF file and the GenBank file from the MGAS315 reference genome were used by snp-search to generate the SNP database. snpsearch -output -all_or_filtered_snps was used to generate an 'alignment' of the concatenated SNPs from each strain in FASTA format. These SNPs represent positions found in positions conserved in all the genomes, sometimes described as the core genome.
To generate a phylogenetic tree, snp-search -tree option was used that runs the phylogenetic tree algorithm FastTree 2 to generate a maximum-likelihood (ML) phylogenetic tree of the 200 genomes based on the SNP alignment from the core genome ( Figure 3A). One immediate observation from this tree was a long branch composed of a few strains within one of the clades (black circles in Figure 1). This raised questions as to whether certain SNPs were producing unrealistic branch lengths either due to false SNPs or because the SNPs are not caused by neutral evolution. Snp-search was re-run using the -ignore_snps_on_annotation phage,transposase,transposon,insertion which queries the SQL database to select only those SNPs not found in genes in mobile elements and likely to undergo recombination. These SNPs were used to remake a multiple alignment file which was then processed by FastTree 2 to produce a ML tree and viewed using Dendroscope [15]. The tree ( Figure 3B) shows that the unusually long branches to the strains represented by black circles are shortened so these strains group more consistently with the other strains from the nearest clade. A detailed comparison of the two phylogenies shows that the removal of SNPs in phage-related genes resulted in other more minor changes in the tree. A further 95 GAS strains were sequenced and added to the existing 200 strains. All 295 sequences were mapped to the MGAS315 reference genome, variants called and snp-search was used to generate the tree. MLST data extracted in silico showed almost all strains were differentiated into three Sequence Types (STs): ST15, ST315 and ST406 ( Figure 4). As shown in Figure 4, ST315 and ST406 strains were confined to one clearly differentiated clonal lineage each, whereas the ST15 strains were present in five differentiated lineages. To evaluate the utility of SNPs for the assessment of genetic diversity between the different STs, we used snp-search to find the SNPs shared only within each clade.
As shown in Figure 4, ST315 clade has 7 unique SNPs, ST406 has 9 and ST15 clades has 1,13 and 33 unique SNPs. The snp-search output for ST315 clade strains is shown in Table 1. For each unique SNP, snp-search identifies whether it is in a coding region or not, and if it is then it reports whether the SNP is synonymous or non-synonymous. Additional information such as the gene and amino acid information is also given.

Commands
Snp-search is implemented as an easy-to-use commandline tool. The user does not need experience in Unix, a programming language, or SQL (but to manually query the database the user does need to know some SQL). Figure 4 SNP phylogeny of 295 GAS strains. snp-search was used to generate the concatenated FASTA file and FastTree was used to generate the tree. Strains are separated according to the number of SNP differences between them. snp-search was also used to output the number of unique SNPs for each distinguishable clade by providing the name of the samples for each clade.
The user can use the example test sets of plasmid samples (suitable for testing purposes as they are small samples) provided in the github repository. (https://github. com/phe-bioinformatics/snpsearch/tree/master/test_data) to test snp-search functionalities. Two files are provided in the example test set directory: plasmid_test_set.vcf and Reference_file_plasmid_test_set.gbk. The following are some commands that can be used when using snp-search (the required fields are in bold): Conclusions SNP data is now commonly generated, at least for smaller genomes, by mapping whole genome sequencing reads to a reference where both known and novel SNPs can be found. To facilitate and understand these findings, we have developed snp-search, a tool that manages SNP data from high throughput haploid sequence datasets and outputs useful information which in turn can be used to test important biological hypotheses. The main advantages of snp-search are: -Parses VCF file and stores SNPs in a database -Designed simple database schema but powerful in answering complex biological queries -Provides detailed information about each SNP -Multiple SNPs filtering steps provided including filtering on the function of genes within which the SNPs are found. -Has no size limit -Only dependent on sqlite3 to build the database.
Sub-populations of SNPs can be extracted from the database to address specific biological questions, such as whether variation in specific genes is linked to phenotypic differences or does the exclusion of SNPs from horizontally transferred elements generate more reliable phylogenetic inferences.