High-Throughput Genetic Testing for Thrombotic Microangiopathies and C3 Glomerulopathies

Subjects

During the 11-month period beginning on January 1, 2014, 193 patients were screened using the Genetic Complement-Mediated Renal Disease Panel (GRP), a TGE&NGS panel that includes CFH, CFI, CFB, C3, CFHR5, CD46, DGKE, ADAMTS13, THBD, PLG, and CNV analysis of CFHR3-CFHR1. The GRP process and variants interpretation are illustrated in Figure 1 and Supplementary Figure 1. We included all patients on whom testing was requested with clinical information and classified patients based on the requesting physician’s clinical impression. Physician-based clinical diagnoses included: (1) TMA: 118 aHUS, six TTP and 11 other TMA patients; (2) C3G: 30 C3GN and five DDD patients; and (3) Other: nine patients with untargeted diseases, such as systemic lupus erythematosus (SLE) and IgA nephropathy (IgAN). Fourteen patients with an ambiguous diagnosis (12 aHUS/TTP patients and two C3GN/DDD patients) were also included (Table 1). The Institutional Review Board of Carver College of Medicine at the University of Iowa approved this study.

Validating the GRP

Prior to clinical implementation, the GRP was validated using 92 archived and de-identified TMA or C3G patients who underwent Sanger sequencing for one or more of the following genes: CFH, CFI, CFB, C3, CFHR5, CD46, and THBD. In aggregate, in the validation process, we compared genotype data from 55 exonic and 14 intronic common/rare variants. For NGS, samples were divided into two batches, each of which was sequenced on one lane of an Illumina HiSeq 2000 (Illumina Inc., San Diego, CA). To test reproducibility, three samples were sequenced twice on the HiSeq instrument by including them in both batches, and 13 samples were sequenced on both the HiSeq and MiSeq (Illumina Inc.) instruments. In all cases, results of duplicate runs were 100% concordant. Average Qvar (Phred-like quality score), Depth, and QD (Phred-like quality score divided by depth) of variant callings were 20596, 1172, and 18.30, respectively.

As shown in Figure 2, in total, 622 high-quality (QD>8) and 26 fair-quality (5<QD≤8) variants were Sanger-validated with 100% and 65% concordance. Four of 19 (21%) poor-quality variants (QD<5) were Sanger-confirmed. The positive predictive value based on this validation process was 96%. We observed a low-quality gap in exons 20 and 21 of CFH. Specifically, c.3138C>T (rs61822181) and c. 3150T>C (rs113347629) in exon 20 were false-positives and had QDs between 1.7 and 6.8; and c.3572C>T, p.Ser1191Leu (rs460897) in exon 22 was a false-negative and was omitted in the filtering process due to a very low QD score of 0.09 (Figure 2). Other low-quality gaps were in introns of CFI and C3, remote from targeted coding sequence and splicing sites, and were considered insignificant. On the basis of these results, we retain all variants with QD>5 and Depth>10 and Sanger confirm positive findings. In addition, we always Sanger sequence exons 20–22 of CFH.

NGS Summary

NGS was completed on 193 samples in 50 runs of varying sample number on the MiSeq. The average number of total reads per sample was 4.01 million, with a range from 2.00 to 11.89 million (Supplementary Figure 2). Average coverage was 723× with 99.37% of the total target region covered by more than 30 reads. The lowest coverage (297×) was associated with low total reads (2.00 million). Homozygosity for del(CFHR3-CFHR1) was associated with a reduction in 30× coverage to less than 94% in all 17 homozygous deletion carriers.

Variants Calling and Filtering

A mean of 217 variants were identified per patient, which then were filtered by quality (QD>5 and Depth>10), minor allele frequency (MAF), and nucleotide effect. In total, we identified 17 novel variants not reported in the literature or any database (see CONCISE METHODS) and 74 rare variants (MAF<1% in all populations) (Supplementary Table 1, Table 2). All novel and rare variants were Sanger-validated (100% confirmation).

Detecting the CFHR3-CFHR1 CNV

The CFHR3-CFHR1 CNV was detected by multiplex ligation-dependent probe amplification (MLPA). Homozygosity for this deletion was identified in approximately 12% of aHUS patients (Table 3), a prevalence significantly higher than that seen in 314 healthy controls (3%; Fisher’s exact test, P<0.01). All del(CFHR3-CFHR1) homozygotes were also identifiable by NGS as total target region coverage was reduced to <94% at 30× (Supplementary Figure 2). We were able to identify del(CFHR3-CFHR1) heterozygous in validation batches sequenced on the HiSeq; however, their unambiguous identification on the MiSeq was not possible due to smaller pool size and lower total coverage (data not shown).

