High-throughput detection of mutations responsible for childhood hearing loss using resequencing microarrays

Overall array performance

The Cincinnati array also contained 8 genes, 3 of which were common to the Harvard array (see Methods). We ran 12 arrays in the pilot batch and characterized array performance. Base call rates for these 12 arrays (315,504 bases) ranged from 82.5% to 96.9% with an average call rate of 91.3% using GSEQ 4.0 (Figures 1 & 2). Dideoxy sequencing was performed for 296,296 of the bases and comparison of this data with array calls gave a call accuracy of 99.23% with nearly 180 false positives per array.

Improved array performance with protocol optimization

When we compared data from Cincinnati and Harvard arrays, the number of no-calls and false positives from the former were found to be higher. Early data from the Harvard arrays had shown that call rate worsened when the fragmentation was incomplete and that the bases most affected were those within the long range PCR fragments. In addition, the reproducibility of the quantity of product from long range PCR was less. Based upon these two factors, the target amplification process for Cincinnati arrays was modified. Instead of using a combination of long and short range PCR (54 fragments ranging from 434 to 13,282 bases in length), as was employed in the pilot batch, all long range PCRs were converted to short range resulting in 180 fragments ranging from 315 to 980 bases in length. The impact of using shorter PCR products was evaluated by comparing array data across the two protocols (Figures 1 & 2). A total of 13 arrays (341,796 bases) were run with the optimized "short range PCR only" protocol and an average call rate of 97.9% (range 96.7% to 98.4%) was obtained with GSEQ (Table 1). Dideoxy sequencing was performed for 336,171 of the bases and comparison of this data with array calls gave an array call accuracy of 99.83%. The average number of false positives dropped from 180 to 42 per array. The modified protocol with shorter PCR products was then adopted for subsequent arrays.

Detection of insertions and deletions

Through dideoxy sequencing, we identified four cases with at least one insertion or deletion (35delG/35delG, 35delG/167delT, and M1V/167delT in GJB2 and homozygous 1180_1187del8ins(β-sat) in TMPRSS3), of which 2 had been previously reported and therefore had probes tiled on the array for their detection (GJB2 35delG and 167delT) (Additional file 1). We analyzed array data to look for no-calls or variant calls in the vicinity of the indel sites but did not observe any such patterns for the 35delG and 167delT alleles. It should be noted that the local high GC-content surrounding the 35delG (stretch of 5 Gs; probe GC 56%) and the 167delT (stretch of 4 Cs; probe GC 64%) would be expected to make detection extremely difficult. On the other hand, there was a continuous stretch of 13 no-calls and a variant call spanning the deleted bases of the TMPRSS3 gene that led to the detection of the mutation. We also analyzed raw feature intensities within fragments to see if indels cause degradation in intensities surrounding the variant site. We observed lower peak intensities surrounding the TMPRSS3 mutation but did not find such evidence for either of the single base GJB2 deletions (data not shown). GDAS/GSEQ are currently not designed for identifying indels so their low detection rate was an expected observation.

Differential impact of high C-content and G-content on probe performance

In agreement with previously reported findings [19, 20], we observed that an increase in probe G- and C-contents have differential impact on performance. We used complementary feature quartets to determine intensities associated with the C- or G-content of a probe. Intensity characteristics varied differently with respect to an increase in C- or G-content and average peak intensity was affected more severely by a high G-content than by an equally high C-content (Figure 3a). In order to assess if our data demonstrated previously reported debilitating effect of stretches of Gs (G-stacks) on probe performance [21], we compared hybridization intensities among probes with the same G-content grouped based on the presence or absence of G-stretches (≥ 4 Gs in a continuous stretch). For the same G-content, probes with G-stretches produced lower peak intensities than probes with C-stretches or without any stretches of Gs or Cs (Figure 3b).

Overall array performance with sPROFILER

sPROFILER was only used on no-calls from GDAS/GSEQ. Examined bases were called wild-type based on single-strand evidence or were left as no-calls depending upon the feature intensity profile. Base calls were subsequently compared against GSEQ and dideoxy sequencing calls. Table 1 and Figures 1 & 2 provide detailed comparison of call rates, number of false positives and false negatives before and after analyzing GSEQ calls with sPROFILER. For Cincinnati arrays run with short and long range PCR fragments, the average call rate increased to 96.7% (99.82% call accuracy). The average number of no-calls dropped from 2350 to 902 per array and the number of false positives dropped from 180 to 44 bases per array. Improvement was obtained at the cost of incorrectly assigning an additional 0.38 true variants per array as wild-type. For the optmized Cincinnati arrays, average call rate increased to 99.6% (range 99.5% to 99.8%) with 99.88% call accuracy and number of no-calls dropped from 563 to 103 (Table 1). The number of false positives dropped from 42 to 30 per array with the application of filters to screen variant calls based on low quality scores and the number of no-calls and variant calls in the vicinity. Improvement was achieved at the cost of an additional 0.4 false negative per array. We also calculated the number of exons per array that would need to be interrogated with dideoxy sequencing in order to resolve no-calls and confirm positive calls. After processing GSEQ calls with sPROFILER, the average number of exons to be sequenced in order to clarify no-calls and confirm variants dropped from 150 (range 132 to 165) to 68 (range 52 to 82) per array.

