NovelSNPer: A Fast Tool for the Identification and Characterization of Novel SNPs and InDels

Typically, next-generation resequencing projects produce large lists of variants. NovelSNPer is a software tool that permits fast and efficient processing of such output lists. In a first step, NovelSNPer determines if a variant represents a known variant or a previously unknown variant. In a second step, each variant is classified into one of 15 SNP classes or 19 InDel classes. Beside the classes used by Ensembl, we introduce POTENTIAL_START_GAINED and START_LOST as new functional classes and present a classification scheme for InDels. NovelSNPer is based upon the gene structure information stored in Ensembl. It processes two million SNPs in six hours. The tool can be used online or downloaded.


Background
One of the main goals of resequencing projects is to detect genomic variation between an individual of interest and the reference genome that is stored in public databases. Typically, the whole genome or a section of the genome sequence encompassing the genomic region of interest is chosen as reference sequence for the resequencing project. The DNA sample is manually enriched for nonrepetitive sequence in the target genomic region using capture oligomers before performing the sequencing reactions. Depending on the technology, reads of length in the range between a couple of dozens and a few hundred nucleotides result. One sequencing reaction typically yields in the order of dozens of millions of reads.
Thus far, several fast and robust programs are available to assemble the reads from an NGS experiment, to compare the resulting consensus sequence with the reference sequence and to extract genomic variations [1][2][3][4]. However, only few tools are available to date that offer support for the further characterization of identified variations. SNP annotation is a largely standardized process performed in software pipelines at the large genome centers [5,6]. Previously, Ensembl made its SNP annotation algorithm accessible via the tool "Variant Effect Predictor" [7]. SNP annotation in more specific contexts is furthermore performed by a series of other software tools [8][9][10][11][12]. However, with the rapid spreading of high-throughput sequencing machines, there is a growing need for high-quality assembly and annotation tools to cope efficiently with the massive amount of sequencing data in respect to annotation and speed data analysis and to extend the analysis to indels.
The motivation to develop NovelSNPer was to offer an easy-to-use and fast tool that allows the distinction between variations that are already classified as reference variations (e.g., SNPs that were assigned an rs-number and that were stored in dbSNP [13]) and previously unknown variations which were discovered in the sequencing data and to characterize the variations. Other tools with a similar scope that were recently developed as, for example, SNPnexus [14] or VarySysDB [15] harness a variety of external data sources to characterize polymorphisms, but are based on static datasets restricted to human and offer only limited batch processing capabilities.
Our classification scheme comprises twenty-one variation classes that were defined largely analogously to those applied by Ensembl [16] (for further detail see Table 1). The usage of NovelSNPer is not restricted to humans. NovelSNPer is applicable to lists of variations deduced from resequencing data of all species whose genome sequence is 2 Advances in Bioinformatics Table 1: List of variations. Name and description of all variation classes used in NovelSNPer. A cross (X) indicates that a variation class can be assigned to an SNP, MNP, or InDel, respectively. Generates a potential start codon in 5 UTR or in an exon of noncoding transcript Start lost X X X Deletion of an existing start codon at the beginning of coding region Delete exon X Deletion of a whole exon Merge exon X Deletion of a whole intron Acceptor X X Deletion of the exonic and intronic splice site upstream of an exon Donor X X Deletion of the exonic and intronic splice site downstream of an exon Splice site X X X Variation near the exon-intron boundary Essential splice site X X X Intronic variation within 3 bp range of an exon-intron boundary stored in the Ensembl database. The only prerequisite is that the resequencing data was assembled to a genomic reference sequence which is stored in the Ensembl database.