Variant Interpretation

To provide a clinically relevant report, a multidisciplinary board (Renal Group Meeting) reviewed all genetic data in the context of the available clinical data (Figure 1, Supplementary Figure 1). Variants with MAF>1% known to be unrelated to disease were classified as benign. Ultra-rare variants reported as pathogenic in the literature with supporting functional evidence were classified as pathogenic. Novel or rare variants that change protein sequence but have an unknown impact on protein function were classified as either likely pathogenic or variant of uncertain significance (VUS), a distinction that reflected two additional calculations. Likely pathogenic variants were also: (1) missense variants with pathogenicity scores ≥5 (based on GERP++,¹⁷ PhyloP,¹⁸ MutationTaster,¹⁹ PolyPhen2,²⁰ SIFT,²¹ and LRT²²); ultra-rare or disease-associated (based on a comparison of reported MAF to disease-associated prevalence in the Renal Variant Database, an in-house variants database of 1085 de-identified patients with TMAs or C3Gs); and found in disease-related functional domains/loci; or (2) novel and caused loss of function. Note that known pathogenic and likely pathogenic variants have similar pathogenicity scores and that likely pathogenic variants have higher pathogenicity scores than VUSs (Supplementary Figure 3). Based on genotypic findings and the clinical phenotype, additional testing was occasionally recommended. The Renal Group Meeting also served as the forum for discussion of all email inquiries regarding patient management of these ultra-rare diseases.

Variant Burden and Distribution

Variant burden per disease group was calculated as the sum of VUS/likely pathogenic/pathogenic variants divided by the number of patients. For example, in the group of 147 patients with TMA, eight VUSs, three likely pathogenic, and three pathogenic variants were identified in CFH, for a rare variant load of 0.095. As shown in Figure 3A, the variant load in CFH, CFI, and C3 was high in both TMA and C3G patients, while for CD46 and CFB, the variant load was high in only TMA or C3G, respectively. ADAMTS13 variants were found mainly in patients with aHUS or TTP (variant load, 0.073), suggesting a common mechanism driving these two diseases. The variant load in DGKE was high in C3G patients but all DGKE variants were present in heterozygosity; the function impact of this finding is not known. A single patient, a newborn with aHUS, was homozygous for a DGKE variant (c.465–2A>G) predicted to disrupt exon splicing.

The distribution of VUS/likely pathogenic/pathogenic variants was significantly different between C3G and TMA patients at the level of functionally grouped genes. As shown in Figure 3B, 16.2% (n=6), 13.5% (n=5), and 5.4% (n=2) of patients with C3G, respectively, carried variants only in C3 convertases genes (C3 and CFB), only in AP regulators (CFH, CFI, CD46 and CFHR5), or in both gene groups, as compared with 5.4% (n=8), 21.1% (n=31), and 0.7% (n=1) of patients with TMA. We also assessed single/combined variants per patient and found that 43.2% C3G and 36.7% of patients with TMA carry at least one VUS/likely pathogenic/pathogenic variant, both prevalences significantly higher (P<0.05) than controls. As show in Figure 3C, the prevalence of multiple VUS/likely pathogenic/pathogenic variants per patient was also higher in the C3G group as compared with the TMA group and controls (13.5% [n=5] versus 6.8% [n=10] versus 4.1% [n=103], respectively; P<0.05). Control data were retrieved and filtered from the 1000 Genomes Project (1000Gs), Phase 3 (n=2504), an unphenotyped control group unlikely (P<0.01) to include patients with TMA or C3G.

Positive Rate

Identification of VUS/likely pathogenic/pathogenic variants was a significant finding. Additionally, we also considered the deletion of CFHR3-CFHR1 as a positive result in patients with aHUS and recommended screening for anti-CFH autoantibodies if that test had not been considered.²³ Across all patient categories (aHUS, TTP, C3GN and DDD), we had a diagnostic rate of approximately 40%; however, if we consider only CFH, CD46, CFI, CFB, and C3 in the aHUS group of patients (118 patients), the positive rate was 27.1%, which is lower than the commonly reported rate of approximately 45%.²⁴ The diagnostic rate was also low in patients with other diagnoses (Figure 3D). More than half of patients had a negative genetic screen, suggesting the presence of unrecognized genetic and/or acquired drivers of disease.

Molecular Modeling

To better understand the functional impact of variants in C3 convertase genes, refined molecular modeling was completed for VUS/likely pathogenic/pathogenic variants identified in C3 and factor B. Using the AMOEBA polarizable force field as a part of the Force Field X software package, four of the five C3 variants and both factor B variants identified in C3G patients were predicted to have molecular interactions that might impact convertase function and/or regulation. Two of the four C3 variants in TMA patients were also predicted to destabilize folding (Supplementary Figure 4, Supplementary Table 2).

