VaRank: a simple and powerful tool for ranking genetic variants

Input and output

To run VaRank, the user specifies the input files and chosen options with a single command. In case of wrong commands, VaRank will print a description of the options and defaults settings.

VaRank reads SNV and indel variant descriptions from a VCF file (variant call format Danecek et al., 2011), the reference format for genomic variants. The VCF file can be gzip compressed. The program can take any combination of a single VCF file with multiple patients and/or multiple VCF files with single patient’s data as input (Fig. 2).

Each variant is checked for consistency (genotype, depth of coverage, variant call especially for indels) and VaRank prints warnings if appropriate (i.e., patients redundancy based on the sample identifiers in the VCFs), information from each analysis step (including status and running time), as well as statistics on the submitted data such as total number of variants and number of patients.

The output from VaRank is presented in an easily-accessible tab-separated file that can be opened in any spreadsheet program. Output files report one variant per line. Nevertheless, two types of rankings are provided: one presenting each variant independently, ordered from most likely pathogenic to least likely pathogenic (files with “byVar” suffix), and another ranking genes from most likely causative to least likely causative (files with “byGene” suffix). For the latter, each gene is scored along two criteria: (i) based on its homozygous most pathogenic variant, or (ii) based on its first two heterozygous most pathogenic variants. In order to make sure that no variants are overlooked by the user (by only displaying the most pathogenic variants) all other gene variants are also reported. The “byVar” file is more appropriate for analyzing patient’s data under the dominant or pseudominant hypothesis, while the “byGene” file is more appropriate for recessive diseases (especially for compound heterozygous cases).

Each of these output files is also available in two versions: one contains all submitted variants (“AllVariants” files) while the other one is prefiltered (“filteredVariants” files). The default filters remove variants: (i) called with a read depth <=10, (ii) with a supporting read count <=10, (iii) with a ratio of supporting reads <=15%, (iv) with a validated status annotation in the dbSNP (Sherry et al., 2001) database (based on at least two supporting evidences) that are not pathogenic (based on the ClinicalSignificance field), and (v) with an allele frequency >1% (extracted from the dbSNP database or the Exome Variant Server). All these parameters can be modified by the user.

Finally, a short report of counts (homozygous, heterozygous and total counts) for each of the variant categories (5′ and 3′ UTR, upstream, downstream, frameshift, in-frame, nonsense, splice site, start loss, stop loss, missense, synonymous, intronic, not annotated) is generated for each sample and for the whole submitted dataset.

Variant annotation