Implementation
2.1. Input. NovelSNPer reads in an ASCII table (text file) that contains the information on all nucleotide polymorphisms between the derived sequence and the reference sequence in six columns: Name, Chromosome, RelStart, RelEnd, AltAlleles, and RefAllele, where Name represents an identifier for the variation chosen by the user, Chromosome is the chromosome number the variation was found on, RelStart and RelEnd are the position (start and end) of the variation in the reference sequence in base pairs (relative position), AltAlleles contains the nucleotide found at that position in the sequence data, and RefAllele is the nucleotide at that position in the reference sequence.
In case of heterozygosity both alleles are specified in AltAlleles separated by a slash or in IUPAC code. Alternatively NovelSNPer reads a common variant reporting format (VCF or Personal Genome SNP format [17,18]). The start position of the reference sequence relative to the chromosome sequence and the species must be specified as parameters on the command line upon execution of NovelSNPer.
There are three modes to use NovelSNPer.
(i) NovelSNPer uses the gene annotations from the Ensembl database to characterize the variations of the input file.
(ii) NovelSNPer uses the gene annotations from a local gene file to characterize the variations. The formats GFF and BED are supported [19,20].
(iii) A genome annotation of a desired species from Ensembl is downloaded and saved at the home desktop by NovelSNPer.
An example for an input file is shown in Table 2. The command line syntax is NovelSNPer.pl-species human Input.txt.
Alternatively, the website can be used to evaluate the input file as shown in Figure 1.

Output.
As output of a NovelSNPer run two tables are produced each in text and HTML format: one table provides a basic summary of the results for each variation, the other provides more detailed information about each transcript that is affected by a variation. The basic table contains the following data: Name, Chromosome, RelPos, ChromPos, FuncClass, NearestGene, Distance, Status, and ID. The meaning of the variables is given in Section 3. This table has the same number of  entries as the input table and AltAllele  RefAllele  var1  12  12483171  12483171  C  G  var2  12  12483200  12483200  G  T  var3  12  12506184  12506184  T  -var4  12  12527105  12527108  GATT  TTAG  var5  12  12550471  12550471  A/G  T  var6  12  12588583  12588583  A  C  var7 12 12588586 12588587 -CA Figure 1: Web interface of NovelSNPer. Input file, species, genetic code, and other options can be entered into a graphical user interface at http://www2.hu-berlin.de/wikizbnutztier/software /NovelSNPer/.
and ProteinDomain. It serves the scientist, who is interested in performing a deeper analysis of individual variations above all by providing the transcript ID in which a variation was assigned a specific functional class. In the following, we will describe the variables of the detailed table that were not yet described earlier: The columns AltAllele1, AltAllele2, and AltAllele3 host the alleles that were found in the sequencing data. In case of homozygosity, only AltAllele1 will have an entry, the other two will be set to NA. In case of heterozygosity, both AltAllele1 and AltAllele2 will have an entry, but AltAllele3 will be set to NA. We provide the column AltAllele3 to account for cases of variations with three alleles that could be produced, for example, if the sequencing experiment was performed with a pool of DNA or in case of ambiguous variation calling. RefAllele is the nucleotide that is found in the reference sequence at the position of the variation. Codon, AminoAcid, and RefAminoAcid are only given in case the variation is in a coding region, otherwise NA is given. The nucleotide at the variable position in the codon is coded as IUPAC symbol. The variable Strand denotes whether the gene associated with a variation is annotated on the plus or minus strand of the reference genome. The user can choose to save the used Ensembl version and Assembly version in the file. The tables in text format can be further processed in a statistical analysis program like R [21] or with a spread sheet like Excel. In the HTML table IDs, transcripts, genes, and protein domains are linked with Ensembl allowing easy access to further data. Furthermore, a report is produced containing key statistical figures like the number of variations processed, a listing of the number of variations in each functional class and the number of cases, in which the reference nucleotide that was specified in the input file did not coincide with the nucleotide of the reference sequence.
An example for the two output files is shown in Tables  3 and 4. Additionally the output can be converted into the GFF format and then be visualized using Artemis [22] or the Integrated Genome Browser [23].