When using sPROFILER on the Harvard arrays, the average call rate increased from 96.9% to 99.6% and number of no-calls dropped from 781 to 101 per array (Table 1). False positive calls could be reduced from 41 to 38 per array with the filters for false positive calls. An additional 0.4 variants were falsely called wild-type using the algorithm. Because all no-calls are followed up by dideoxy sequencing, the implementation of the sPROFILER resulted in substantially fewer exons (average of 153 down to 52) needing follow-up.

Clinical sensitivity of arrays

After excluding synonymous and/or common variants, 18/37 patients examined by the Cincinnati array and 9/24 patients examined on the Harvard array had at least one variant of potential or likely clinical significance (excludes two cases known to have GJB2 biallelic mutations). See Additional file 1 for a full list of all variants found on a per patient basis. To further assess the likelihood that each variant may be disease-causing, efforts are currently underway to examine these variants in further control studies and using in silico algorithms to predict protein impact. After these studies are completed, variants with potential clinical significance will be confirmed in a CLIA environment with results returned to patients under an IRB approved protocol.

Patient enrolment

All 74 patient samples and clinical information were collected under IRB-approved protocols. Most cases (66) were recruited from clinical centers: 25 patients from Ohio (J. Greinwald), 20 from Massachusetts (H. Rehm and M. Kenna), 12 from Belgium (G. Van Camp), 5 from Israel (K. Avraham), 2 from Iowa (R. Smith), 2 from Nebraska (P. Kelley), and the remaining 8 represented anonymized samples obtained through clinical testing at the Laboratory for Molecular Medicine at the Partners Center for Personalized Genetic Medicine. All patients had previously tested negative for biallelic GJB2 mutations, except for 3 patients who had heterozygous mutations in the GJB2 gene. All patients had bilateral SNHL and were either recessive cases (e.g. with affected siblings) or singletons. DNA and clinical information was collected on most patients (de-identified clinical information for those cases that were anonymized) and included demographic information, audiometric profiles, family history, neuro-otologic history and physical examinations.

Gene selection and array design