Sample Preparation, Target Genomic Enrichment, and Next Generation Sequencing

For each sample, genomic DNA was extracted from peripheral blood using the Gentra Puregene Kit (Qiagen Inc., Valencia, CA) and integrity was evaluated by 1% agarose gel electrophoresis. The absorbance at 230:260:280 was measured using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) to ensure DNA samples met our minimal quality metrics of 1.8 for 260/280 and 260/230 ratios. DNA concentration was determined using the Qubit dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA).

For each GRP gene, the coding sequence and flanking splice sites were captured using the Agilent SureSelect Target Enrichment System (Agilent Technologies, Santa Clara, CA). Library preparation was performed with SureSelect TGE baits and SureSelectXT Reagent Kits (Agilent Technologies) following the manufacturer’s protocol. Preparation was automated using a Zephyr Workstation (PerkinElmer, Waltham, MI). Library quality and concentration were evaluated using a Bioanalyzer 2100 (Agilent Technologies). Libraries passing this quality control step were pooled and sequenced with a 100-bp paired-end module in one lane on a HiSeq 2000 (Illumina Inc.) or a 150-bp paired-end module on a MiSeq Sequencer (Illumina Inc.).

Next-Generation Sequencing Data Analysis

NGS data storage and analysis were performed on dedicated computing resources maintained by the Iowa Institute of Human Genetics at the University of Iowa. Sequencing data were archived as fastq files on a secured storage server and then analyzed using locally implemented open-source Galaxy software on a high-performance computing cluster. The workflow for variant calling integrated publicly available tools: reads were mapped using Burrows–Wheeler Alignment against human reference genome GRCh37/hg19; duplicates were removed by Picard; realignment, calibration, and variant calling were performed with GATK; variant annotation was performed with a CLCG Annotation and Reporting Tool developed by our bioinformatics team.

Variant Prioritization and Sanger Validation

Total number of reads per sample varied as a function of the number of samples per run and DNA input per sample. Low-quality variants (Depth<10 or QD<5) were filtered out by quality control. Common variants with MAF>1% in any population were excluded (based on the NHLBI GO Exome Sequencing Project [ESP, evs.gs.washington.edu], the 1000 Genomes Project [1000Gs, www.1000genomes.org], and most recently, the Exome Aggregation Consortium [ExAC, exac.broadinstitute.org]). Variants also were filtered based on predicted effect, retaining nonsynonymous single nucleotide variants, canonical splicing changes, and indels. These variants were prioritized based on MAF, nucleotide conservation, reported functional/expressional impact, and phenotype correlation. Other reference databases routinely queried included the aHUS Mutation Database (www.fh-hus.org), Human Gene Mutation Database, and our in-house Renal Variant Database. GERP++,¹⁷ PhyloP,¹⁸ MutationTaster,¹⁹ PolyPhen2,²⁰ SIFT,²¹ and LRT²² were used to calculate variant-specific pathogenicity scores, which we based on the sum of tools predicting a given variant to be deleterious. All reported variants were Sanger-validated.

Copy Number Variation

CNVs across the CFHR3-CFHR1 region were identified using an MLPA set of 13 probes and six control probes, all designed following the MRC-Holland synthetic probe design protocol. Patients with deletion in CFHR3-CFHR1 were further tested for CNV of CFH, CFHR4, CFHR2, and CFHR5. MLPA was performed using the SALSA MLPA Reagent Kit (MRC-Holland, Amsterdam, The Netherlands), resolved on a 3130xl Genetic Analyzer (Life Technologies), and analyzed with GeneMapper software (Life Technologies). In each MLPA run, we included eight control samples (five normal controls [three females, two males], two controls with a heterozygous deletion of CFHR3-CFHR1 [one female, one male], and one control with a homozygous deletion for CFHR3-CFHR1).

Molecular Modeling

Homology models were acquired and refined using the AMOEBA polarizable force field as a part of the Force Field X software package.^30,31 The model refinement consisted of local minimization followed by rotamer optimization around the mutation and then a second minimization step. The first minimization step eliminated obvious steric clashes in the protein; rotamer optimization allowed side chain atoms of residues near the mutation to be altered into a specific set of discrete conformations (rotamers) with low energy,³² and the final minimization step allowed rigid conformations in side chains to relax. This protocol was used to refine the wild-type model to remove model bias before modeling the identified variants; wild-type and mutant models were superimposed using the PyMOL molecular visualization program.

Statistical Analyses

Data were analyzed using R (v2.15.1). Fisher exact test was used to analyze categorical data. All tests were two-tailed, and P values less than 0.05 were considered significant.