Algorithm and Working
Process. First, NovelSNPer connects to the Ensembl database and downloads the strand and the position of all transcripts and exons. Alternatively, NovelSNPer reads a GFF-file where the gene positions are annotated. Then NovelSNPer creates a hash function, which is referenced by the transcript IDs. The hash values are a list containing start, end, coding start, coding end, strand, and a hash function. This second hash function is referenced by the exons of the transcripts and its hash values are start and end of the exon.
Afterwards, the genetic code is loaded in the cache. Then the following steps are performed for each variation in the input file (see Figure 2 for the flowchart). (ii) A sequence object ("slice") of length one nucleotide is generated from the chromosomal sequence at the global position of the variation.  (iii) The reference nucleotide from the input is compared with the nucleotide obtained from the chromosomal sequence. If these do not coincide this variation is not processed and all output variables are set to NA, but it is kept in the output for reference purposes. (iv) The alternative allele sequence is compared with the reference sequence. Depending on this comparison, the variation is classified as single nucleotide polymorphism (SNP), multinucleotide polymorphism (MNP), inversion or indel. (v) If the variation is an indel, the sequence of the variation is compared with 13 bp of the flanking sequence to determine, if it is a potential transposon [24].
(vi) The nearest gene (NearestGene) and the distance (Distance) to it in base pairs are determined. The distance is given as negative number if the nearest gene is found upstream of the variation and as positive number if it is downstream.
(vii) If the sequence slice harbours a variation the variation to be examined is recognized as known variation, otherwise the variation is considered as novel. The variable Status is set to KNOWN or NOVEL, respectively. If an input variation has the same coordinates as a known variation, but its alleles do not coincide with the alleles in the database, its status is set to ALTERNATIVE.
(viii) If the distance to the nearest gene exceeds 5 kb, the variation is classified as INTERGENIC, if the distance is greater than 0 but less or equal 5 kb, the variation is classified as UPSTREAM or DOWNSTREAM, respectively. Upstream regions of this extent usually harbour the promoters [25]. The nearest gene and the distance to it are given for both novel and known variations. If the variation is within a gene (distance = 0) the following steps are performed for all its transcripts since the effect of a variation can be dependent on its context. First, it is decided if the variation is located in an exon, an intron, or if it overlaps both. In the Advances in Bioinformatics 5 second case the variation is classified as INTRONIC. If the variation is in an exon, it is determined if it is in the noncoding or in the coding part of an exon. If it is in the noncoding part of an exon, the variation is classified as 5PRIME UTR and 3PRIME UTR, depending on whether it is located before the start codon or after the stop codon. If it is in the coding part of an exon, it is determined if the variation changes the length of the protein or not. If it changes the length, it is FRAMEKEEP, respectively, FRAMESHIFT. If the variation does not change the length of the protein, it is checked if the variation entails an amino acid exchange in the resulting protein sequence or not; the variation is classified as NON SYNONYMOUS CODING or SYNONYMOUS CODING, respectively. If one of the SNP alleles is part of a stop codon the SNP is considered as STOP LOST or STOP GAINED depending on whether the reference allele or the alternative allele constitutes the stop codon. If an SNP is located in the first two or in the last two nucleotides of an intron it is classified as ESSENTIAL SPLICE SITE. If an SNP is within the first three to eight or last three to eight nucleotides of an intron or within the first or last three nucleotides of an exon, it is classified as SPLICE SITE. If an SNP is located within a mature miRNA or within a non-coding gene, it is classified as WITHIN NON CODING TRANSCRIPT.
We label an SNP as START LOST if it is a nonsynonymous SNP in a start codon. Finally, we label an SNP as POTENTIAL START GAINED if the SNP is located in the 5' UTR of a protein-coding gene or within the exon of a noncoding transcript and if the variation results in the creation of an ATG start codon.
(ix) Each variation in the coding part of an exon is checked if the corresponding amino acid is part of a Pfam protein domain [26]. Also the pattern of the sequence shows if the variation is an inserted transposon or if the variation occurs because of a transposon which already shifted away.

