Abstract

Introduction

Exome sequencing has become a popular and effective strategy to identify variants causing Mendelian diseases or extreme phenotypes. Non-synonymous single nucleotide variants (nsSNVs) are the major candidate variants in such studies. Typically, exome sequencing will discover a large number of nsSNVs, among which many are novel (e.g. not reported in dbSNP). Researchers will rely on various functional predictions and annotations to filter and prioritize those nsSNVs to shorten the list for further (experimental) validation. To facilitate this process, we developed dbNSFP v1.0 (Liu et al. 2011).

Based on the Consensus Coding Sequence (CCDS) version 20090327 (Pruitt et al. 2009), dbNSFP v1.0 compiled a collection of 75,931,005 nsSNVs in the human genome, including both previously reported and potentially novel ones. For each nsSNV, four functional prediction scores: SIFT (Kumar et al. 2009), Polyphen2 (Adzhubei et al. 2010), LRT (Chun and Fay 2009) and MutationTaster (Schwarz et al. 2010); and one conservation score: phyloP (Siepel et al. 2006), were collected. Since the publication of dbNSFP v1.0, it has attracted much attention and been broadly used by human geneticists and sequencing centers as well. It has been recommended by the Faculty of 1000 and adopted by various software (e.g. Lindenbaum et al. 2011; Li et al. 2012; San Lucas et al. 2012; Chang and Wang 2012; Sifrim et al. 2012; Zhang et al. 2013) and databases (e.g. Li et al. 2011).

To fulfill the increasing demand for better functional annotation for SNVs discovered in exome sequencing studies, we have upgraded the dbNSFP to v2.0. Now the database is separated into two parts, dbNSFP_variant and dbNSFP_gene. As their names indicate, the former focuses on variant annotations (including prediction scores and conservation scores), and the latter focuses on gene annotations. As to variant annotation, the database has expanded its SNV collections not only based on a more up-to-date GENCODE 9 annotation (Harrow et al. 2012), but also included all potential essential splice site SNVs (ssSNVs), which are another type of candidate variants in exome sequencing studies. Its core score collection has added two more functional prediction scores: MutationAssessor (Reva et al. 2011) and FATHMM (Shihab et al. 2013); and two more conservation scores: GERP++ (Davydov et al. 2010) and SiPhy (Garber et al. 2009; Lindblad-Toh et al. 2011). To facilitate filtering common SNVs observed in human populations, allele frequencies from the 1000 Genomes Project phase 1 data (Abecasis et al. 2012) and the NHLBI Exome Sequencing Project data (Fu et al. 2013) were added. As to gene annotation, dbNSFP v2.0 has collected rich functional annotations for all genes in the database, including various IDs for different databases, functional descriptions, and gene expression and gene interaction information, among others. Details of the upgrade and preliminary analyses of the core scores are reported in the following sections.

New and Updated Functional Annotations

As CCDS might be over-conservative on gene annotation and human reference build NCBI36/hg18 has been replaced by GRCh37/hg19, we completely rebuilt the backbone SNVs of dbNSFP based on GENCODE 9. In short, at each coding site and essential splice site (defined as the first two and last two nucleotide sites of an intron) we arbitrarily “mutated” the reference allele to the other three alternative alleles and collected them into the database. The total backbone SNVs now reached 89,617,785 (including 87,347,043 nsSNVs (Table 1) and 2,270,742 ssSNVs), which is an 18% increase compared to dbNSFP v1.0. Among them, 89,572,881 contain hg18 coordinates, obtained via the liftOver tool from the UCSC Genome Browser (Meyer et al. 2012), to facilitate queries based on hg18.

Three functional predictions scores and one conservation score of dbNSFP v1.0 (SIFT, Polyphen2, MutationTaster and PhyloP) have been updated. To avoid confusion, the original scores from the algorithms (not the re-scaled scores as in v1.0) were used, except LRT, for which we believe our monotone re-scaled score is easier to interpret than its original score; and missing scores were no longer imputed (as in v1.0) and were presented as a single dot (“.”). In the dbNSFP v2.0, we collected Polyphen2 version 2.2.2 scores based on both HDIV and HVAR training sets. As the same nsSNV may have multiple HDIV (or HVAR) predictions and scores according to different amino acid positions of different transcripts of the same Uniprot (The UniProt Consortium 2011) gene, all scores were included (separated by “;”) and their orders corresponded to the Uniprot accession numbers and amino acid positions in the Uniprot_acc column and the Uniprot_aapos column, respectively. As to Polyphen2 version 2.2.2, the score thresholds separating “probably damaging”, “possibly damaging” and “benign” predictions are 0.956 and 0.453 for HDIV, and 0.908 and 0.447 for HVAR, respectively. Throughout this paper, we regard score 0.5 as the threshold for Polyphen2's binary predictions (i.e. “deleterious” versus “tolerated”).