The annotation of variants is performed by the annotation engine (Fig. 2). It is composed of several parts including the main annotation software, which can be either Alamut Batch or SnpEff, the barcode annotation (see Fig. 3 and corresponding section) and optionally PolyPhen-2 predictions. As an example, Alamut Batch collects among others (Table 2): the gene symbol, the OMIM ID, the transcript ID (i.e., RefSeq), the protein IDs (RefSeq and UniProt), the HGVS nomenclature (genomic, cDNA and proteic), information from public variation databases such as dbSNP and the Exome Variant Server (EVS, http://evs.gs.washington.edu/EVS/) and predicted effects at both nucleotide and protein levels. When available, known mutations are highlighted by extracting either reported SNVs/indels flagged as “probably-pathogenic”/“pathogenic” in the field “Clinical significance” introduced since dbSNP134 or using the HGMD database (Stenson et al., 2014).

The choice of transcript is a critical task that can lead to misannotation (McCarthy et al., 2014). To avoid underestimating variant effects, they are annotated on all transcripts available (i.e., one variant can be either intronic or exonic depending on the isoform) and the most pathogenic effect is retained. Given that VaRank is compatible with 2 annotation software there are small differences. Using Alamut Batch, each variant is scored for each transcript. By default we report the longest transcript for each gene except if any variation is more pathogenic in another transcript. In the case of SnpEff, the annotations are already sorted from the most to the least pathogenic.

In order to further enrich the annotation for each variant and each gene, VaRank can integrate (using the option -extann) external annotations provided by the user as a tab separated file. One could for instance associate private expression data or transmission mode for each gene of interest.

Annotating and analyzing several individuals together can be very computationally effective since most of the variants are common polymorphisms. As an example, looking at 180 patients sequenced for 217 genes, a total of 204,625 variations could be identified where only 9,378 were non redundant. In this case, the separate annotation of each patient’s variant set would have required ∼20× times the computational cost of the combined analysis. Although the total number of non-redundant variants does not plateau, each new sample adds only a very limited number of new variants to the analysis (Fig. 4).

Variant ranking

The observed variants (SNVs/indels) can be characterized at different levels (DNA, RNA and protein levels) that VaRank aims at summarizing into a single score. This score is then used to rank variants based on their predicted pathogenicity and thus accelerates identification of relevant ones by biologists. The aim of this score is not to provide yet another score to assess the pathogenicity of each variant but a rationale to present the most relevant variants according to the biologist common use and interpretation rules. Thus, the relative weights of each score components were determined experimentally to best separate categories. VaRank uses the variation type (i.e., substitution, deletion, insertion, duplication) and the coding effect to score. The VaRank scoring is categorized from the most likely to the less likely pathogenic state as follows (score in parenthesis): known mutation (110), nonsense (100), frameshift (100), start loss (80), stop loss (80), missense (50), in-frame (40) and synonymous coding (10) (Table 1).

In addition when using Alamut Batch, each variant is assessed for any potential effect on the nearest splice site. Following the guidelines from Houdayer et al. (2012) and our own tests (data not shown), we selected three assessment programs: MaxEntScan (Yeo & Burge, 2004), NNSplice (Reese et al., 1997) and Splice Site Finder (based on Shapiro & Senapathy, 1987). A variant is considered to affect splicing when at least two out of the three programs indicate a significant score change between the wild type and the mutated sequences (respectively −10%, −5% and −15%). A VaRank score is then attributed depending on the following three splice site categories (score into parenthesis): (i) essential splice site: two first intronic bases up- or downstream the exon (90), (ii) intron-exon boundary: donor site from −3 to +6, acceptor site from −12 to +2 (70) and (iii) deep intronic changes (25). It should be noted that a SNV affecting the first or last base of an exon can either have a coding effect as a missense or an effect on splicing. VaRank reports the most pathogenic of these two possible effects.

If appropriate the variant score is further adjusted using additional information as nucleotide-level conservation (phastCons Siepel et al., 2005) and protein-level pathogenicity predictions (SIFT and PolyPhen-2) (Adzhubei et al., 2010; Kumar, Henikoff & Ng, 2009) that are used to compute an adjustment score (0 or +5) to be added to the relevant category (Table 1).

Barcode

Comparing variations of several individuals, related and/or unrelated, affected or not, has proven to be a very effective strategy for distinguishing polymorphisms from variants causing or increasing the likelihood of disease (Ng et al., 2010; Ng et al., 2009). In order to take advantage of this, VaRank introduces a barcode that allows a quick overview of the presence/absence status of each variant within all samples and their zygosity status (“0” representing homozygous wild type, “1” heterozygous and “2” homozygous for the variant) (Fig. 3A). In Fig. 3B, three variants are reported for one specific patient out of a cohort of 32 samples analyzed together. For example, the third variation c.601G>A in TTC21B is heterozygous for this patient. In light of the presented barcode, one can immediately notice that the variant is also present in 28 other samples from the cohort, of which in total 12 are homozygous and 17 heterozygous.

In order to allow inheritance analysis, a second barcode (the family barcode) representing only user selected samples can be defined. As an example in trio exome sequencing, we have represented two typical pedigrees (Figs. 3B and 3C), one consanguineous family on the left and one sporadic case on the right. In the case of consanguinity, homozygous mutations are often the cause of the disease in the family. This could be highlighted by selecting the “121” barcode indicating homozygous variants (“2”) in the proband inherited from heterozygous parents (“1”). In the sporadic case, several hypotheses could be tested including a de novo variant which could be highlighted using the barcode “010.”

Together with the barcode, simple counts on the individuals (homozygous, heterozygous and total allelic counts) are also available and can easily be used to further filter variants. Indeed, in rare diseases such as the Bardet-Biedl syndrome (BBS, OMIM# 209900), mutations are often private (i.e., one mutation found only in one family) (Muller et al., 2010) meaning that their frequency in the population is very low. Counts can be used to estimate the frequency of a known variant in the user cohort and add significant value to variants not yet reported in any public variant database but for which a frequency can be estimated based on the user’s cohort. As an example, looking at 2,888 non redundant SNVs observed in 107 patients with moderate to severe intellectual disability, 979 did not have any frequency information in the dbSNP and EVS databases. Such information could be directly retrieved from the VaRank output.

The observed frequency of variants in public databases but also in private cohorts can be a powerful filtering strategy. Using the same data (Fig. 3B), variant c.7911dup in ALMS1 is present only once in the cohort of patients at the homozygous state and is very likely the disease causing mutation in this patient.

Bardet-Biedl syndrome (BBS) dataset

The Bardet-Biedl syndrome (BBS; OMIM# 209900) is a pleiotropic recessive disorder, part of the ciliopathies, characterized by extensive genetic heterogeneity counting to date 19 genes (Aldahmesh et al., 2014; Scheidecker et al., 2014). We applied targeted high-throughput sequencing for 30 ciliopathy related genes to 52 patients with clinical features compatible with BBS (Redin et al., 2012). VaRank was used to annotate and rank the variants identified in those patients. Thirty-two cases could be solved by this approach leading to frameshift, missense and splice site mutations all validated by Sanger sequencing (we excluded Copy Number Variations). Sequencing data from the 32 positive samples have been reanalyzed using Alamut Batch version 1.1.11 and PolyPhen-2 v2.2.2 installed on our local servers. A total of 784 non redundant variants have been annotated resulting on average into 167 private variants per sample. We extracted the validated mutations and highlighted the ranking position in the output files (Table 3). In 30/32 samples, mutations were ranked first, while in the remaining ones they were present in the top five among ∼170 variants per patient. This result clearly shows the effectiveness of the ranking in such situation.

As mentioned in the original paper, one variant is always ranked in first position and represents a false positive described in BBS2 as a third allele mutation according to the triallelic hypothesis (NM_031885.3:c.209G>A, rs4784677). It is flagged as pathogenic in dbSNP, but it is too frequent to be a fully penetrant mutation according to the observed frequency in the Exome Variant Server (EVS) (0.77%). Interestingly, using the barcode this variant is very easily filtered out since it is present in almost all patients in our cohort (31/32 samples).

Intellectual disability (ID) dataset

Intellectual disability is a common neurodevelopmental disorder (∼2% of children and adolescents) (Ellison, Rosenfeld & Shaffer, 2013) that can be initiated by either environmental, genetic or multifactorial causes. The monogenic forms are very heterogeneous genetically, with several hundred genes identified so far. We present here the results of 203 patients affected with moderate to severe intellectual disability.

A first cohort of 107 patients was already analyzed by VaRank (Redin et al., 2014) and led to the identification of 25 causal mutations (Table 4A). A second cohort of 96 additional patients is reported here for which 12 causative mutations could be identified (Table 4B). The identified mutations were all validated by Sanger sequencing (we excluded pathogenic CNV). High-throughput sequencing data of protein-coding exons of 217 genes (first cohort) and 275 genes (second cohort) (either on the X-chromosome or associated to autosomal dominant or recessive forms of ID) were collected for the 25 positive known patients and the 96 additional ones. A total of 8,388 non redundant variants have been reannotated by VaRank using Alamut Batch version 1.1.11 and PolyPhen-2 v2.2.2, resulting on average into 1,129 (±189) private variants per sample. Each sample has been reanalyzed and a summary of the ranking results using the default filtration strategy for this dataset defined by the biologists (frequency in dbSNP or EVS >1% and presence in <3 samples) is available in Table 4. Among the 37 patients with causative mutations identified almost all were ranked first or within the top five among 1,129 variants on average demonstrating the usefulness of this strategy.

As an example, in a boy with autism spectrum disorder, attention deficit and autoaggressive behavior, one mutation (c.797_798delinsTT, p.C266F) has been observed in the MAOA gene. A total of 688 variants have been annotated by VaRank and the mutation described for this patient could be ranked without filtering at position 8. When the usual filters are applied (frequency in dbSNP or EVS >1%), the mutation is ranked at the first position in the “filteredVariants.rankingByVar.txt” file. This result was the first mutation report in the MAOA gene since 20 years (Piton et al., 2014).

Whole Exome Sequencing (WES) dataset

WES of a single patient with a clear BBS phenotype from a consanguineous Italian family (Scheidecker et al., 2014) revealed for the first time mutations in the BBIP1 gene counting since as the 18th BBS gene. A total of 50,569 non redundant variants have been annotated by VaRank using Alamut Batch version 1.1.11 and PolyPhen-2 v2.2.2.

The nonsense mutation c.173T>G (p.Leu58*) in BBIP1 has been ranked at the 50th position/50569 in the “allVariants.rankingByVar” file. Applying the default filtration strategy (i.e., frequency in dbSNP or EVS >1% and sequence quality filters) changed the ranking to the 6th position/4,908 (“filteredVariants.rankingByVar” file). This exome sequencing dataset was analyzed together with 29 other exomes from unrelated patients and pathologies for a BBS patient for which we could further filter variations using for example the barcode. Given that the most frequent known disease causing mutation in BBS is the c.1169T>G (p.M390R) mutation in BBS1 (Mykytyn et al., 2002), found in EVS at the frequency of 26/12,694 (0.2%) at the heterozygous state, we used the barcode statistics to further reduce the total number of variations. Being very tolerant, we kept variants present less than four times at the heterozygous state or two times at the homozygous state out of 29 samples. The final ranking placed the mutation at the forth position (and first homozygous) out of 3,493 remaining variants.

The integration of these three datasets of increasing complexity (BBS with 30 genes consolidated on 188 patients, ID with more than 200 genes in 121 patients and the WES data consolidated for 35 samples) highlights major directions for further developments. Considering the distribution of the non-redundant variation by functional category, one can first observe, despite a different gene composition, that the distribution of the categories is similar among the datasets (Fig. 5). The vast majority of the identified and annotated variations are either intronic or synonymous. Those categories contain known variations in dbSNP (∼85% have an rs#) and are either polymorphisms or rare variant. Thus they are often considered as non-pathogenic. Nevertheless, some of these could potentially affect the correct splicing of the closest gene. Although efforts are being taken to enhance the predictions on splicing especially for the consensus sequences (Houdayer et al., 2012; Jian, Boerwinkle & Liu, 2014), little is done in more distant regions including enhancers sequences, branch point or even promoter regions. The 5’ and 3’ UTRs are also sources of variations that are often overlooked but contain important functional signals such as miRNA binding sites and polyadenylation signal (Chatterjee & Pal, 2009). Missenses are one of the major variations sources and also one of the more difficult to interpret. There is a high number of predictions methods and tools available aiming at predicting the pathogenicity of this category of variant but there is still improvement to be made (Flanagan, Patch & Ellard, 2010; Hicks et al., 2011; Tavtigian et al., 2008). Recent approaches aimed at combining knowledge based information such as structural information for missense variants (Luu et al., 2012) or focusing only on orthologous sequences (Wong & Zhang, 2014). Other interesting annotation sources might be included as the assessment of non-frameshift indels (Bermejo-Das-Neves et al., 2014; Hu & Ng, 2012) and the consideration of CNVs to improve the decision-making. Those issues will be amplified using whole genome sequencing (WGS).