Start Gained/Start Lost.
To estimate how strong a variation in genotype affects the phenotype, it is helpful to look at the proteome. A variation, which leads to a change in the protein, can have a profound effect on the phenotype. There exist variations which alter the initiation code and thus lead to a possible suppression of the protein [29]. Also an altering of the protein due to a variation in the translation start is well discussed [30][31][32]. Hence, in addition to the traditional variations, NovelSNPer also searches for changes in the start codon and determines if a variation in the 5 -UTR leads to a new potential start codon.

Subclasses of Complex Indels.
Deletions can overlap the coding region and the noncoding region. Depending on their location and length they can have different functions. They can delete a whole exon, merge two exons by deleting the intron, delete the acceptor, or delete the donor. If the variation deletes a whole exon or merges two exons, it is possible to determine the alternative protein exactly. If the deletion only deletes the splice site it is not possible to calculate the alternative protein, because the emerging splice pattern is unknown. This is one reason why it is essential to subclassify the complex indels. NovelSNPer uses the self-describing classes MERGE EXONS, DELETE EXONS, ACCEPTOR, and DONOR ( Figure 3).

Evolutionary Genesis.
As an option, NovelSNPer can analyse if a given indel is a duplication of the adjacent sequence, an inversion, or a shifting element. Although it does not reveal anything about the function of the indel, it gives a hint of the evolutionary genesis of the indel and the mechanism how this indel has derived.
Also it calculates the repeat-length of an indel. A difference in the length of the indel and the period-length of the indel indicates that there were multiple mutations which led to this indel.

Genetic Code.
NovelSNPer loads the standard genetic code at startup per default. But it is optionally possible to use the genetic code of prokaryotes, yeasts, mycoplasma, or mitochondria (as given in [33][34][35]).

Reliability.
We wrote the Perl script NovelSNPer.pl that allows for easy and fast analysis of variation lists produced in a next-generation sequencing project. This script first determines if a variation identified in the sequence data corresponds to an already known variation that is stored in the Ensembl database or if it is a previously unknown variation. Unknown variations are classified into one or several out of a total of 21 functional classes. The functional class of an variation permits the scientist to judge the consequences that a variation has for gene expression or protein conformation. Variations are furthermore assigned to a gene if they are located within that gene or the distance  Figure 2: The workflow of NovelSNPer. Each variation from a list is checked if it is a previously known variation or a novel variation, then the variation is attributed to one or several out of twenty-one functional classes.
to the nearest gene is given. Variations can be assigned to two genes if these overlap on the two strands of the genome.
To test our tool, we compared the annotation attributed by NovelSNPer with that given by Ensembl (version 61). For this purpose, we tested all 30 million human variations represented in Ensembl. The classes POTENTIAL START GAINED and START LOST could not be tested because they were uniquely used by NovelSNPer.  Figure 4: Example of a multi-class SNP. SNP rs11540005 can be assigned to seven functional classes in the human NDUFV1 gene. The NDUFV1 gene seems to be related with dilated cardiomyopathy [36]. Shown are seven transcripts from the diverse transcriptional landscape at this genomic site. Graphical visualization was done with fancyGene [37]. SOAT1 by NovelSNPer); 10% of the differently-classified variations are classified as FRAMESHIFT by Ensembl while NovelSNPer classified them as SYNONYMOUS CODING or NON SYNONYMOUS CODING. An analysis of these variations revealed that the length of the alternative allele and the length of the reference allele are the same, but greater than 1 (e.g., the alternative allele and the reference allele of rs71585892 has a length of 5 bases). An SNP was considered as classified in agreement with Ensembl if the classes coincided. In cases where NovelSNPer assigned more than one class, the class given by Ensembl had to be among the classes that were assigned by NovelSNPer in order to be considered as agreement.
Interestingly, it is not an exception that an SNP is assigned to more than one functional class, but about half of the SNPs belong to more than one functional class. One SNP (rs11540005) belonged to seven classes. This can be understood considering the fact that one gene is typically made up of several transcripts or splice variants, and the role of an SNP can be dependent on the transcript variant ( Figure 4). For example, the reading frame can vary between two transcripts, and, hence, one and the same SNP can be SYNONYMOUS CODING in one transcript and NON SYNONYMOUS CODING in the other. Even within one transcript, an SNP can take on more than one functional class, for example, an SNP can be simultaneously of type SPLICE SITE and of type 3PRIME UTR. Last, but not least, there are rare cases where two genes overlap, and an SNP at such a locus occupies two different roles in these two genes.

