Received: 23 January 2019 | Revised: 17 April 2019 | Accepted: 24 April 2019 DOI: 10.1002/humu.23772 INFORMATICS eDiVA—Classification and prioritization of pathogenic variants for clinical diagnostics Mattia Bosio1,2,3 | Oliver Drechsel4 | Rubayte Rahman5 | Francesc Muyas1,2 | Raquel Rabionet 1,2,6 | Daniela Bezdan 1,2 | Laura Domenech Salgado 1,2 | Hyun Hor 7 | Jean‐Jacques Schott 8,9 | Francina Munell 10 | Roger Colobran10 | Alfons Macaya 10 | Xavier Estivill 11,12 | Stephan Ossowski 1,2,13 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Barcelona, Spain 2 Universitat Pompeu Fabra (UPF), Barcelona, Spain 3 Barcelona Supercomputing Center (BSC), Barcelona, Spain 4 Bioinformatics Unit (MF1), Department for Methods Development and Research Infrastructure, Robert Koch Institute, Berlin, Germany 5 NKI Netherlands Cancer Institute, The Netherlands 6 Institut de Recerca Sant Joan de Déu, University of Barcelona, Barcelona, Spain 7 Department of Neurology, University Hospital Zurich, Zurich, Switzerland 8 L'Institut du Thorax, INSERM, CNRS, Univ Nantes, Nantes, France 9 Service de Cardiologie, L'institut du thorax, CHU Nantes, Nantes, France 10 Vall d'Hebron Institut de Recerca (VHIR), Barcelona, Spain 11 Sidra Medicine, Doha, Qatar 12 Women's Health Dexeus, Barcelona, Spain 13 Institute of Medical Genetics and Applied Genomics, University of Tübingen, Tübingen, Germany Correspondence Stephan Ossowski, Institute of Medical Abstract Genetics and Applied Genomics, University of Mendelian diseases have shown to be an and efficient model for connecting Tübingen, Calwerstr. 7, 72076 Tübingen, Germany. genotypes to phenotypes and for elucidating the function of genes. Whole‐exome Email: Stephan.Ossowski@med.uni-tuebingen. sequencing (WES) accelerated the study of rare Mendelian diseases in families, de allowing for directly pinpointing rare causal mutations in genic regions without the Funding information need for linkage analysis. However, the low diagnostic rates of 20–30% reported Spanish Ministry of Economy and Competitiveness; “la Caixa” Foundation; for multiple WES disease studies point to the need for improved variant ‘Centro de Excelencia Severo Ochoa 2013‐ pathogenicity classification and causal variant prioritization methods. Here, we 2017; European Union ‐ H2020 research and innovation programme, Grant/Award Number: present the exome Disease Variant Analysis (eDiVA; http://ediva.crg.eu), an 635290 ‐ PanCanRisk; CERCA Programme ‐ automated computational framework for identification of causal genetic variants Generalitat de Catalunya (coding/splicing single‐nucleotide variants and small insertions and deletions) for rare diseases using WES of families or parent–child trios. eDiVA combines next‐ generation sequencing data analysis, comprehensive functional annotation, and causal variant prioritization optimized for familial genetic disease studies. eDiVA features a machine learning‐based variant pathogenicity predictor combining [Correction added on 11 June 2019, after first online publication: Co‐author's name changed from Hyun‐Gyu Hor to Hyun Hor] --------------------------------------------------------------------------------------------------------------------------- This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2019 The Authors. Human Mutation Published by Wiley Periodicals, Inc. Human Mutation. 2019;40:865–878. wileyonlinelibrary.com/journal/humu | 865 866 | BOSIO ET AL. various genomic and evolutionary signatures. Clinical information, such as disease phenotype or mode of inheritance, is incorporated to improve the precision of the prioritization algorithm. Benchmarking against state‐of‐the‐art competitors demonstrates that eDiVA consistently performed as a good or better than existing approach in terms of detection rate and precision. Moreover, we applied eDiVA to several familial disease cases to demonstrate its clinical applicability. KEYWORDS disease variant prioritization, machine learning, NGS diagnostics, rare genetic disease, whole‐ exome sequencing 1 | INTRODUCTION although higher success rates have been reported for specific disease or inheritance types (Sawyer et al., 2016) and for studies using parent– Rare genetic diseases are classical models for studying gene function and child trios (Yang et al., 2013). While some of the unsolved cases might linking genotypes to disease phenotypes. Although each of these diseases be explained by intergenic or intronic regulatory variation or only affects a small number of patients, the sum of people affected by one unidentified structural variants, the low detection rate also indicates of the more than 7,000 rare diseases exceeds 30 million individuals in the the need for development of better prioritization strategies for coding US alone (Cutting, 2014). Whole‐exome sequencing (WES), and more variants and robust classifiers comprehensively integrating the available recently whole‐genome sequencing (WGS), are routinely applied to amount of prior omics and the knowledge of the disease. identify variants causing rare Mendelian diseases in studies of families or Many computational algorithms have been developed to assess parent–child trios (Choi et al., 2009; Ng et al., 2010; Louis‐Dit‐Picard pathogenicity of genetic variants. Tools such as SIFT (Kumar, Henikoff, et al., 2012; Rabbani, Mahdieh, Hosomichi, Nakaoka, & Inoue, 2012). & Ng, 2009), CADD (Kircher et al., 2014), PolyPhen‐2 (Adzhubei et al., Usually, each exome sequencing experiment yields tens of 2010), or Eigen (Ionita‐Laza, McCallum, Xu, & Buxbaum, 2016) are thousands of genetic variants in coding and splicing regions that require commonly used in clinical practice to help variant interpretation. They thorough functional annotation and filtering to allow identification of derive a functional impact score based on amino acid or nucleotide the causal variant. Several tools have been published performing variant conservation, and biochemical properties of the amino acid changes as annotation, including Annovar, VEP, or SNPeff, which augment the features. While some algorithms additionally categorize variants into sequencing information with a comprehensive set of current omics, various categories such as “neutral,” “benign,” “deleterious,” “damaging,” population genomics, and clinical knowledge (Cingolani et al., 2012; “probably‐damaging,” or “pathogenic” (e.g., SIFT, Condel, PolyPhen‐2, McLaren et al., 2016; Wang, Li, & Hakonarson, 2010). These tools utilize and Mutation Assessor), scores of other methods need to be interpreted a large selection of available databases containing gene annotations, by using (often arbitrary) cutoffs for pathogenicity (e.g., CADD). These various genomic features, variant allele frequencies in different predicted pathogenicity labels are an integral part of the American populations, functional impact prediction, and evolutionary conservation College of Medical Genetics and Genomics standards and guidelines for (Bao et al., 2014). Other methods, such as eXtasy (Sifrim et al., 2013), the interpretation of sequence variants (Richards et al., 2015). Methods PhenoDB (Sobreira, Schiettecatte, Boehm, Valle, & Hamosh, 2015), combining multiple classifiers, such as MetaLR, have been shown to Phen–Gen (Javed, Agrawal, & Ng, 2014), VarSifter (Teer, Green, produce better results than single classifiers (Dong et al., 2015). Mullikin, & Biesecker, 2012), KGGseq (M.‐X. Li, Gui, Kwan, Bao, & Recently, specialized ensemble learning methods for estimating Sham, 2012), and SPRING (Wu, Li, & Jiang, 2014), focus on prioritization pathogenicity of rare variants have been published: Mendelian Clinically of potentially causal variants using both functional annotation and Applicable Pathogenicity (M‐CAP; Jagadeesh et al., 2016), using clinical information. These tools systematically filter, evaluate, and gradient‐boosting trees on pathogenicity scores and conservation prioritize thousands of variants, taking into account knowledge found in features, and Revel (Ioannidis et al., 2016), using an RF to integrate genome annotation databases (Rhead et al., 2010), disease gene several pathogenicity predictors. repositories (OMIM, Online Mendelian Inheritance in Man; Landrum To combine an intuitive user interface with comprehensive variant et al., 2014), and patient pedigree information, as well as phenotype prediction, annotation, pathogenicity classification, and causal variant descriptions and disease definitions provided for example, as Human prioritization we developed eDiVA (exome Disease Variant Analysis), Phenotype Ontology (HPO) terms (Köhler et al., 2014). Finally, methods http://www.ediva.crg.eu. The eDiVA pipeline is composed of four main such as Endeavour (Tranchevent et al., 2008) and GeneDistiller (Seelow components: (a) eDiVA‐Predict, where sequencing results are processed et al., 2008) prioritize disease genes, not individual variants, by to predict the presence of genomic variants; (b) eDiVA‐Annotate, that integrating diverse genomic data sources. enriches variants via a domain‐knowledge database; (c) eDiVA‐Score, Detection rates of causal variants using WES have been reported to which estimates variant pathogenicity using a random forest model; and be as low as 20–30% of cases (H. Lee et al., 2014; Yang et al., 2013), (d) eDiVA‐Prioritize, in which variants from small groups of related BOSIO ET AL. | 867 samples (i.e., families or parent–child trios) are analyzed jointly. eDiVA (https://www.broadinstitute.org/gatk/guide/best‐practices as of June returns a shortlist of candidate variants compatible with the selected 2017) to extract genetic variants from raw reads. Reads in fastq disease inheritance model and the pedigree information. Using the format are aligned using bwa‐mem (H. Li, 2013), alignments are post‐ pathogenicity probability computed by eDiVA‐Score, variants are ranked processed using samtools (H. Li, 2011), GATK (McKenna et al., 2010), such that better candidates appear on top of the result list. eDiVA has Picard (Picard Tools—By Broad Institute), and custom quality filters been developed with specific emphasis on usability, automation, and (details provided in Supporting Information Material). Finally, VCF reproducibility of results and is available as a web service with a graphical files are generated using GATK HaplotypeCaller. Subsequent user interface (see Supporting Information Material), or as an open‐ regenotyping of all positions harboring a single‐nucleotide variant source repository with Docker containers. eDiVA can be run using the (SNV) or small insertions and deletions in at least one family member NextFlow (Di Tommaso et al., 2017) pipeline management system to yields a complete matrix of variants for the whole sample set (family) ensure its compatibility with most standalone or cloud‐computing in multisample VCF format. Due to the computational resources platforms as well as to guarantee reproducibility on any system. required for read alignment and variant calling, eDiVA‐Predict is eDiVA has been optimized for two common clinical diagnostics currently not enabled on the eDiVA web server, but can be used with scenarios, parent–child trios comprised of healthy parents and one the standalone version of eDiVA on a local or remote computing affected child (tested for recessive, compound heterozygous, and X‐ infrastructure (e.g., Amazon Cloud). Alternatively, variant prediction linked inheritance or dominant de novo variants) and families with can be performed using any tool able to produce one multisample multiple affected relatives (additionally tested for dominant inheri- VCF file reporting genotype quality and coverage information for all tance). We demonstrate that eDiVA outperforms competing approaches variable positions (e.g., GATK, McKenna et al., 2010; freebayes, in a semisynthetic benchmark study introducing thousands of known Garrison & Marth, 2012). disease variants from ClinVar (Landrum et al., 2014) or HGMD (Stenson et al., 2017) into real WES data from the 1000 Genomes Project CEPH 2.3 | eDiVA‐Annotate: Functional variant parent–offspring trio of European ancestry (NA12878, NA12891, and annotation NA12892). We, furthermore, report summary statistics on eDiVA and Phen–Gen results for 35 unreported disease cases, composed of 15 Using the eDiVA‐Annotate module each variant is individually cases of spinocerebellar ataxia, 16 cases of primary immunodeficiency, linked with public information sources to integrate multiple and four cases of congenital myasthenia. knowledge domains, and to provide a comprehensive annotation profile. First, ANNOVAR (Wang et al., 2010) is applied to relate each variant to its corresponding gene (choosing among UCSC, 2 | MATERIALS AND METHODS Ensembl, or Refseq gene annotations), and to its functional consequence at the protein level. Next, functional, population 2.1 | eDiVA pipeline genomics, and evolutionary data relevant for variant prioritization eDiVA consists of a Python pipeline combined with an SQL Database are added to each variant. To this end we created a MySQL back‐end composed of four components: variant prediction, variant database, eDiVA‐DB, containing all relevant positional information annotation, pathogenicity estimation, and variant prioritization obtained from UCSC table browser (Rhead et al., 2010) and other (Figure S1). The main functionality of eDiVA is to process next‐ sources. Each variant is annotated with population allele fre- generation sequencing (NGS) data for small sets of samples (e.g., quency information from the dbSNP (Sherry et al., 2001), families or parent–child trios) and to output a shortlist of potentially discovEHR (Dewey et al., 2016), 1000 Genomes Project causal variants for the diagnosed disease. eDiVA is available as an open‐ (1000GP; The 2015 Genomes Project Consortium, 2015), Exome source repository, https://github.com/mbosio85/ediva, with a Docker Variant Server (Exome Variant ), and GnomAD exomes (Lek et al., container composition wrapped within a NextFlow (Di Tommaso et al., 2016) databases. The latter three databases also provide informa- 2017) interface to guarantee exact reproducibility on the most common tion on specific populations (e.g., Caucasian, Asian, African computing platforms (including several cloud platforms) and as a freely American, etc.), which can be selected for improved causal variant accessible web server: http://www.eDiVA.crg.eu. The modular nature of prioritization. Information on evolutionary conservation is incor- eDiVA allows for easy integration of specific parts, for example, the porated from PhyloP (Rhead et al., 2010), and PhastCons (Hubisz, eDiVA‐Score module for pathogenicity estimation, in other pipelines or Pollard, & Siepel, 2011), including conservation scores for tools. Comprehensive examples for the use of eDiVA and example input primates, mammals, and vertebrates. Precalculated scores for files are included in the repository and on the website. functional impact of variants have been integrated based on the algorithms SIFT (Kumar, Henikoff, & Ng, 2009), PolyPhen‐2 (Adzhubei et al., 2010), Mutation Assessor (Reva, Antipin, & 2.2 | eDiVA‐Predict: WES or WGS processing and Sander, 2011), Condel (González‐Pérez & López‐Bigas, 2011), variant calling Eigen (Ionita‐Laza et al., 2016), and CADD (Kircher et al., 2014). The eDiVA‐Predict module performs sample‐wise variant calling Furthermore, eDiVA‐DB includes information on genomic features according to the recent GATK (McKenna et al., 2010) best practices like segmental duplications and simple sequence repeats provided 868 | BOSIO ET AL. by UCSC table browser (Rhead et al., 2010). Finally, eDiVA‐DB (eDiVA‐Score, see above); (b) removal of all variants that do not provides clinical data from ClinVar (Landrum et al., 2014) and segregate according to the selected inheritance mode; (c) OMIM related to each variant and affected gene. filtering based on functional and population genomic features; eDiVA‐Annotate uses multisample VCF files and returns a file and (d) prioritization based on user defined clinical phenotypes with annotated variants in comma‐separated value format. This step (as HPO IDs). Filtering based on segregation requires the user to can be performed on the eDiVA web server. submit a simple pedigree file defining the relationship between samples and their disease state (i.e., affected or unaffected), and to choose the most likely inheritance pattern for the disease (or 2.4 | eDiVA‐Score: Estimating variant to run all modes). eDiVA‐Prioritize can process variants following pathogenicity five types of inheritance patterns: (a) dominant de novo, (b) eDiVA's prioritization algorithm relies on accurate estimation of autosomal dominant inherited, (c) autosomal recessive homo- pathogenicity for each variant. We therefore developed eDiVA‐ zygous, (d) autosomal recessive compound heterozygous, or (e) Score, a machine learning classifier, which assigns a pathogenicity X‐linked. probability to each variant based on its annotation characteristics Optionally, eDiVA removes variants that are improbable of obtained from eDiVA‐Annotate. eDiVA‐Score is built by training a being damaging, are likely false‐positive calls or do not have random forest (RF) model using the R “randomForest” package sufficient read coverage in all family members to reliably with 1000 binary classification trees (Breiman, 2001; Hastie, estimate segregation patterns. By default, eDiVA applies a Tibshirani, & Friedman, 2009) and five‐fold cross validation. Eleven lenient filter setting defined in Table S1. Finally, eDiVA allows features were selected to train the RF model: (a) the maximum the user to specify a list of HPO terms (Köhler et al., 2014) minor allele frequency (MAF) of 1000Genomes and GnomAD relevant for the disease as an additional source of information to databases; (b) four conservation measures (conservation in prioritize variants in genes. eDiVA highlights all variants in genes primates and mammals using the PhastCons (Hubisz et al., 2011) related to the submitted phenotypic traits using a custom and PhyloP (Pollard, Hubisz, Rosenbloom, & Siepel, 2010); (c) four algorithm to estimate the HPO‐gene association (detailed in the functional impact predictors: Condel (González‐Pérez & López‐ Supporting Information material). Bigas, 2011), Phred‐scaled CADD score (Kircher et al., 2014), Eigen (Ionita‐Laza et al., 2016), and Mutation Assessor (Reva et al., 2.6 | Performance evaluation using semisynthetic 2011); (d) the likelihood to be in a segmental duplication, which cases correlates with false‐positive variant calls (Ho, Tsai, Chen, & Lin, 2011); and (e) an in‐house estimator of systematic sequencing To assess the performance of eDiVA and several competing methods, errors called ABB‐score (Muyas et al., 2019). Note that Condel, we implemented a semisynthetic benchmark based on real WES data Eigen and CADD are combination scores integrating several from a trio in which we spiked‐in known pathologic variants from the features also included in eDiVA‐score, namely evolutionary ClinVar database (Landrum et al., 2014). We chose a publicly available conservation (PhastCons and PhyloP in mammals and primates) CEPH trio sequenced within the framework of the 1000 Genomes and Mutation Assessor scores. The RF model has been trained Project composed of samples with European ancestry NA12878 using 15,000 random pathogenic and likely pathogenic variants (daughter), NA12891 and NA12892 (parents), downloadable from from the ClinVar database (Landrum et al., 2014) as positive cases. https://public_docs.crg.es/sossowski/MicrobeGenomes/human/eDiVA/ We then built a control set composed of 15,000 nonpathogenic insilico_simulation_data/, and we called variants and generated a variants from ClinVar, and 100,000 random variants from multisample VCF file using eDiVA‐Predict. For the purpose of this GnomAD (Lek et al., 2016) not contained in ClinVar. The vast benchmark study, all 138,705 variants found in the original trio are majority of variants in both positive and negative training set are considered true negatives, that is, variants not associated with the rare (allele frequency [AF], <1%; Figure S2a,b), thus circumventing disease. that AF dominates the classification model. Following the neutral Next, we embedded known disease variants in the trio following theory of molecular evolution (Kimura, 1983) missing data is segregation patterns expected for Mendelian diseases. This positive generated using expected values for nonpathogenic (neutral) set, containing variants associated with the diseases, consists of all variants (Figure 1). The only exception is AF, as missing data in variants from ClinVar (Landrum et al., 2014) database labeled as the context of AF means that the SNV is novel, that is, has AF of “pathogenic” or “likely pathogenic”, having an OMIM reference in the zero. Variants used for training of the RF have been excluded in all database and that had not been used for training of eDiVA‐Score. For benchmarking tests performed in this study. each pathogenic variant, we extracted: chromosome, position, reference and alternative nucleotides, dbSNP identifier, gene name, inheritance mode of the associated disease (where available, 2.5 | eDiVA‐Prioritize: Causal variant prioritization randomly assigned otherwise), and HPO terms for the disease. Causal variant prioritization consists of four steps, (a) ranking by Variants without HPO annotation have been excluded from the estimated probability of variants to cause a phenotypic change benchmark set. BOSIO ET AL. | 869 (a) (b) AF CADD Eigen Phred 1.00 40 AF 60 0.75 30 CADD 40 0.50 20 0.25 20 Eigen_ 10 0.00 0 0 Not pathogenic Pathogenic Not pathogenic Pathogenic Not pathogenic Pathogenic Mammals PhyloP Mammals PhyloP Mutation Assessor Condel 3 1.00 Mutation 4 Assessor 0.75 0 2 Primates 0.50 PhastCons −3 0 0.25 Condel −6 −2 0.00 −4 Mammals Not pathogenic Pathogenic Not pathogenic Pathogenic Not pathogenic Pathogenic PhastCons Primates PhastCons Mammals PhastCons Primates PhyloP 1.00 1.00 0 Primates PhyloP 0.75 0.75 −2 ABB_score 0.50 0.50 −4 0.25 0.25 Segmental Duplication −6 0.00 0.00 0 2000 4000 Importance Not pathogenic Pathogenic Not pathogenic Pathogenic Not pathogenic Pathogenic F I G U R E 1 eDiVA‐Score random forest model. (a) Estimated importance of features used in the model (extracted with varImp command). (b) Distribution of values for the top‐9 features used in the model, comparing ClinVar pathogenic against ClinVar nonpathogenic variants. AF: allele frequency; eDiVA: exome Disease Variant Analysis We have simulated three inheritance patterns: autosomal the alternative allele (i.e., the fraction of reads showing the recessive homozygous, autosomal recessive compound heterozygous, alternative allele) introduced in the original VCF file has been and dominant de novo, as these are the most likely patterns found in obtained using a beta distribution and a binomial distribution for parent–child trio based rare‐disease diagnostics. To create realistic homozygous and heterozygous variants, respectively. A total of 6,811 disease genotypes, each pathogenic variant was introduced into the disease‐associated variants from ClinVar not previously used in the exomes of the daughter and the parents, if applicable according to training of eDiVA‐Score were used for benchmarking: 3,353 the inheritance mode. The read distribution of reference and recessive homozygous, 2,592 dominant de novo, and 866 recessive alternative reads was simulated depending on the inheritance mode compound heterozygous disease‐causing variants (see Table S2 for and the original coverage data. The variant allele frequency (VAF) of additional information on simulated genotypes). 870 | BOSIO ET AL. 2.7 | Benchmarking of variant pathogenicity recessive compound heterozygous”, “autosomal recessive homozy- estimation methods gous”, or “autosomal dominant new mutation”. We locally installed Phen–Gen and launched it with the corresponding setups: (a) We evaluated the ability of eDiVA‐Score and six competing methods, “Recessive”, “allow_de_novo = 0” for recessive and compound inheri- namely CADD, Eigen, DANN, Revel, M‐CAP, and MetaLR (Dong et al., tance, and (b) “Dominant”, “allow_de_novo = 1” for the dominant de 2015, Ioannidis et al., 2016, Ionita‐Laza et al., 2016, Jagadeesh et al., novo inheritance model. We locally installed Exomiser and analyzed 2016, Kircher et al., 2014, Quang, Chen, & Xie, 2015), to prioritize all trio cases using PhenIX prioritization mode (details in Supporting pathogenic over benign variants. We generated a receiver operating Information Material). We tested eDiVA in two configurations, (a) characteristic (ROC) curve for each tool and benchmark set and without phenotype description, and (b) using HPO IDs describing the measured performance by area under the curve (AUC). disease phenotype for disease‐specific prioritization of candidate First, we evaluated the performance of each method on the variants. ClinVar test set (containing only variants not used for model To benchmark the ability of eDiVA, Exomiser, PhenoDB, and training), using variants labeled “pathogenic” as true positives (TP) Phen–Gen to distinguish disease‐causing from benign variants we and variants labeled “benign” as true negatives (TN; Figure 2a–c). compared three quality metrics, (a) recall (i.e., did the causal variant Second, we benchmarked using variants from the HGMD and appear in the output list or not), (b) average number of false positives GnomAD databases (not used in model training or present in across all benchmarked cases as a proxy for precision, and (c) ranks ClinVar) as TP and TN, respectively (Figure 2d–f). Third, we of causal variants reported for each mode of inheritance using violin measured the performance of all methods on HGMD data only, plots (Figure 3a–c). To compare ranks, variants reported by eDiVA using the categories for damaging and likely damaging mutation (DM are sorted by eDiVA‐Score, Phen–Gen results are sorted by DCOD‐ and DM?) as TP and any other HGMD category as false positives (FP) score (“Probability of deleteriousness based on genic predictor”), and (Figure 3g–i). Functional impact values for the benchmarked methods Exomiser results are sorted by “Exomiser Gene Combined Score”. have been obtained from the respective publications. CADD, DANN, Results of PhenoDB are presented in the default order (chromosome and eDiVA provide damage estimates for all positions of the genome, and position), as no prioritization score is provided. and Eigen for close to 70% of all positions, whereas Revel, M‐CAP, and MetaLR are trained specifically for rare (AF, <1%) or known variants and are only available for a subset of ClinVar, HGMD, and 3 | RESULTS GnomAD. We, therefore, performed three separate performance tests for each of the three benchmark sets, applying the following 3.1 | eDiVA: A platform for pathogenicity criteria (a) using only variants having Revel and M‐CAP scores estimation and causal variant prioritization available (ClinVar: 3,887 TP and 10,494 TN; HGMD/GnomAD: eDiVA is a disease variant prioritization tool optimized for NGS‐ 63,712 TP and 100,000 TN; HGMD: 63,712 TP and 1,892 TN); (b) based genetic disease diagnostics in families and parent–child trios. It random subset of all variants, assigning a default value of 0 to missing is composed of four components: eDiVA‐Predict handles read values (ClinVar: 19,888 TP and 16,694 TN; HGMD/GnomAD: 96,569 alignment and variant prediction, eDiVA‐Annotate performs func- TP and 100,000 TN; HGMD: 96,569 TP and 7376 TN); and (c) using tional annotation of variants, eDiVA‐Score estimates the probability only rare variants (AF, ≤0.01) from the previous pool of variants of variants to be pathogenic, and eDiVA‐Prioritize filters and ranks (ClinVar: 16,531 TP and 15,531; HGMD/GnomAD: 90,004 TP and variants according to various quality criteria, proper segregation, and 97,828 TN; HGMD: 96,004 TP and 2,817 TN). likelihood to cause phenotypic changes. eDiVA is available as Furthermore, we studied five variant sets provided by Grimm standalone software at https://github.com/mbosio85/ediva, and as et al. (2015), forming a collection of data sets for benchmarking a web service providing access to functional annotation, pathogeni- pathogenicity classifiers published in independent studies. Finally, we city classification and causal variant prioritization modules (www. combined these five sets to form a combined benchmark (see ediva.crg.eu). The eDiVA web service facilitates analysis of families or Supporting Information Materials for details). parent–child trios in a few clicks, requiring only a VCF file, and optionally a set of HPO IDs describing the disease phenotype. eDiVA returns a shortlist of candidate variants and genes, ranked by 2.8 | Benchmarking of disease variant pathogenicity score (together with gene relatedness to the specified prioritization methods HPO IDs if available), and including all annotation features in comma‐ We compared eDiVA with three commonly used tools for variant separated value (.csv) and Microsoft Excel (.xlsx) format. annotation and prioritization: Exomiser (Robinson et al., 2014), PhenoDB (Sobreira et al., 2015), and Phen–Gen (Javed et al., 2014) 3.2 | Benchmarking eDiVA and competing methods on a set of 6,811 semisynthetic parent–child trios (see above). PhenoDB was executed from the https://phenodb.org/ website using To comprehensively evaluate eDiVA's performance in finding standard parameters (a) AF, <0.01; (b) including variants which are disease‐causing variants, and to compare it to previously published present in dbSNP, and (c) analysis type chosen among “autosomal tools, we performed a benchmark in two categories. First, we BOSIO ET AL. | 871 F I G U R E 2 Benchmarking of the pathogenicity classifiers eDiVA‐Score, CADD, Eigen, Revel, and M‐CAP using ROC for (a) set of 10,494 ClinVar pathogenic variants (TP) and 3,887 ClinVar “benign” variants (TN); (b) set of 16,694 ClinVar pathogenic variants (TP) and 19,888 ClinVar “benign” variants (TN), setting missing values to benign, (c) subset of rare variants (AF, <1% from set c); (d) set of 63,712 variants from HGMD (TP) and 100,000 from GnomAD (TN) for which values from all tools are available; (e) set of 96,569 variants from HGMD (TP) and 100,000 from GnomAD (TN), setting missing values to benign; (f) subset of rare variants (AF, <1% from set e); (g) set of 63,712 HGMD variants (“DM” and “DM?”) as TP, and 1,892 HGMD variants (other categories) as TN for which values from all tools are available; (h) set of 96,569 variants from HGMD (“DM” and “DM?”) as TP, and 7,376 HGMD (other categories) as TN, setting missing values to benign; and (i) subset of rare variants (AF, <1% from set h). AF: allele frequency; eDiVA: exome Disease Variant Analysis; M‐CAP: Mendelian clinically applicable pathogenicity; ROC: receiver operating characteristic; TN: true negative; TP: true positive 872 | BOSIO ET AL. F I G U R E 3 Receiver operating characteristic curves comparing pathogenicity classifiers on five independent data sets (and the combined set) composed of pathogenic and neutral variants. Revel, M‐CAP, and eDiVA show a similarly strong performance, with the exception of the PredictSNP and Varibench sets, on which Revel and M‐CAP outperform eDiVA‐Score. eDiVA: exome Disease Variant Analysis; M‐CAP: Mendelian clinically applicable pathogenicity; SNP: single‐nucleotide polymorphism evaluated the ability of eDiVA‐Score to distinguish disease‐causing specialized method for estimating functional impact or pathogenicity from benign variants compared to four publicly available methods for (Figure S5). estimating deleteriousness. Second, we benchmarked the perfor- We benchmarked the ability of eDiVA‐Score, CADD, DANN, mance of eDiVA, PhenoDB, and Phen–Gen on identification of causal Eigen, Revel, M‐CAP, and MetaLR to predict the deleteriousness of variants using semisynthetic parent–child trios analyzed by WES, variants and to distinguish pathogenic from benign variants in nine optionally allowing for the use of clinical phenotype descriptions for setups (Section 2). We first compared the performance on classifying causal variant prioritization. pathogenic and benign variants from ClinVar (Figure 2a), on distinguishing disease variants from HGMD (Stenson et al., 2017) from 100,000 random variants from GnomAD (Figure 2d), for which 3.3 | Benchmarking of eDiVA‐Score and other scores are available for all methods. Note that Revel and M‐CAP variant pathogenicity classifiers have been trained on a subset of the HGMD variants (e.g., using class We developed eDiVA‐Score, a machine learning‐based method for “DM” as positive training set), giving them an advantage due to estimating variant pathogenicity (deleteriousness) independent of potential overfitting in any of the following benchmark tests using any prior clinical information (see Section 2). Feature‐selection HGMD variants (for an in‐depth discussion of the interplay between identified population allele frequency, functional impact, and con- overfitting and circularity in training and benchmarking data (Grimm servation in placental mammals as the most important features et al., 2015). Using ROC analysis, we found that eDiVA‐Score (Figure 1a). The correlation matrix for all features is shown in distinguishes disease‐causing and benign variants with high sensitiv- Figure S3. Features selected for inclusion in the RF show distinct ity and recall in both benchmark sets (AUC of 0.95 and 0.90), distributions for pathogenic variants compared to benign variants in considerably better than CADD (AUC of 0.91 and 0.74), DANN (AUC ClinVar (Figure 1b), random coding variants reported in GnomAD of 0.89 and 0.82), Eigen (AUC of 0.87 and 0.77), Revel (AUC of 0.91 (Figure S4b). All integrated conservation scores (PhyloP and and 0.89), M‐CAP (AUC of 0.84 and 0.90), and MetaLR (AUC of 0.88 PhastCons scores for vertebrates, mammals and primates) classify and 0.87). Of note, eDiVA‐Score showed better precision‐recall pathogenic variants better than random, but perform worse than any curves than competing methods (Figure S6). BOSIO ET AL. | 873 Disease variant prioritization tools depend on pathogenicity more general the functional impact) of variants, whereas eDiVA‐ values for any position of the exome, since de novo mutations can Score, Revel, and M‐CAP have been trained to identify pathogenic occur randomly and novel ultra‐rare variants are still being variants, partly explaining the divergent performance levels across discovered. Therefore, we next benchmarked the methods on the different benchmark sets. Moreover, eDiVA‐score, MetaLR, and random variants chosen from the complete ClinVar and HGMD/ M‐CAP use CADD as one of many features, explaining the better GnomAD benchmark sets, whereas setting missing data to benign performance of the derived scores. (Section 2). As expected, the recall of Revel, M‐CAP, and MetaLR In summary, our benchmark results demonstrate the good decreased substantially due to missing information, whereas the performance of eDiVA‐Score as pathogenicity classifier, comparable other methods performed slightly better than in the previous tests to and often better than state‐of‐the‐art methods available to date. (Figure 2b,e). Finally, we tested how the methods perform on Furthermore, eDiVA‐Score outperforms other general‐purpose classification of rare variants (AF, <0.01), otherwise following the methods not restricted by variant AF (i.e. CADD, DANN, and Eigen), same criteria for selection of the test sets as in the previous while showing competitive results when compared with specialized benchmark (Figure 2c,f). Again, eDiVA‐Score shows the best tools such as MetaLR, M‐CAP, and Revel, which are only available for performance of all methods. known (rare) SNVs. We wondered if the use of random GnomAD variants as TN (nonpathogenic) set might bias the results of the HGMD/GnomAD 3.4 | Causal variant prioritization in parent–child benchmark due to for example, overfitting onto the allele frequency trios feature. Therefore, we next measured the performance of all methods on HGMD data only, using the categories for highly likely We benchmarked the performance of eDiVA and three widely used pathogenic (“DM” and “DM?”) as TP set and less likely pathogenic tools, PhenoDB, Phen–Gen, and Exomiser, on identification of causal (any other HGMD category) as TN set (Section 2). We performed the variants for rare Mendelian diseases in parent–child trios. To this same three tests as discussed above for the ClinVar and HGMD/ end, we simulated three scenarios typically encountered in parent– GnomAD benchmark sets. On the subset of variants for which scores child trio diagnostics, (a) autosomal dominant de novo, (b) autosomal are available for all methods (Figure 2g) eDiVA's performance (AUC recessive homozygous, and (c) autosomal recessive compound 0.77) was found to be slightly lower than MetaLR's (0.80), Revel's heterozygous Mendelian inheritance modes. In total, we simulated (AUC 0.82), and M‐CAP's (AUC 0.85), but substantially better than 6,811 semisynthetic parent–child trios by integrating reported the performance of the other general‐purpose methods CADD (AUC pathogenic variants from ClinVar into real WES data of a trio 0.67) and Eigen (AUC of 0.70). However, eDiVA still outperformed all obtained from 1000GP (see Section 2 and Table S2). other methods on the complete HGMD variant set (missing scores Figure 3a shows violin plots with the rank distribution of causal set to benign), as well as on the rare variant set (Figure 2h,i). variants in the output lists of 6,811 analyzed trios. The optimal result Finally, we compared the performance of all methods on a is a skewed distribution close to zero, meaning that the causal variant benchmark set compiled by Grimm et al. (2015), consisting of is reported as first or very close to the top of the list in the majority mutually exclusive subsets of the previously published benchmark of cases. Here, comparison with PhenoDB is not meaningful, as sets Varibench, HumVar, ExoVar, predictSNP, and SwissVar (see PhenoDB (unlike Phen–Gen, Exomiser, and eDiVA), offers no ranking Supporting Information Material for details). These popular bench- based on pathogenicity scores (but sorts by chromosome and mark data sets differ in the way they define pathogenic and neutral position). Compared to Exomiser and Phen–Gen, eDiVA's ranking variants, for example, the maximum AF for pathogenic variants can method shows the best performance for recessive homozygous differ dramatically, allowing us to benchmark diverse challenges. inheritance, eDiVA and Exomiser show best performance for Furthermore, Grimm et al. filtered these benchmark sets to minimize dominant de novo inheritance, and all tools show similarly good overlap between them, reducing the likelihood that tools are performance for compound heterozygous inheritance. eDiVA con- benchmarked on variants they have been trained on and hence sistently reported causal recessive homozygous variants and com- reducing the impact of overfitting on the benchmark results (Grimm pound heterozygous variants within the top five candidates et al., 2015). We found that none of the methods consistently (median = 1), and dominant de novo variants within the top 25 of performs better than other methods, but that eDiVA‐Score, M‐CAP, reported candidates (median = 4; Figure 4a). Considering that the Revel, and MetaLR show comparably high performance, except on CEPH trio has been sequenced as part of the 1000GP we finally PredictSNP and Varibench, for which MetaLR, Revel, and M‐CAP tested if the use of 1000GP allele frequency information for filtering show a better performance than eDiVA‐Score. PredictSNP incorpo- biases the performance estimates of eDiVA. However, we found no rates HGMD variants in the positive and negative control sets; see difference when not using the 1000GP AF database (Figure S7). Tables 2 and 3 of Grimm et al. (2015), likely to be giving a strong Nonetheless, we cannot exclude the possibility that eDiVA (or the advantage to Revel and M‐CAP, which have been trained on HGMD. other methods) show reduced performance in understudied popula- CADD, DANN, and Eigen performed significantly worse than the tions. other three methods on all benchmark sets. Note that CADD, DANN, The use of HPO IDs for prioritization further reduced FPs Eigen, and MetaLR have been trained to predict deleteriousness (or reported by eDiVA (label eDiVA‐HPO in Figure 4c). Overall, we 874 | BOSIO ET AL. (a) (b) (c) F I G U R E 4 Benchmark of the causal variant prioritization tools eDiVA, Exomiser, Phen–Gen, and PhenoDB. (a) Violin plots showing the rank of disease‐causing variants within the reported candidate lists for the three tested inheritance types: “recessive homozygous”, “compound heterozygous”, and “dominant de novo”; (b) Recall values for 6,811 semisynthetic trio cases, representing the fraction of identified causal variants (i.e., “solved cases”). (c) Average number of false positives reported per case as a proxy for precision. eDiVA has been tested in two configurations, with HPO‐based gene prioritization (eDiVA_HPO) and with the default configuration not using HPO terms (eDiVA). Adding HPO filtering reduces false positives at the cost of a slightly reduced Recall. HPO: Human Phenotype Ontology observed a two‐fold reduction in FPs across all inheritance modes simulation with incomplete phenotyping show a small reduction in tested. However, filtering by in silico gene panels also resulted in a causal gene ranking efficiency (Figure S8). However, even imperfect reduction in recall (Figure 4b,c). Finally, we observed improved phenotypic information improved the performance as compared to prioritization ranks under all inheritance types, with the strongest complete exclusion of such information. impact seen for de novo variants (Figure 4a). In summary, the benchmark analyses show that eDiVA achieves We also investigated the impact of incomplete or imperfect highly competitive causal variant prioritization performance with phenotyping on eDiVA's performance by altering the HPO annota- respect to ranking, precision and recall, while requiring no fine tuning tion of genetic variants imported from ClinVar (see Supporting of parameters by the user for specific inheritance types. When Information Methods). Benchmarking results on the semisynthetic disease‐specific HPO term descriptors are available, eDiVA's BOSIO ET AL. | 875 precision is further enhanced to the point at which complete overlooked due to unrealistic simulated benchmark tests not automation of causal variant identification is feasible for recessive reflecting the level of noise found in real family or trio NGS data. homozygous and compound heterozygous segregation. We have addressed these problems by developing eDiVA, a pipeline that combines multisample variant calling of family data, QC and filtering, extensive functional annotation, machine learning‐ 3.5 | eDiVA results on clinical cases based classification of deleterious variants, and prioritization of eDiVA has successfully been used in published case studies on mitral causal variants optimized for various clinical scenarios. Furthermore, valve prolapse (Durst et al., 2015), cystic fibrosis (Ramos et al., 2014), we developed a highly realistic benchmark test combining real WES phenylketonuria (Trujillano et al., 2014), arthrogryposis (Wambach data of a parent–child trio with thousands of pathogenic ClinVar et al., 2017), and Opitz‐C (Urreizti et al., 2017), among others, variants to generate 6,811 semisynthetic disease trios. Using these identifying both known, as well as novel rare‐disease genes. We cases, we have demonstrated that eDiVA's pathogenicity estimator recently assessed the performance of eDiVA for the diagnosis of rare (eDiVA‐Score) as well as eDiVA's prioritization algorithm perform congenital genetic diseases using WES of 35 parent–child trios, favorably compared to existing state‐of‐the‐art methods. eDiVA has including 15 cases of congenital ataxia, four cases of congenital been able to find disease‐causing variants with higher recall, fewer myasthenia, and 16 cases of primary immunodeficiency. Here we false positives and better ranking than competing tools in three report general statistics on the number of candidate genes per case, benchmarked modes of inheritance. Finally, we evaluated the use of while case reports and novel candidate genes will be published in phenotypic descriptors for optimizing the prioritization process. separate papers. Across the 35 studied trios, eDiVA on average We found that adding HPO ID‐based prioritization introduces a reported a median of 11, 3, and 10 candidate genes per trio for trade‐off between recall and the number of false positives in the recessive homozygous, recessive compound heterozygous, and output list. Despite the marginal reduction in recall, focusing on dominant de novo inheritance mode, respectively, using default known disease genes is often the preferred choice for diagnostic parameters. In comparison, Phen–Gen reported a median of 36 and purposes. Our knowledge of genetic factors playing a role in disease 52 candidate genes for recessive (including compound heterozygous) is constantly growing, reflected in a rapid increase of genotype– and dominant (including de novo) inheritance mode, respectively. phenotype relations stored in various databases. Hence, it would be Histograms of reported candidate gene numbers for eDiVA and beneficial to reanalyze WES data sets once in a while (e.g., every Phen–Gen are shown in Figures S9 and S10. eDiVA found causal 6–12 months) to benefit from new knowledge and to facilitate variants in known genes for the respective disease in none cases, and identification of previously unknown/unreported causal variants. variants in genes associated with closely related disease phenotypes Moreover, combined reanalysis of the growing cohorts of WES data in seven cases. Screening of Phen–Gen results did not reveal stored in many institutes would allow to identify matching causal additional candidates missed by eDiVA. The function of a novel genes across multiple families or cases. However, most analysis disease variant for congenital ataxia has been described in pipelines require substantial hands‐on time and long candidate‐ Bahamonde et al. (2015), and reports for other candidate genes are variant lists have to be screened by experts, making regular in preparation. reanalysis of data sets impractical. eDiVA has been developed with a specific emphasis on high reproducibility of results and complete automation of the analysis using artificial intelligence‐based methods. 4 | D IS C U S S IO N Machine learning classifiers are used to perform candidate ranking and prioritization, reducing hands‐on time of clinical experts to a Despite the massive increase in sequencing capacity and the minimum. Integration with NextFlow, moreover, guarantees repro- availability of highly optimized analysis tools, multiple large‐scale ducibility of results at later time points and on most computing rare‐disease studies reported that in only 20–30% of cases a causal platforms. Therefore, eDiVA is a dedicated solution for regular variant can be identified using WES. Several reasons might explain reanalysis of large disease cohorts or collections of diagnostic cases. the inability of WES analysis to identify causal variants in a majority Additional steps can be taken to improve the identification of of cases, including for example, the inability to identify regulatory disease‐related variants from WES data. The availability of custom allele variants (Claussnitzer et al., 2015), our limited knowledge of the frequency databases with geographical specificity would help to reduce function of noncoding RNAs, generation of new exon donor or the number of false‐positive genotype–phenotype associations due to acceptor sites by intronic variants (Y. Lee et al., 2012), small copy population specific variants. To this end, institutes and hospitals with number variations (Krumm et al., 2012), incomplete penetrance, and access to large cohorts of sequenced exomes may use in‐house data to unknown function of coding genes, among others. However, we argue filter population specific variants, an approach we have pursued our self that the potential of WES has not been exhausted and that causal by collecting thousands of Iberian cases in an aggregated allele coding variants are often missed due to inappropriate correction of frequency database (http://geevs.crg.eu/, unpublished). Identification noise in the data, insufficient use of clinical (phenotypic) data, or of extended homozygosity regions could in addition help to diagnose reporting of long unranked candidate lists, requiring tedious screen- causal homozygous variants in consanguineous cases. Moreover, the ing by clinicians. We further claim that these shortcomings are often integration of structural and copy number variants (SVs and CNVs) has 876 | BOSIO ET AL. been shown to increase recall rates substantially (Gambin et al., 2017). R E F E R E N CE S Despite their frequent involvement in rare diseases (McCarroll & Adzhubei, I. A., Schmidt, S., Peshkin, L., Ramensky, V. E., Gerasimova, A., Altshuler, 2007), CNVs are often disregarded in WES analyses, and are Bork, P., … Sunyaev, S. R. (2010). A method and server for predicting rarely processed in combination with point mutations. Prioritization damaging missense mutations. Nature Methods, 7, 248–249. algorithms will have to be adapted to consider compound heterozygotes Bahamonde, M. I., Serra, S. A., Drechsel, O., Rahman, R., Marcé‐Grau, A., composed of a point mutation in one and a CNV in the other allele. CNV Prieto, M., … Fernández‐Fernández, J. M. (2015). A single amino acid deletion (ΔF1502) in the S6 segment of CaV2.1 domain III associated analysis is currently being integrated in eDiVA and will be available in with congenital ataxia increases channel activity and promotes Ca2+ the near future. Influx. PLOS One, 10, 10. Better use of phenotypic descriptors has the potential to improve Bao, R., Huang, L., Andrade, J., Tan, W., Kibbe, W. A., Jiang, H., & Feng, G. both precision and recall of causal variant prioritization methods. We (2014). Review of current methods, applications, and data manage- ment for the bioinformatics analysis of whole exome sequencing. observed that HPO ID‐based prioritization dramatically improved Cancer Informatics, 13, 67–82. the precision of eDiVA. However, incomplete maps of known Breiman, L. (2001). Random forests. Machine Learning, 45, 5–32. genotype–phenotype (or gene–phenotype) relations in public data- Choi, M., Scholl, U. I., Ji, W., Liu, T., Tikhonova, I. R., Zumbo, P., … Farhi, A. bases led to a mild reduction in recall. Robinson et al. (2014) (2009). Genetic diagnosis by whole exome capture and massively parallel DNA sequencing. Proceedings of the National Academy of proposed a method to overcome this limitation, tapping into the Sciences of the United States of America, 106, 19096–19101. genotype–phenotype associations from mouse data to solve causal Cingolani, P., Platts, A., Wang, L., Coon, M., Nguyen, T., Wang, L., … Ruden, D. variant identification for corresponding human phenotypes. Other M. (2012). A program for annotating and predicting the effects of single methods based on image analysis, for example, Hadj‐Rabia et al. nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso‐2; iso‐3. Fly (Austin), 6(2), 80–92. (2017) or face2gene (http://suite.face2gene.com/), have also shown Claussnitzer, M., Dankel, S. N., Kim, K.‐H., Quon, G., Meuleman, W., promising results for diagnosis of patients with visible phenotypic Haugen, C., … Liu, J. (2015). FTO obesity variant circuitry and features. Finally, an important step in the evaluation of newly adipocyte browning in humans. New England Journal of Medicine, 373, discovered genotype–phenotype associations is the identification of 895–907. Cutting, G. R. (2014). Annotating DNA variants is the next major goal for additional cases with a similar phenotype and mutations in the same human genetics. American Journal of Human Genetics, 94, 5–10. gene. Several approaches for gene matching have been published, for Dewey, F. E., Murray, M. F., Overton, J. D., Habegger, L., Leader, J. B., example, GeneMatcher (Sobreira et al., 2015), which have been Fetterolf, S. N., … Pendergrass, S. A. (2016). Distribution and clinical connected via the Matchmaker Exchange platform. Integration of impact of functional variants in 50,726 whole‐exome sequences from approaches using image analysis, cross‐species phenotype–genotype the DiscovEHR study. Science, 354, aaf6814. Di Tommaso, P., Chatzou, M., Floden, E. W., Barja, P. P., Palumbo, E., & correlation, and gene matching has the potential to further improve Notredame, C. (2017). Nextflow enables reproducible computational AI‐based variant prioritization methods such that they can rival the workflows. Nature Biotechnology, 35, 316–319. diagnostic precision of clinical experts in the future. Durst, R., Sauls, K., Peal, D. S., deVlaming, A., Toomer, K., Leyne, M., … Jett, In summary, we have shown that eDiVA is a step towards full C. (2015). Mutations in DCHS1 cause mitral valve prolapse. Nature, 525, 109–113. automation of causal variant identification in family and parent–child Dong, C., Wei, P., Jian, X., Gibbs, R., Boerwinkle, E., Wang, K., & Liu, X. trio data using machine learning‐based approaches. eDiVA can be (2015). Comparison and integration of deleteriousness prediction used as a support tool for clinicians to find disease‐causing variants, methods for nonsynonymous SNVs in whole exome sequencing or as a fully automated solution for periodic reanalysis of large WES studies. Human Molecular Genetics, 24, 2125–2137. Exome Variant Server (n.d.). http://evs.gs.washington.edu/EVS/ (or WGS) cohorts. eDiVA is able to identify known causal disease Gambin, T., Akdemir, Z. C., Yuan, B., Gu, S., Chiang, T., Carvalho, C. M. B., … variants with high precision and recall, and facilitates identification of Bayram, Y. (2017). Homozygous and hemizygous CNV detection from novel disease variants with minimal hands‐on time. exome sequencing data in a Mendelian disease cohort. Nucleic Acids Research, 45, 1633–1648. Garrison, E., & Marth, G. (2012). Haplotype‐based variant detection from A C K N O W L E D GM E N T S short‐read sequencing. ArXiv12073907 Q‐Bio, https://arxiv.org/abs/ 1207.3907 We acknowledge support of the Spanish Ministry of Economy and González‐Pérez, A., & López‐Bigas, N. (2011). Improving the assessment of Competitiveness, Centro de Excelencia Severo Ochoa 2013–2017, the outcome of nonsynonymous SNVs with a Consensus Deleterious- and the CERCA Programme/Generalitat de Catalunya. This project ness Score, Condel. American Journal of Human Genetics, 88, 440–449. has received funding from the “la Caixa” Foundation, the CRG Grimm, D. G., Azencott, C.‐A., Aicheler, F., Gieraths, U., MacArthur, D. G., Samocha, K. E., … Borgwardt, K. M. (2015). The evaluation of tools emergent translational research award and the European Union's used to predict the impact of missense variants is hindered by two H2020 Research and Innovation Programme under the grant types of circularity. Human Mutation, 36, 513–523. agreement No 635290 (PanCanRisk). Hadj‐Rabia, S., Schneider, H., Navarro, E., Klein, O., Kirby, N., Huttner, K., … Grange, D. K. (2017). Automatic recognition of the XLHED phenotype from facial images. American Journal of Medical Genetics. OR CID Part A, 173, 2408–2414. Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical Mattia Bosio http://orcid.org/0000-0001-8989-4220 learning: data mining, Inference, and Prediction, New York: Springer. Francesc Muyas http://orcid.org/0000-0002-7857-0623 https://www.springer.com/kr/book/9780387848570 BOSIO ET AL. | 877 Ho, M.‐R., Tsai, K.‐W., Chen, C., & Lin, W. (2011). dbDNV: A resource of McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., duplicated gene nucleotide variants in human genome. Nucleic Acids Kernytsky, A., … DePristo, M. A. (2010). The Genome Analysis Toolkit: Research, 39, D920–D925. A MapReduce framework for analyzing next‐generation DNA Hubisz, M. J., Pollard, K. S., & Siepel, A. (2011). PHAST and RPHAST: sequencing data. Genome Research, 20, 1297–1303. Phylogenetic analysis with space/time models. Briefings in Bioinfor- McLaren, W., Gil, L., Hunt, S. E., Riat, H. S., Ritchie, G. R. S., Thormann, A., … matics, 12, 41–51. Cunningham, F. (2016). The Ensembl variant effect predictor. Genome Ionita‐Laza, I., McCallum, K., Xu, B., & Buxbaum, J. D. (2016). A spectral Biology, 17, 122. approach integrating functional genomic annotations for coding and Muyas, F., Bosio, M., Puig, A., Susak, H., Domènech, L., Escaramis, G., … noncoding variants. Nature Genetics, 48, 214–220. Ossowski, S. (2019). Allele balance bias identifies systematic genotyping Ioannidis, N. M., Rothstein, J. H., Pejaver, V., Middha, S., McDonnell, S. K., errors and false disease associations. Human Mutation, 40, 115–126. Baheti, S., … Teerlink, C. C. (2016). REVEL: An ensemble method for Ng, S. B., Buckingham, K. J., Lee, C., Bigham, A. W., Tabor, H. K., Dent, K. predicting the pathogenicity of rare missense variants. American M., … Bamshad, M. J. (2010). Exome sequencing identifies the cause of Journal of Human Genetics, 99, 877–885. a mendelian disorder. Nature Genetics, 42, 30–35. Jagadeesh, K. A., Wenger, A. M., Berger, M. J., Guturu, H., Stenson, P. D., Pollard, K. S., Hubisz, M. J., Rosenbloom, K. R., & Siepel, A. (2010). Cooper, D. N., … Bejerano, G. (2016). M‐CAP eliminates a majority of Detection of nonneutral substitution rates on mammalian phyloge- variants of uncertain significance in clinical exomes at high sensitivity. nies. Genome Research, 20, 110–121. Nature Genetics, 48(12), 1581–1586. Quang, D., Chen, Y., & Xie, X. (2015). DANN: A deep learning approach for Javed, A., Agrawal, S., & Ng, P. C. (2014). Phen‐Gen: Combining phenotype annotating the pathogenicity of genetic variants. Bioinforma Oxf Engl, and genotype to analyze rare disorders. Nature Methods, 11, 935–937. 31, 761–763. Kimura, M. (1983). The Neutral Theory of Molecular Evolution. Cmbridge, Rabbani, B., Mahdieh, N., Hosomichi, K., Nakaoka, H., & Inoue, I. (2012). Next‐ UK: Cambridge University Press. generation sequencing: Impact of exome sequencing in characterizing Kircher, M., Witten, D. M., Jain, P., O'Roak, B. J., Cooper, G. M., & Shendure, J. Mendelian disorders. Journal of Human Genetics, 57, 621–632. (2014). A general framework for estimating the relative pathogenicity of Ramos, M. D., Trujillano, D., Olivar, R., Sotillo, F., Ossowski, S., Manzanares, J., human genetic variants. Nature Genetics, 46, 310–315. … Vazquez, C. (2014). Extensive sequence analysis of CFTR, SCNN1A, Krumm, N., Sudmant, P. H., Ko, A., O'Roak, B. J., Malig, M., Coe, B. P., … SCNN1B, SCNN1G and SERPINA1 suggests an oligogenic basis for Eichler, E. E. (2012). Copy number variation detection and genotyping cystic fibrosis‐like phenotypes. Clinical Genetics, 86, 91–95. from exome sequence data. Genome Research, 22, 1525–1532. Reva, B., Antipin, Y., & Sander, C. (2011). Predicting the functional impact Kumar, P., Henikoff, S., & Ng, P. C. (2009). Predicting the effects of coding of protein mutations: Application to cancer genomics. Nucleic Acids non‐synonymous variants on protein function using the SIFT Research, 39, e118–e118. algorithm. Nature Protocols, 4, 1073–1081. Rhead, B., Karolchik, D., Kuhn, R. M., Hinrichs, A. S., Zweig, A. S., Fujita, P. Köhler, S., Doelken, S. C., Mungall, C. J., Bauer, S., Firth, H. V., Bailleul‐ A., … Pheasant, M. (2010). The UCSC genome browser database: Forestier, I., … Eppig, J. T. (2014). The Human Phenotype Ontology Update 2010. Nucleic Acids Research, 38, D613–D619. project: Linking molecular biology and disease through phenotype Richards, S., Aziz, N., Bale, S., Bick, D., Das, S., Gastier‐Foster, J., … Rehm, data. Nucleic Acids Research, 42, D966–D974. H. L. (2015). Standards and guidelines for the interpretation of Landrum, M. J., Lee, J. M., Riley, G. R., Jang, W., Rubinstein, W. S., Church, sequence variants: A joint consensus recommendation of the D. M., & Maglott, D. R. (2014). ClinVar: Public archive of relationships American College of Medical Genetics and Genomics and the among sequence variation and human phenotype. Nucleic Acids Association for Molecular Pathology. Genet Med Off J Am Coll Med Research, 42, D980–D985. Genet, 17, 405–424. Lee, H., Deignan, J. L., Dorrani, N., Strom, S. P., Kantarci, S., Quintero‐ Robinson, P. N., Köhler, S., Oellrich, A., Wang, K., Mungall, C. J., Lewis, S. E., Rivera, F., … Fogel, B. L. (2014). Clinical exome sequencing for genetic … Haendel, M. (2014). Improved exome prioritization of disease genes identification of rare Mendelian disorders. Journal of the American through cross‐species phenotype comparison. Genome Research, 24, Medical Association, 312, 1880–1887. 340–348. Lee, Y., Gamazon, E. R., Rebman, E., Lee, Y., Lee, S., Dolan, M. E., … Lussier, Sawyer, S. L., Hartley, T., Dyment, D. A., Beaulieu, C. L., Schwartzentruber, Y. A. (2012). Variants affecting exon skipping contribute to complex J., Smith, A., … Warman Chardon, J. (2016). Utility of whole‐exome traits. PLOS Genetics, 8, e1002998. sequencing for those near the end of the diagnostic odyssey: Time to Lek, M., Karczewski, K. J., Minikel, E. V., Samocha, K. E., Banks, E., Fennell, address gaps in care. Clinical Genetics, 89, 275–284. T., … Birnbaum, D. P. (2016). Analysis of protein‐coding genetic Seelow, D., Schwarz, J. M., & Schuelke, M. (2008). GeneDistiller—Distilling variation in 60,706 humans. Nature, 536, 285–291. Candidate Genes from Linkage Intervals. PLoS ONE, 3(12), e3874. Li, H. (2011). A statistical framework for SNP calling, mutation discovery, https://doi.org/10.1371/journal.pone.0003874 association mapping and population genetical parameter estimation Sherry, S. T., Ward, M. H., Kholodov, M., Baker, J., Phan, L., Smigielski, E. from sequencing data. Bioinforma Oxf Engl, 27, 2987–2993. M., & Sirotkin, K. (2001). dbSNP: The NCBI database of genetic Li, H. (2013). Aligning sequence reads, clone sequences and assembly variation. Nucleic Acids Research, 29, 308–311. contigs with BWA‐MEM. ArXiv13033997 Q‐Bio Sifrim, A., Popovic, D., Tranchevent, L.‐C., Ardeshirdavani, A., Sakai, R., Li, M.‐X., Gui, H.‐S., Kwan, J. S. H., Bao, S.‐Y., & Sham, P. C. (2012). A Konings, P., … Moreau, Y. (2013). eXtasy: Variant prioritization by comprehensive framework for prioritizing variants in exome sequencing genomic data fusion. Nature Methods, 10, 1083–1084. studies of Mendelian diseases. Nucleic Acids Research, 40, e53–e53. Sobreira, N., Schiettecatte, F., Boehm, C., Valle, D., & Hamosh, A. (2015). Louis‐Dit‐Picard, H., Barc, J., Trujillano, D., Miserey‐Lenkei, S., Bouatia‐ New tools for mendelian disease gene identification: PhenoDB Naji, N., Pylypenko, O., … Jeunemaitre, X. (2012). KLHL3 mutations Variant Analysis Module; and GeneMatcher, a web‐based tool for cause familial hyperkalemic hypertension by impairing ion transport linking investigators with an interest in the same gene. Human in the distal nephron. Nature Genetics, 44(4), 456–460. https://doi.org/ Mutation, 36, 425–431. 10.1038/ng.2218 Teer, J. K., Green, E. D., Mullikin, J. C., & Biesecker, L. G. (2012). VarSifter: McCarroll, S. A., & Altshuler, D. M. (2007). Copy‐number variation and Visualizing and analyzing exome‐scale sequence variation data on a association studies of human disease. Nature Genetics, 39, S37–S42. desktop computer. Bioinformatics (Oxford, England), 28, 599–600. 878 | BOSIO ET AL. Stenson, P. D., Mort, M., Ball, E. V., Evans, K., Hayden, M., Heywood, S., … Wang, K., Li, M., & Hakonarson, H. (2010). ANNOVAR: Functional Cooper, D. N. (2017). The Human Gene Mutation Database: Towards annotation of genetic variants from high‐throughput sequencing data. a comprehensive repository of inherited mutation data for medical Nucleic Acids Research, 38, e164–e164. research, genetic diagnosis and next‐generation sequencing studies. Wu, J., Li, Y., & Jiang, R. (2014). Integrating multiple genomic data to Human Genetics, 136, 665–677. predict disease‐causing nonsynonymous single nucleotide variants in The 1000 Genomes Project Consortium (2015). A global reference for exome sequencing studies. PLOS Genetics, 10, e1004237. human genetic variation. Nature, 526, 68–74. Yang, Y., Muzny, D. M., Reid, J. G., Bainbridge, M. N., Willis, A., Ward, P. A., Tranchevent, L.‐C., Barriot, R., Yu, S., Van Vooren, S., Van Loo, P., Coessens, … Person, R. (2013). Clinical whole‐exome sequencing for the B., … Moreau, Y. (2008). ENDEAVOUR update: a web resource for gene diagnosis of mendelian disorders. New England Journal of Medicine, prioritization in multiple species. Nucleic Acids Research, W377–W384. 369, 1502–1511. 36 (Web Server). https://doi.org/10.1093/nar/gkn325 Trujillano, D., Perez, B., González, J., Tornador, C., Navarrete, R., Escaramis, G., … Estivill, X. (2014). Accurate molecular diagnosis of SU PP ORT IN G IN FOR M ATI O N phenylketonuria and tetrahydrobiopterin‐deficient hyperphenylalani- nemias using high‐throughput targeted sequencing. European Journal Additional supporting information may be found online in the of Human Genetics EJHG, 22, 528–534. Supporting Information section at the end of the article. Urreizti, R., Cueto‐Gonzalez, A. M., Franco‐Valls, H., Mort‐Farre, S., Roca‐ Ayats, N., Ponomarenko, J., … Hecht, J. (2017). A de novo nonsense mutation in MAGEL2 in a patient initially diagnosed as Opitz‐C: Similarities between Schaaf‐Yang and Opitz‐C syndromes. Scientific How to cite this article: Bosio M, Drechsel O, Rahman R, et al. Reports, 7, srep44138. eDiVA—Classification and prioritization of pathogenic Wambach, J. A., Stettner, G. M., Haack, T. B., Writzl, K., Škofljanec, A., Maver, variants for clinical diagnostics. Human Mutation. 2019;40: A., … Baldridge, D. (2017). Survival among children with “Lethal” congenital contracture syndrome 11 caused by novel mutations in the 865–878. https://doi.org/10.1002/humu.23772 gliomedin gene (GLDN). Human Mutation, 38, 1477–1484.