Two new functional predictions scores, MutationAssessor and FATHMM, have been added. Precomputed scores from MutationAssessor release 2 (hg19) were downloaded from http://mutationassessor.org/. We used its functional impact combined score as our MutationAssessor score, which ranges from -5.545 to 5.975; the larger the score the more likely it will be deleterious. MutationAssessor provides four types of predictions (to be functional): high, medium, low and neutral. To form binary predictions, we treat high and medium predictions as “deleterious” and low and neutral predictions as “tolerated”. FATHMM v2.1 database were downloaded from http://fathmm.biocompute.org.uk/ and installed on our local server, and its default scores (weighted for human inherited-disease mutations with Disease Ontology) for all the backbone SNVs were then retrieved. These score ranges from -18.09 to 11.0; the smaller the score the more likely it will be deleterious. Binary prediction is also provided and the threshold separating “deleterious” and “tolerated” is -1.5. In case there are more than one FATHMM scores for the same nsSNV due to isoforms, we took the smallest score (most deleterious) as our FATHMM score.

We have also added two new conservation scores: GERP++ and SiPhy. GERP++ base-wise scores were downloaded from http://mendel.stanford.edu/SidowLab/downloads/gerp/. Although GERP RS scores were typically used to measure the conservation of a nucleotide site in Mendelian disease studies (e.g. Cooper et al. 2010), an alternative measure might be a scaled RS score with the corresponding neutral rate (NR) of the site (i.e. RS/NR ratio). Therefore, both NR and RS scores were included in the database. SiPhy scores based on 29 mammalian genomes were downloaded from http://www.broadinstitute.org/mammals/2x/siphy_hg19/. We used its logOdds scores as our SiPhy scores. We also included the SiPhy estimated stationary distribution of A, C, G and T of the site (29way_pi) in our database to facilitate alternative conservation measures. For both the RS and logOdds scores, the larger the score the more conserved the site.