New Classes of Variants.
Since we were not aware of any systematic genome wide analysis of start gained and start lost SNPs, we were interested to study the abundance of such cases. We found that 1,418 of 23,621 human protein coding genes (6.0%) dispose of at least one alternative start gained SNP and 218 (1.0%) of at least one start lost SNP. We furthermore identified 5,585 genes of a total of 49,509 (protein-coding and nonprotein-coding) genes Table 5: Ambiguous short indels. Multiple deletion annotation on human chromosome 12. The deleted nucleotide A is underlined in the reference sequence. In the new sequence a gap at the position of the deletion is shown.

Indels with Ambiguous Identification.
In a test run we checked all 4 million insertions and deletions of the human genome represented in Ensembl. We found over 20,000 indels with ambiguous annotation. This means several different indels in the Ensembl database are the same indel in reality. A list of some of these multiple indels is available as Additional File 4. As an example the deletions rs71918324, rs71702364, rs5796236, rs35226411, rs72397401, rs35471040, and rs61916428 on chromosome 12 between position 6551619 and 6551638 represent one and the same indel. The details of this indel can be seen in Table 5. An ambiguous deletion of length six and its 6 bp-sequence is shown in Table 6(a). An ambiguous deletion which is 3 bp long and has five different identifiers is shown in Table 6(b).
It is not sufficient to compare only the position of indels. It is not even certain to have two different indels, if the alleles are different. To decide whether two indels are the same or not the sequence between this two alleles has to be considered. Even though the percentage of multiple annotated indels is very low NovelSNPer offers an effective way to identify these misleading annotations.

Bovine Datasets.
In a second run NovelSNPer annotated 2.3 million variations of the bovine genome, which were published in [28]. The variation data-set consists of 2,307,509 SNPs and 683 indels. 639 of these indels are a duplication of the previous nucleotide and 44 indels are an insertion of a new nucleotides. In Table 7 the distribution of the variations across genomic regions is shown. The distribution into functional classes of all variations in coding regions is shown in Table 8.