To determine what genes to place on a finite resequencing platform that was initially capable of accommodating less than 27,000 bp, we undertook a comprehensive literature search to estimate the likely frequency of mutations in genes known to cause childhood SNHL. Based on the prevalence of mutations in published populations and family studies, the inheritance pattern(s) of mutations (i.e. preference for recessive childhood hearing loss) and the impact of the gene on patient outcome (e.g. blindness developed due to Usher Syndrome), an initial set of genes were selected for inclusion on one or both hearing loss arrays (e.g. GJB2, GJB6, CDH23, KCNE1, KCNQ1, MYO7A, OTOF, USH1C, and TMPRSS3). Due to space limitations only the most conserved portion of CDH23 was tiled. Most other genes had extremely limited data on relative contribution and therefore several genes were selected based upon other minor characteristics (e.g. contribution to both dominant and recessive hearing loss (e.g. MYO6), interest in discovering a larger undocumented role in hearing loss (e.g. SLC26A5) and/or small gene size to fill remaining capacity on array (e.g. TMIE). It should be noted that current studies suggest that SLC26A5 mutations may not be a cause of hearing loss [34]. Details of the sequences included are in the Methods and Additional file 1.

Two arrays were designed, one at University of Cincinnati/Children's Hospital Medical Center (Cincinnati, OH) and one at Harvard Medical School (Boston, MA) with 8 genes on each array representing a total of 13 unique genes assessed. The Cincinnati array contained: GJB2, GJB6, CDH23 (59 out of 69 exons; 80.3% of the coding sequence), KCNE1, KCNQ1, MYO7A, OTOF, and SLC26A4 genes and the Harvard array contained: GJB2, MYO6, MYO7A, OTOF, SLC26A5, TMIE, TMPRSS3, and USH1C genes. For all genes, both the coding sequence and splice junctions were assessed. For the Cincinnati array, 2 bp of each flanking splice site were tiled for most exons. For the Harvard array, 10 bp of each flanking splice site were tiled for all exons. Additional file 1 contains an overview of the genes tiled on the arrays. Each array contained probes to interrogate roughly 26,000 bases of DNA (8 probes representing all 4 bases to assess the forward strand and all 4 bases to assess the reverse strand). In addition, probes designed to detect 17 previously reported insertions and deletions were also tiled on the Harvard array (Additional file 1).

Array protocol

Initially the arrays were run according to the manufacturer's protocol (Affymetrix, Santa Clara, CA). Briefly, long range PCR conditions for the LA TaKaRa Polymerase (Takara, Japan) were: TaKaRa LA Taq 0.05 U/ul, 1× LA PCR Buffer II, 400 uM (each) dNTPs, 0.3 uM (each) primers, 4 ng/ul genomic DNA in a 25 ul reaction volume. Short range PCR conditions are described below for dideoxy sequencing. Cycling conditions for most reactions were 94°C 2 min, (94°C 15 sec, 68°C 9 min) × 30, 68°C 14 min, 4°C. Modifications using standard approaches to PCR optimization were made for some difficult reactions. All PCR assays were either quantified using PicoGreen (Molecular Probes, Eugene, OR) or absorbance spectrometry and then pooled in equimolar amounts. The PCR products were then purified, fragmented, labeled and hybridized to the array. Finally, the arrays were washed and scanned and the data were analyzed using the GeneChip DNA Analysis Software (GDAS 2.0, Affymetrix) or the GeneChip Sequence Analysis Software (GSEQ 4.0, Affymetrix). After initial experience, the Cincinnati protocol was modified to include only short range PCR in order to increase PCR and fragmentation consistency and facilitate the use of automation. In addition, PCR products were cleaned robotically using magnetic beads, 7.5× less fragmentation reagent was used than suggested (0.02 U vs the recommended 0.15 U) and fragmentation products were analyzed using an Agilent Bioanalyzer (Santa Clara, CA).

Dideoxy DNA sequencing protocol

Complete dideoxy sequencing of all exons was performed for 14 of the Harvard patient arrays and all of the 25 Cincinnati arrays (enabling assessment of false negative rates) and partial sequencing was performed for the remaining 35 patients in order to confirm or clarify bases assigned as variants or "no calls" by Affymetrix GDAS or GSEQ software. PCR conditions: AmpliTaq Gold (Applied Biosystems, Foster City, CA) 0.05 U/ul, 1× ABI PCR Buffer, 2.5 mM MgCl2, 400 uM (each) dNTPs, 0.4 uM (each) primers, 2 ng/ul genomic DNA in a 25 ul reaction volume. Modifications using standards approaches to PCR optimization have been made for some difficult reactions. Cycling conditions are 94°C 2 min, (94°C 15 sec, 60°C 30 sec, 72°C 30 sec) × 30, 72°C 10 min, 4°C. Sequencing reactions are performed using the same primers for short-range PCR (or internal primers for long-range) in the following conditions: 0.25 μl Big Dye 3.1 (Applied Biosystems), 3.75 μl 5× Big Dye Buffer, 0.50 μl DMSO, 2 μl 10 μM primer, 2 ul template (1-5 ng) in 10 μl total volume. Cycle sequencing conditions are 94°C 4 min, (98°C 30 sec, 50°C 5 sec, 60°C 4 min) × 30, 4°C. Reactions are cleaned up using CleanSEQ magnetic beads (Agencourt Biosciences Corp, Beverly, MA) and then run on an ABI 3730 DNA Analyzer (Applied Biosystems). Chromatograms are analyzed using an automated Phred analysis program to check for quality followed by analysis using Sequencher 4.0 (GeneCodes, Ann Arbor, MI) and/or Mutation Surveyor (SoftGenetics, State College, PA) software.

Data analysis

Affymetrix GDAS 2.0 with optimized settings (Additional file 1) was used for assignment of base calls for Harvard arrays and GSEQ 4.0 was used for Cincinnati arrays. The algorithm settings were determined after evaluating different quality score thresholds and comparing array calls with dideoxy sequencing results to characterize coverage and accuracy. The settings chosen were optimized for the highest call rate without increasing the false negative rate. We also developed a novel algorithm (sPROFILER) for strand-specific probe cell intensity comparison for filtering GDAS/GSEQ base calls. The algorithm was developed to re-examine GDAS/GSEQ output to reduce no-calls and, in some cases, false positives. sPROFILER attempts to resolve no-calls based on intensity signature from a single strand when a base cannot be called due to poor hybridization on one of the strands. The algorithm uses all wild-type bases within the array to determine threshold for peak to next highest intensity ratio on either strand and uses the base call at the position of interest across all arrays to determine the proportion of wild-type calls that are made on that position by GSEQ. The latter is used for scaling the threshold ratio and thus, in effect, making the algorithm more conservative while attempting to assign a wild-type call to a position that is being called variant in a large number of samples and vice-versa. We also adopted two additional bioinformatics filters developed for the Francisella tularensis whole-genome resequencing platform [17]. The filters were designed for reduction in false positives by screening variant calls that 1) are in regions rich in variant calls and no-calls, and 2) have low quality scores for the corresponding base call. With the exception of wild-type calls within no-call stretches, sPROFILER does not attempt to re-examine any base calls conforming to the reference sequence because GDAS and GSEQ are conservative in making wild-type calls and thus achieve low false negative rates. The output calls were compared against sequencing results and against GDAS/GSEQ calls for validation. sPROFILER was implemented in MATLAB. A detailed description of sPROFILER and accompanying MATLAB code are provided in Additional files 1 & 2.