VarAFT: a variant annotation and filtration system for human next generation sequencing data.

With the rapidly developing high-throughput sequencing technologies known as next generation sequencing or NGS, our approach to gene hunting and diagnosis has drastically changed. In <10 years, these technologies have moved from gene panel to whole genome sequencing and from an exclusively research context to clinical practice. Today, the limit is not the sequencing of one, many or all genes but rather the data analysis. Consequently, the challenge is to rapidly and efficiently identify disease-causing mutations within millions of variants. To do so, we developed the VarAFT software to annotate and pinpoint human disease-causing mutations through access to multiple layers of information. VarAFT was designed both for research and clinical contexts and is accessible to all scientists, regardless of bioinformatics training. Data from multiple samples may be combined to address all Mendelian inheritance modes, cancers or population genetics. Optimized filtration parameters can be stored and re-applied to large datasets. In addition to classical annotations from dbNSFP, VarAFT contains unique features at the disease (OMIM), phenotypic (HPO), gene (Gene Ontology, pathways) and variation levels (predictions from UMD-Predictor and Human Splicing Finder) that can be combined to optimally select candidate pathogenic mutations. VarAFT is freely available at: http://varaft.eu.

INTRODUCTION

Massively parallel sequencing, also called NGS (next generation sequencing), led to a genetic revolution with the ability to sequence any human genome in a few hours. Nevertheless, despite the thousands of exomes and genomes that have been studied (Genome Aggregation Database (gnomAD) (1)), we still have only a limited vision of the human genome variability especially in the context of rare human genetic disease. Indeed, most disease-causing mutations are private, and the availability of functional tests is limited. Therefore, distinguishing neutral mutations from disease-causing ones is challenging. This is even more challenging for rare diseases, defined in Europe as conditions with a frequency below 1: 2000, most of them being very rare. A review from Orphanet (2), revealed that the majority of rare diseases are defined by a handful of published reports describing a few individuals with a previously unidentified genetic syndrome. It is now accepted that the limitation is no longer the sequencing of one, many or all genes but rather the data analysis. In addition, while scientists were previously experts for a limited number of genes, they are now facing the ‘all genes data deluge’. This revolution has therefore resulted in a dependency on bioinformatics tools and methods to gather, store, analyze and mine the data flow. Indeed, NGS technologies typically result in the production of hundreds of millions to billions of reads per exome or genome, respectively. The analysis of these raw data can be divided into three steps as described by Gargis et al. (3). The primary analysis includes the production of sequence reads and assignment of base quality scores; the secondary analysis includes de-multiplexing, alignment of reads to a reference genome and variant calling; and the tertiary analysis is dedicated to the identification of disease-causing mutations. It involves the annotation and filtration of identified sequence variations. As reported by Salgado et al. (4) and Eilbeck et al. (5), the annotation includes various layers that should be combined in the filtration step to rapidly select a handful of candidate mutations. This filtration step can benefit from the combination of data from multiple samples as reported by Sawyer et al. (6) with Whole Exome Sequencing (WES) success rates ranging from 23% for singletons to 34% for families. To simplify this tedious process, various systems have been released such as QueryOR (7), VarElect (8), VCF-Miner (9) and BierApp (10). To annotate and prioritize mutations, these systems include multiple annotations either captured through global systems such as ANNOVAR (11) and VEP (12) or individually retrieved. However, they only partially respond to users' needs and may require a preliminary annotation step performed by bioinformaticians. In addition, for web-based solutions, confidentiality issues may arise depending on national legislation (13). In this context, we designed a new system called VarAFT (Variant Annotation and Filtration Tool), that provides a full graphical interface and includes unique features to improve mutation annotation and prioritization. It combines classical data (phylogenetic, conservation and protein structures) with additional information at variant, gene and phenotype levels. In addition, it is one of the few systems able to combine small (single nucleotide variations, small insertion/deletions) and large rearrangements (copy number variations) to get a comprehensive picture of the individual genome.

With VarAFT, users can easily annotate, filter and perform breadth and depth of coverage analysis from their data without computer programming skills and with limited hardware requirements, to efficiently identify disease-causing mutations as demonstrated in various situations (14–21).

MATERIALS AND METHODS

VarAFT is a freely available application written in Java and can therefore be used on most computers. Various binaries are available to download for Mac, Windows and Linux operating systems. VarAFT does not require any specific hardware configuration, however performances are dependent on the number of CPU cores and the amount of available memory. For example, the annotation of one sample containing 58,783 variations takes 19 min with 1 CPU core and drops to 9 min with 4 CPU cores. The breadth and depth of coverage analysis of the same sample takes 21 min with 1 CPU core and 7 min with 4. This time could be reduced to, respectively, 27 and 6 s if the user limits the breadth and depth of coverage analysis to the ACMG actionable genes list (22). Specific versions have been created to allow the installation of VarAFT on Windows machines without administrative rights. The graphical user interface was created using the Java Swing library.