Discussion
Whole genome resequencing experiments are performed to systematically identify genomic variations. For example, the complete genome of a single bos taurus animal was sequenced to identify millions of previously unknown cattle SNPs [28]. In another work, artificial mutations that are responsible for phenotypes in caenorhabditis elegans could be identified thanks to whole genome sequencing [38]. Distilling the huge quantity of information into meaningful lists of SNPs is a multi-step bioinformatics process. NovelSNPer is an easy to use tool that helps scientists with the analysis of next-generation sequencing data. Lengthy lists of SNPs from  Non-synonymous SNPs can also entail a modification of a protein's conformation depending on how dissimilar the exchanged amino acids are. Whereas these two classes of SNPs can have a more or less direct effect on a protein (see, e.g., [39]), SNPs in untranslated regions (UTRs), introns, and up-and downstream regions have the potential to alter the binding behaviour of transcription factors or splice factors and, thus, to alter gene expression, indirectly. Finally, variations in the immediate neighbourhood of exon-intron boundaries could influence the splicing process of a transcript. One important type of genomic variation are frame shifts, that is, disruptions of the reading frame. This variation is mediated by indels of one nucleotide or of a stretch of nucleotides whose length is not a multiple of three into the coding part of a gene.
Two more variation features permit life scientists to assess the severeness of the consequences a variation can have: (a) the protein domain [26] that is affected by an indel or a nonsynonymous SNP and (b) the difference between observed and expected conservation score [27]. A great difference between observed and expected conservation score suggests that there is negative selection pressure more than expected to maintain the nucleotide at that position. A variation at such a position should therefore be likely to have negative consequences for the phenotype. Protein domains can be considered as protein building blocks that occur in various proteins within a species and across species. They are more conserved than those parts of the protein that are not organized in protein domains. A non-synonymous SNP that leads to the exchange of an amino acid that is part of a protein domain should therefore have more severe consequences for the phenotype than a non-synonymous SNP that affects an amino acid outside a protein domain.
NovelSNPer is one of the few tools that explicitly classifies variations as POTENTIAL START GAINED or START LOST if they mediate between a start codon and a non-start codon in the assembly and the reference sequence. We found such a classification only in the Mouse SNP database by the Center for Genome Dynamics of the Jackson Laboratory http://cgd.jax.org/cgdsnpdb/.
A reason why other tools or databases are obviously reluctant to use these two classes might be that a start codon alone is not sufficient to predict the translation start for sure. A variation which is classified as POTENTIAL START GAINED is a potential new start codon that competes with the previously existing start codon (if this has remained unaltered). In contrast, a START LOST variation can be considered as unambiguous signal that the translation start is shifted as compared to the wild type. Alternative translation start sites were discussed, for example, in [30,31]. The figures that we presented in Section 3 suggest that the protein synthesis of quite a substantial number of genes can be altered by POTENTIAL START GAINED or START LOST variations. One percent of the genes appear to be susceptible to protein synthesis suppression altogether.
By querying the Ensembl database online, NovelSNPer does not need to be maintained regularly. Its great value is furthermore that resequencing data from all species whose genome sequence is available in Ensembl can be processed. By using NovelSNPer offline it is possible to annotate variations for all species as long as the genome-structure of the transcripts is available as BED file or GFF file [19,20].
NovelSNPer completes the functionality of existing tools in many ways: variation calls are checked against the reference sequence to guarantee the correctness of the results, the usage of relative SNP positions is possible, the novel functional classes POTENTIAL START GAINED and START LOST are introduced, heterozygous variation calls can be treated, there is no limit in the number of variations to be processed, and information about the degree of conservation and protein domains is presented. SNPs and indels can be processed alike.

Conclusions
To make full usage of the data that are generated in next-generation-sequencing experiments, the lists of called variations must be efficiently screened. NovelSNPer is a fast tool that exploits the annotations about genes and genomes from the Ensembl database to classify each called variation as novel or as previously existing. Each variation is classified into one or more of twenty-one functional variation classes. Two of these classes, START LOST and POTENTIAL START GAINED, were thus far rarely included in SNP analysis programs. However, we showed that these two types of functional classes are quite frequent and could play a significant role in protein synthesis. The great number of species that can be analysed, the structured and detailed outputs, and the integration of new features like additional variation classes or the calculation of the conservation score for each variation make NovelSNPer a versatile analysis tool. Also the possibility to predict the effect of variations on a newly discovered transcript which is not in a public database is very useful. We thus anticipate a wide range of applications in the biological, medical and agricultural sciences.

Availability and Requirements
NovelSNPer is implemented in Perl using Bioperl and Ensembl's Perl APIs. The latest version of the program can be downloaded free of charge from http://www2.hu-berlin.de /wikizbnutztier/software/NovelSNPer/. The source code is also available at the website.
NovelSNPer runs under Linux and Windows and needs at least 512 MB RAM. MySQL, Perl, BioPerl and the API modules ensembl, ensembl-compara, and ensemblvariation (application protocol interfaces http://www.ensembl.org/info/docs/api/index.html) have to be installed locally. Alternatively, the web interface can be used.