Additional SNV annotations include: observed alternative allele frequencies in the 1000 Genomes Project phase 1 data; observed alternative allele frequencies in the NHLBI Exome Sequencing Project ESP6500 data set (from ANNOVAR (Wang et al. 2010)); rs numbers from a cleaned version of dbSNP build 129 (from UniSNP, http://research.nhgri.nih.gov/tools/unisnp/); ancestral allele (from the 1000 Genomes Project phase 1 data); reference amino acid (aaref); alternative amino acid (aaalt); coding sequence strand (+ or -); reference codon; SNV position in the codon (1, 2, or 3); codon degenerate type (0, 2 or 3); SLR test statistic of the codon, which is a measure of natural selection acting on the codon (Massingham and Goldman 2005); the protein domain(s) the SNV resides on, according to the InterPro database (Hunter et al. 2012); and amino acid position(s) (aapos) as to Ensembl transcript(s). If there are multiple amino acid position(s) for the same SNV, the order of the positions corresponds to the multiple Ensembl transcript IDs in the Ensembl_transcriptid column. ssSNVs will have missing (“.”) in the aaref and aaalt columns and “-1” in the aapos column.

To facilitate gene-based SNV prioritization, we have strengthened the functional annotation of genes and added the following information: gene IDs from multiple databases (HGNC (Gray et al. 2013), Uniprot, Entrez Gene (Maglott et al. 2011), CCDS, Refseq (Pruitt et al. 2012), UCSC, and MIM (Amberger et al. 2011)); pathway information, function descriptions, disease descriptions and MIM phenotype IDs from Uniprot; trait association from the GWAS catalog (Hindorff et al. 2009); eGenetics (Kelso et al. 2003) and GNF/Atlas (Su et al. 2002) gene expression data from BioMart (Guberman et al. 2011); gene interaction data from IntAct (Kerrien et al. 2011) and BioGRID (Chatr-aryamontri et al. 2012); and estimated probability that the gene is haploinsufficient (Huang et al. 2010) and recessive (MacArthur et al. 2012).

A Summary of Functional Prediction Scores and Conservation Scores

As different functional prediction scores and conservation scores are derived from different methods using different information (summarized in Table 2), they have different ranges and distributions. Evenly dividing each score into 100 bins between its minimum and maximum, Figure 1 shows the frequency of each bin collected in dbNSFP v2.0. SIFT, Polyphen2 (both HDIV and HVAR) and MutationTaster have U-shape distributions with majority of the scores close to either 0 or 1. LRT (rescaled) presents an L-shape distribution with majority of its scores close to 0. FATHMM and GERP++ show skewed unimodal distributions, while MutationAssessor, PhyloP and SiPhy show complex multimodal distributions.

Open in a separate window

Figure 1

Distributions of the prediction and conservation scores.

Understanding the correlation between the scores is important for SNV prioritization. Typically consensus prediction by multiple scores is considered more reliable than a prediction by a single score. However, some scores are more correlated than others; therefore an agreed prediction of two less correlated scores carries more weight than that of two highly correlated scores. Based on absolute Spearman's rank correlation coefficient (aRCC), correlation strengths between the scores are mostly low to moderate (0.25-0.65) (Table 3). Exceptions include high correlations between the two Polyphen2 scores (0.97) and the three conservation scores (0.78-0.90), and very low correlations between FATHMM and other scores (0.13-0.21) (Table 3). The pairwise agreement of the binary predictions between the scores ranges from 42.05% (between FATHMM and Polyphen2-HDIV) to 76.82% (between MutationTaster and LRT), when excluding the 88.97% agreement between the two Polyphen2 scores (Table 3). Figure 2 shows the UPGMA dendrogram of the scores when using 1-aRCC as a measure of distance between scores. It is worth notice that LRT and MutationTaster are more closely related to the conservation scores than other functional prediction scores. As to LRT, this observation is not surprising because it is based on the comparison of aligned coding sequences of 32 vertebrate species. As to MutationTaster, on the other hand, this observation suggests that conservation information is heavily weighted in its prediction model. Another interesting observation from Table 3 and Figure 2 is that FATHMM has low correlation with any other score.

Open in a separate window

Figure 2

UPGMA dendrogram of the prediction and conservation scores.

The numbers of missing data in the prediction scores and conservation scores per chromosome can be derived from Table 1. In general, conservation scores have low missing data rates (0.07%-1.50%), while prediction scores' missing rates range from 11.17% (SIFT) to 22.07% (LRT). The missing rates are higher than those of dbNSFP v1.0, partially due to the fact that the GENCODE annotation of a gene is less stringent than that of CCDS. Figure 3 shows the percentages of dbNSFP v2.0 entries having 0 to 6 (non-missing) prediction scores, and the percentages of those having 0 to 3 (non-missing) conservation scores.

Open in a separate window

Figure 3

Percentages of the dbNSFP v2.0 entries having 0 to 6 (non-missing) prediction scores (A), and the percentages of those having 0 to 3 (non-missing) conservation scores (B).

New Features of the Companion Search Program

The companion search program written in Java (search_dbNSFP) has also been upgraded. The human reference sequence build GRCh37/hg19 is now the default choice. To search SNVs based on hg18, users need to use the “-v hg18” option. The search program now supports vcf format for the input file. If the file name has an extension of “vcf”, the program will automatically query the database by the “chr pos ref alt” format. Alternatively, users can also search SNVs by “chr pos” and “chr pos ref alt refAA altAA” formats, by gene name or various database IDs, and by Uniprot ID/access number and protein position. By default all columns of dbNSFP2.0_variant and most columns of dbNSFP2.0_gene (except the first three columns) will be written to a user specified output file. Users can also specify the columns to output using the new “-w” option. All queries that do not have a match in the database will be written to an “.err” file. More details can be found in the readme file of the companion search program included in the database distribution package.

Suggested Usage

The main purpose of the dbNSFP is to facilitate the SNV filtering/prioritizing step in exome-sequencing based Mendelian disease studies. To filter/prioritize SNVs based on the contents of the database, we have the following suggestions: (1) filter in (i.e. retain) nsSNVs and ssSNVs from a list of SNVs discovered through sequencing by searching the dbNSFP using the companion search program; (2) filter out common nsSNVs and ssSNVs by removing any SNV with a MAF larger than a curtain threshold (e.g. 5%) in any of the 1000 Genomes Project populations and NHLBI Exome Sequencing Project populations; (3) prioritize (i.e. rank) SNVs by the number of prediction scores supporting a “deleterious” prediction (from all Polyphen2 scores, pick one score, e.g. the most deleterious one, to represent the Polyphen2 score); (4) prioritize SNVs by a conservation score (since the three conservation scores in the dbNSFP v2.0 are highly correlated); (5) filter out SNVs in genes that are not expressed in the disease related tissue(s); (6) highlight SNVs in genes that cause related MIM diseases or are associated with related phenotypes based on GWAS; (7) highlight SNVs in genes that interact with disease causing genes; and (8) highlight SNVs in genes that are in disease related pathways. At last, we want to give some friendly warnings to the users of dbNSFP or other functional prediction tools: all methods have their own limitations; do not blindly trust any single method; using consensus prediction or majority vote might be a good practice in general but not a silver bullet.

Acknowledgments

References