The ‘coverage module’ use BEDTools (23) to compute breadth and depth of coverage data for any sequencing experiment using BAM files. Tables and charts are respectively generated with the Swing JTable library from Oracle (http://oracle.com) and the JFreeChart library (http://jfree.org). The breadth and depth of coverage analysis is performed at the genomic level unless a BED file is provided to limit the analysis to regions of interest.

The ‘annotation module’ can collect variant data from various file formats including VCF/gVCF (single or multisamples) or tabulated files. It combines all information from the dbNSFP (24) and ANNOVAR with unique features including OMIM (25), HPO (26), Gene Ontology (27), pathways (Reactome (28), KEGG (29), PID (30)) and predictions from UMD-Predictor (31) and HSF (Human Splicing Finder) (32). Note that ANNOVAR, KEGG, UMD-Predictor and HSF require a user registration to comply with their license.

Once annotated, data from any sample can be combined through the ‘filtration module’. It allows users to combine data from multiple sources and layers. The display mode can be parametered to display a subset of available columns. Interactive filtration features allow the progressive reduction of the list of candidate mutations by combining the various annotations. An in-house mutation database can be generated by VarAFT or provided by users, to exclude frequent mutations reported in a specific population and/or platform-dependent artefacts. Once filtration steps have been defined and validated, they can be saved, reapplied and shared for subsequent analysis to ensure filtration standardization in a clinical diagnosis context or for large research networks. At any filtration step, selected data can be exported for downstream analysis or reporting. Moreover, the quality of each selected mutation can be viewed in its sequencing context using IGV (33) directly from VarAFT.

RESULTS

A highly integrative system to easily pinpoint candidate disease-causing variants

As previously reported, the ability to efficiently filter genetic variation to select candidate disease-causing mutations is improved by combining data at the variant, gene and phenotypic levels (4,5). Although multiple information are available at each level, no system was able to collect and combine all this information (4). VarAFT was therefore designed to aggregate a substantially larger amount of information including the ability to combine small and large genetic variants (Table 1). In parallel, to simplify the combination of data from multiple samples, the user is able to combine samples according to predefined transmission modes, autosomal recessive or dominant, and to take into account pedigree structure (consanguinity, de novo mutations) (Figure 1). To accommodate other types of scenario, such as analysis of somatic mutations or population genetics, a custom module is also available.

Table 1.

Different annotations available within VarAFT and their availability for filtration

LEVEL	ANNOTATION	FILTRATION
VARIANT	Localization
	RefGene (Gene, Transcript, Function, Nomenclature)	X
	Ensembl (Gene, Transcript, Function, Nomenclature)	X
	Frequency
	1000 Genomes	X
	DbSNP	X
	Genome Aggregation Database (gnomAD)	X
	Known VARiants (KaViar)	X
	Haplotype Reference Consortium (HRCR1)	X
	Great Middle East Database (GME)	X
	Prediction
	UMD-Predictor	X
	Human Splicing Finder	X
	SIFT	X
	Polyphen 2	X
	LRT	X
	Mutation Taster	X
	Mutation Assessor	X
	FATHMM
	PROVEAN	X
	VEST3
	MetaSVM
	MetaLR
	M-CAP	X
	CADD	X
	DANN	X
	FATHMM-MKL
	Eigen	X
	GERP++	X
	Conservation
	phyloP100way
	phyloP20way
	phastCons100way
	phastCons20way
	SiPhy_29way
	wgRNA
	predicted microRNA targets (target ScanS)
	genomicSuperDups
	gwasCatalog
	wgEncodeBroadHmmGm12878HMM
GENE	Expression
	GTEx	X
	Pathways/GO
	KEGG	X
	Pathway Interactome Database (PID)	X
	REACTOME	X
	Gene Ontology	X
	Score Tolerance
	Residual Variation Intolerance Score (RVIS)	X
	Gene Damaging Index (GDI)	X
	Lost Of Function Tool (LOF)	X
	genome-wide haploinsufficiency score (GHIS)	X
DISEASE PHENOTYPE	Ontology
	Human Phenotype Ontology (HPO)	X
	Database
	Online Mendelian Inheritance in Man (OMIM)	X
	Catalogue Of Somatic Mutations In Cancer (Cosmic)	X
	ClinVar	X

X = annotation available for filtration.

Figure 1.

Analysis and Filtration module. (A) Top: main screen for filtration and analysis of variants. Top part: easy combination of samples (singleton, trio or any combination) for AD, ARD or other modes of inheritance. Middle: filtration parameters. Bottom: list of annotated variants. This list is dynamically updated based on filtration criteria. (B) Filtration criteria are divided into five sections (variant type, frequency, pathogenicity predictions, gene information and others) including multiple parameters that can be combined by the user. (C) The selection of one variant from the list gives access to additional data related to: general information linked to the variant itself, the presence of this variant in the analysed samples/patients, its impact on transcripts and related HGVS nomenclature (RefGene and Ensembl), its reported frequency in general populations, the prediction of its pathogenicity from multiple tools, the tolerance to loss of function of the related gene, the tissue expression pattern of the gene/transcripts, additional information related to chromosomal region as promotors, regulatory regions, etc. and access to various useful external websites.

Experiment quality control compatible with clinical use

The ‘coverage analysis module’ was designed to evaluate experiment quality. It provides the breadth and depth of coverage for any transcript or exon at the nucleotide level, either through dynamic histograms or tables (Figure 2). A report can be generated to rank genes and exons according to their breadth of coverage at a depth of 1, 5, 10, 20 or 30× and evaluate their quality in accordance with international guidelines (EuroGentest: www.eurogentest.org). In a clinical diagnosis context, BED files can be provided or generated through VarAFT to restrict the exome analysis to some transcripts or genes. Indeed, as indicated in the Eurogentest guidelines, to limit incidental findings it is recommended to focus on genes of interest for which a relationship between genotype and phenotype has been published and confirmed (34).

Figure 2.

Breadth and depth of coverage analysis of a WES experiment. (A) Top: main screen displaying the list of gene symbols and associated RefSeq transcripts, number of exons and coding sequence size with the following statistical values: mean depth ± standard deviation and breadth of coverage at 30×, 20×, 10×, 5× and 1× depth. Bottom-left: breadth of coverage for each exon of the selected transcript at a selected depth (20× for transcript NM_001244910 of the FCGR1B gene); Bottom-right: breadth of coverage at the nucleotide level of a selected exon (exon 1 of the NM_001244910 transcript). (B) Histogram displaying the percentage of transcripts with a breadth of coverage superior or equal to a depth of 1× (red), 5× (blue), 10× (green), 20× (yellow) and 30× (purple). (C) Breadth of coverage representation of the various exons from all transcripts of the selected gene (MEN1). Color code as for Figure 1A and C: red = <10× depth of coverage; yellow = 10–20× depth of coverage; blue ≥ 20× depth of coverage.

Use case demonstrating VarAFT efficiency in various situations

VarAFT has been extensively used in both clinical diagnosis and research contexts. Disease-causing mutations were identified from trio, cohorts and individual cases from autosomal dominant and recessive diseases, such as dystonia, neuromuscular disorders, mental retardation and premature aging. For example, VarAFT was recently used to analyze 306 genes in a cohort of distal myopathy patients (35). The software was also evaluated as a prioritization system to highlight mutations involved in cancers (14) and other situations (14–21).

To demonstrate VarAFT usefulness and efficiency, we chose the following use cases:

Use case #1. The first dataset was extracted from Kamphans et al. and contains VCF files from four individuals (464, 465, 466 and 467) (36) from a family with an autosomal recessive disease. Each sample contained respectively 20 708, 20 560, 20 552 and 20 547 variants. The VarAFT processing for this family included the following five steps: (i) combination of data from the various samples taking into account the mode of inheritance (AR) to select compound heterozygous in the affected individuals (464 and 465) and present in only one parent (466 and 467); (ii) selection of variants localized in exons or bordering introns (±4 bp); (iii) exclusion of variants with a frequency in general populations above 1%; (iv) selection of variants predicted as pathogenic by CADD (37) and (v) selection of variants predicted as pathogenic by the UMD-Predictor system. Note that steps can be conducted in any order and will lead to the same result. The two remaining variants correspond to the two mutations, c.2869C>T (p.Leu957Phe) and c.2355dupC (p.Gly785fs), in the PIGO gene, identified by the authors as the disease-causing mutations in this family (Figure 3).

Figure 3.

Use case#1. The identification of the disease-causing mutations from the AR family described by Kamphans et al. (36) including data from four members. It was performed in five steps using mainly data from the variant annotation layer. * gnomAD; 1000 genomes; KaViar; HRCR1 and GME databases; ** polymorphism and probably polymorphism were excluded.

Use case #2 was extracted from Miltgen et al. and contains VCF files from four patients from a multigenerational family from Flemish origin with Craniocervical Dystonia and an autosomic dominant mode of inheritance (20). Each sample contained, respectively, 66 215, 58 783, 58 959 and 59 495 variants. The VarAFT processing for this family included the following seven steps: (i) combination of data from the various samples taking into account the mode of inheritance (AD) to select heterozygous variants in the affected individuals (D7, D8 and D11) and absent in D10; (ii) selection of variants localized in exons or bordering introns (±4 bp); (iii) exclusion of variants with a frequency in general populations above 1%; (iv) selection of variants predicted as pathogenic by CADD; (v) selection of variants predicted as pathogenic by the UMD-Predictor system; (vi) selection of genes expressed in the affected tissue (brain) and (vii) selection of genes associated with at least one (‘OR’ option) HPO term describing symptoms found in patients (dystonia, Blepharospasm and torticollis). Note that steps can be conducted in any order and will lead to the same result. The two remaining variants correspond to the mutation c.240+1G>T from the SURF1 gene and c.1969G>A (p.Ala657Thr) from the ANO3 gene. The SURF1 mutation was excluded as this gene is only involved in autosomal recessive diseases: Charcot-Marie-Tooth disease, type 4K (MIM #616624) and Leigh syndrome, due to COX IV deficiency (MIM #516000). In contrast, mutations from the ANO3 gene have been reported in the autosomal dominant Dystonia 24 and this mutation was identified by the authors as the disease-causing mutations in this family (Figure 4).

Figure 4.

Use case#2. The identification of the disease-causing mutations from the AD family described by Miltgen et al. (20) including data from four members. It was performed in seven steps using data from the variant annotation and the phenotype layers. * gnomAD; 1000 genomes; KaViar; HRCR1 and GME databases; ** polymorphism and probably polymorphism were excluded.

Use case #3 is an artificial VCF created by inserting the SH3TC2 compound heterozygous c.[279G>A]; [805+2T>C] disease-causing mutations identified by Piscosquito et al., 2016 (38) into a personal exome (https://personalgenomics.zone). The sample contained 37 694 variants. The VarAFT processing for this sample included the following five steps: (i) mode of inheritance (AR) to select only compound heterozygous variants; (ii) selection of variants localized in exons or bordering introns (±4 bp); (iii) exclusion of variants with a frequency in general populations above 1%; (iv) selection of variants predicted as pathogenic by the UMD-Predictor system; and (v) selection of genes associated with at least one (‘OR’ option) general HPO term associated with the pathology (‘Decreased number of large peripheral myelinated nerve fibers’ HP:0003387 or ‘Segmental peripheral demyelination’ HP:0007107). Note that steps can be conducted in any order and will lead to the same result. The two remaining variants correspond to the mutations c.2860C>T (p.Arg954*) and c.279G>A (p.Lys93Lys) from the SH3TC2 gene and identified by the authors as the disease-causing mutations in this family (Figure 5). Note that this mutation was predicted as pathogenic only by the UMD-Predictor system as it impacts the donor splice site. This was confirmed by predictions from the HSF system.

Figure 5.

Use case#3. The identification of the disease-causing mutations from an artificial single VCF corresponding to a proband with an AR. It was performed in five steps using data from the variant annotation and the phenotype layers. * gnomAD; 1000 genomes; KaViar; HRCR1 and GME databases; ** polymorphism and probably polymorphism were excluded.

As illustrated in the three use cases, VarAFT was flexible enough to apply optimal filtration criteria taking into account the mode of inheritance and the available phenotypic information. In each case, the process resulted in the identification of the disease-causing mutations (Figures 3–5). Only a subset of information was used in the processes and additional features (pathways, tissue expression, etc.) are available for more complex situations.

DISCUSSION

VarAFT is a multiplatform freely available software that allows the simultaneous annotation, filtration, and breadth and depth of coverage analysis of WES, WGS and targeted sequencing experiments from any sequencing platform. Its graphical user interface, various modules and unique features, such as pathogenicity predictions from UMD-Predictor and HSF, allow untrained users to rapidly highlight disease-causing mutations in multiple genetic scenarios. In addition, VarAFT allows visualization of data quality (breadth and depth of coverage) through direct access to BAM files. The nucleotides and genotypes can also be easily accessed through IGV.

As reported by Salgado et al. (4), on one hand, automatic prioritization systems are now available to ensure a homogeneous treatment of samples. However, these systems are based on previously established links between genotype and phenotype and can only solve a limited number of diagnosis and research problems. On the other hand, manual systems are numerous and heterogeneous in their content and filtration features. VarAFT was tailored to overcome identified limitations such as the access to multiple layers of information, the combination of small (SNV) and large (CNV) mutations in a single analysis and the accessibility for all scientists, regardless of bioinformatics skill level. So, users can rapidly end up with shorter and more accurate lists of candidate disease-causing mutations, facilitating downstream validation, gene discovery and genetic counseling. As a standalone application, it can be used for clinical diagnosis as data are processed locally avoiding network privacy issues.

VarAFT uses recognized resources, formats and ontologies making it suitable for integration in any NGS environment. VarAFT was presented to the scientific community in various international training courses (RD-Connect, 3Gb-Test, Variant Effect Predictor training course from the Human Variome Project, ELIXIR training course on variants analysis) and was rapidly adopted by 800 users from more than 50 countries. At last, VarAFT was instrumental in the creation of the RD-Connect Genome-Phenome analysis platform (39).