Genetic diagnosis by whole exome capture and massively parallel DNA sequencing

The Roche/NimbleGen whole-exome array capture protocols were developed for DNA sequencing on the 454 platform (11); because the cost of sequencing on the Illumina platform is potentially considerably lower, we adapted hybrid capture using the NimbleGen 2.1M Human Exome Array to the Illumina DNA sequencing platform (see Methods). These arrays tile oligonucleotides from approximately 180,000 exons of 18,673 protein-coding genes and 551 micro-RNAs and comprise 34.0 Mb of genomic sequence. Resulting sequence data were processed using an automated pipeline: quality filtered sequence reads were aligned to the reference human genome (hg18) using Maq software (12). Single nucleotide variants (SNVs) were detected using Samtools (13) and further filtered (see Methods). For heterozygote calls we required at least 10× coverage by reads with different start sites, base call with Phred-like (14) quality scores greater than 45, and the probability of the frequency of major and minor alleles deviating from the binomial distribution of at least 10⁻⁷. Rare SNVs that cluster within 1 kb were tagged and evaluated for mapping errors. SNVs were annotated for effect on the encoded protein and for conservation by comparison versus sequences of 43 vertebrate species (15) and orthologues in fly and worm (see Methods).

Major changes in the protocol included shortening the length of the genomic DNA used for hybridization from 250 to 150 bases, which significantly increased the percentage of sequenced bases mapping to targeted bases (from 29% to 37%). We varied stringency by increasing the wash temperature from 47.5 °C to 53.5 °C. The higher stringency increased the percentage of on-target bases to 56%, but reduced breadth of coverage: at mean 30× per base coverage, the percentage of bases read fewer than 10 times increased from 7% to 23%.

In our current protocol, we captured and sequenced samples from 5 Caucasian subjects to a mean depth of 44×; these samples were also genotyped on the Illumina 370K SNP platform. The targeted bases constituted approximately 37% of all bases read (Table 1). An additional 12% of bases were within 100 base pairs of targeted sequences. The distribution of reads across all targeted bases is shown in Fig. 1A, which demonstrates that nearly all targeted bases are captured. Poorly captured bases, however, are correlated across different capture experiments, with approximately 352 kb (1.0%) read fewer than 10 times in each capture experiment. These underrepresented bases are predominantly from segments with exceptionally high G-C composition [supporting information (SI) Fig. S1]. We have captured 10 additional samples with this protocol with similar enrichment. With a single lane of Illumina sequencing using 75 base paired-end sequence using Illumina pipeline v. 1.4, we were able to produce an approximate mean 30× per-base coverage of the targeted bases, with 93% of targeted bases read at least 10 times; at 60× coverage with 2 lanes, 97% were read at least 10 times.

To assess the effect of coverage depth on sensitivity to detect sequence variants, one of these samples, GIT 264–1, was further sequenced to a mean depth of 99× and the genotype calls at known SNP positions were compared versus the results obtained with genotype calls on the Illumina 370K chip. At this depth of mean coverage, the sensitivity of exome sequencing for detecting homozygous and heterozygous variants was 99.8% and 98.2%, respectively (Fig. 1C). Sensitivity increased steeply as coverage increased from 5× to 20×, then more gradually thereafter, reaching a final plateau at approximately 50×. At any mean depth of coverage, there is a range of coverage at each base, and the ability to call heterozygous bases varies with the exact read coverage. Thus, the sensitivity to detect heterozygous variants with 10 reads is 78.6%, but increases to 95.2% at 20× and approximately 100% at 30× and greater (Fig. 1D). The overall per-base error rate for this sample was 0.75%, with error rate increasing approximately 10-fold from base 20 to base 75 (Fig. 1B).

Among the 5 test samples, we found a mean of 20,356 exomic variations from the reference sequence per subject (Table 1). An average of 1,055 per exome are not previously reported to our knowledge [i.e., not in the reference genome or 4 individual published genomes (16–19)] and were classified as novel. The small fraction of novel variants is consistent with a very high specificity for variant detection. These variants included a mean per exome of 12,372 coding SNVs (cSNVs; 642 novel), of which there were 5,341 missense variants (368 novel), 32 premature termination codons (7 novel), 81 canonical splice site variants (7 novel), and 6,815 synonymous substitutions (240 novel). Among the novel missense variants, 147 per exome introduce substitutions that are highly conserved among vertebrates, and among these, 40 are at positions that are also conserved in Drosophila melanogaster and/or Caenorhabditis elegans. Of the previously unreported exomic variants, 166 were shared among 2 or more subjects.

In subjects such as GIT 264–1 who are the product of consanguineous union (as described later), large segments of the genome are homozygous by descent; in these segments, all variations from the major allele can be considered errors, and heterozygous calls in such segments provide a further estimate of specificity (a small fraction of these could be bona fide de novo mutations). In this subject, 462 Mb of the genome were homozygous by descent as determined by genotyping (defined as segments of ≥100 consecutive SNPs spanning ≥1 Mb with ≥98% homozygosity; Table S1), including 5.3 Mb of the exome, comprising 2,459 genes. In these homozygous segments there was 1 heterozygous call at 99× coverage. This leads to an estimated false heterozygote discovery rate of 6 per exome and specificity for heterozygous calls of 99.9%. The specificity remains 99.9% at mean coverage of 30×. Interestingly, the false-positive call is present in the SNP database, so if one considers only novel variants in these calculations, specificity would be 100% at both levels of coverage.

These findings demonstrate a protocol of sufficient sensitivity and specificity to be highly useful for detecting rare sequence variants across the whole exome.

Subject GIT 264–1, a Turkish male, presented at age 5 months for evaluation of failure to thrive and dehydration. The history was notable for premature birth at 30 weeks and parental consanguinity, with 2 spontaneous abortions and death of a premature sibling on day 4; the parents were healthy (Fig. 2). Blood pressure was 90/55 mmHg and laboratory evaluation was notable for serum K⁺ level of 2.8 mmol/L, HCO₃⁻ level of 36 mmol/L; plasma renin activity and serum aldosterone levels were markedly elevated (Table 2). Despite volume depletion, his diapers were noted to be wet. The differential diagnosis of the admitting physician and consultants was broad, including a renal defect such as Bartter syndrome; however, a neurologic defect or an infectious process was not excluded. A venous blood sample was obtained and genomic DNA was prepared.

From the history, it was inferred that the patient likely had a recessive genetic trait. Genome-wide SNP genotyping was performed as described earlier. To search for copy number variants, the SNP data were interrogated by using a likelihood ratio-based copy number variant detection algorithm (see Methods). Within the homozygous-by-descent segments, 20 homozygous deletions were identified. All these were previously reported in the Database of Genomic Variants (20); none of these alter protein coding sequences (Table S2). No homozygous duplications were identified. Interestingly, none of the known loci for Bartter syndrome were within segments of homozygosity.

Given the diagnostic uncertainty, this patient's genomic DNA was subjected to whole-exome sequencing as described earlier. A total of 143 million reads passed quality assessment and 95.0% aligned to the human reference sequence; 38.4% of the total bases mapped to the targeted bases with mean coverage of 99×. At this depth of coverage, 99.3% of the bases were read at least 5 times and 98.4% of the bases were read at least 10 times (Table S3).

We anticipated finding a homozygous disease-causing mutation within a segment that is homozygous by consanguineous descent. We consequently sought novel and rare mutations that alter the encoded protein within the 5.3M base pairs and 2,495 genes in these segments. We identified 2,405 homozygous variations, including 1,493 homozygous cSNVs, from the reference sequence in these segments. These cSNVs included 668 non-synonymous substitutions (29 of which were novel), 791 synonymous coding variants (16 novel), 12 canonical splice site variants (0 novel) and 19 coding region indels (0 novel; Table S4). There were 3 premature termination codons in these segments (0 novel). The remaining 931 variants were in introns, UTRs, or intergenic regions. All of the novel missense variants were confirmed by Sanger sequencing of PCR-amplified segments, confirming the high specificity for detection of rare variation.

Substitutions at positions that are completely conserved from invertebrates to humans are highly likely to disrupt normal protein function (8), and we ranked the novel missense variants by conservation scores to identify the most likely functional mutations using the phyloP conservation score (Table S4). Ten of the homozygous variants were at highly conserved positions (i.e., score ≥2).

Among the novel homozygous mutations at extremely conserved positions, one strongly stood out: a single-base substitution that introduced a missense variant in SLC26A3 (Fig. 3). This mutation was identified in 172 of 173 reads covering this base; reads of this base had divergent start and end points and were read from both strands; the remaining read of this base was neither WT nor this mutant (Fig. 3A). The identity and homozygosity of this mutation was confirmed by PCR amplification and direct Sanger sequencing (Fig. 3B). This mutation introduces a D652N substitution at a position that is completely conserved among all invertebrates and vertebrates studied (Fig. 3C). This residue is also identical in 8 of 9 human paralogues (Fig. 3D); the only exception, SLC26A11, has glutamate rather than aspartate at the corresponding position. This mutation is noteworthy because SLC26A3 encodes an epithelial Cl⁻/HCO₃⁻ exchanger comprising 764 aa (Fig. 3E); recessive loss of function mutations in this gene are known to cause congenital chloride-losing diarrhea [CLD; i.e., OMIM 214700 (21)] in humans and mice (22). The D652N mutation lies in a β-pleated sheet in the STAS domain (Sulfate Transporters and bacterial Anti-Sigma factor antagonists) of the protein (Fig. 3E), which is known to be required for both the activity and biosynthesis/stability of the transporter. Indeed, other mutations at conserved positions in the STAS domain have been shown to cause congenital chloride diarrhea (23). This mutation was absent among 190 control chromosomes. We identified no other compelling homozygous or heterozygous candidates for disease-causing mutations outside homozygous chromosome segments.

Our genetic findings strongly suggested that this patient had congenital chloride diarrhea; the clinical findings of volume depletion, hypokalemia, and metabolic alkalosis are all consistent with this diagnosis. This prompted clinical follow-up with the referring physician to determine whether the volume loss was from renal or gastrointestinal sources. Follow-up revealed that the patient indeed had diarrhea that occurred every 1 to 2 h, which had been persistent. Urinary electrolyte measurements with the patient volume replete did not reveal hypercalciuria or other features of Bartter syndrome (urinary Ca²⁺:Cre ratio, 0.0076 mg:mg), indicating a primary diagnosis of CLD.

These findings led us to consider whether additional subjects referred with a presumptive diagnosis of Bartter syndrome might instead have CLD caused by mutation in SLC26A3. This mistaken diagnosis is particularly likely at the acute presentation: after volume and electrolyte resuscitation, urinary electrolyte levels can be appropriately high (when in balance, net intake and output are equal) but mistakenly interpreted as inappropriate, and volume loss from the urinary tract versus the gastrointestinal tract may not be immediately discriminated in infants.

We consequently screened 39 subjects referred with a suspected diagnosis of Bartter syndrome in whom mutations in NKCC2, ROMK, ClC-Kb, and Barttin had not been identified. We found 5 additional patients who had homozygous mutations in SLC26A3 (Table 2 and Fig. 4). These mutations included 2 unrelated subjects with a G187X premature termination codon, one with a frameshift mutation at codon 454 resulting in premature termination at codon 458, one with the same D652N mutation found in the index case, and one with a Y520C mutation (Y520 is completely conserved among all vertebrates and invertebrates studied; Fig. 4F). All these mutations are considered novel except for the G187X mutation, which is a known founder mutation for congenital chloride diarrhea in Saudi Arabia (24). None of these mutations were found in sequencing of 190 control chromosomes from unrelated Caucasian subjects. All these patients had hypokalemia, metabolic alkalosis, and evidence of salt wasting. Follow-up with the referring physicians in 3 cases documented that volume depletion was caused by gastrointestinal losses from watery diarrhea and not renal losses; in 2 cases, high stool chloride level was documented (GIT 366–1, 112 mM; GIT 409–1, 148 mM; value should be >90 mM in CLD), a finding that is pathognomonic for congenital chloride diarrhea. In each of the 2 remaining cases, the patient was not available for further evaluation.

The study protocol was approved by the Yale Human Investigation Committee. Consent for participation was obtained in accordance with institutional review board standards. Patients were referred for studies of Bartter syndrome, and consanguineous kindred subjects in whom no mutation in genes causing Bartter syndrome had been identified were chosen for further analysis. Genomic DNA was prepared from venous blood by standard procedures.

Genomic DNA was captured on a NimbleGen 2.1M human exome array following the manufacturer's protocols (Roche/NimbleGen) with modifications at the W. M. Keck Facility at Yale University. DNA was sheared by sonication and adaptors were ligated to the resulting fragments. The adaptor-ligated templates were fractionated by agarose gel electrophoresis and fragments of the desired size were excised. Extracted DNA was amplified by ligation-mediated PCR, purified, and hybridized to the capture array at 42.0 °C using the manufacturer's buffer. The array was washed twice at 47.5 °C and 3 more times at room temperature using the manufacturer's buffers. Bound genomic DNA was eluted using 125 mM NaOH for 10 min at room temperature, purified, and amplified by ligation-mediated PCR. The resulting fragments were purified and subjected to DNA sequencing on the Illumina platform. Captured and non-captured amplified samples were subjected to quantitative PCR to measure the relative fold enrichment of the targeted sequence.

Captured libraries were sequenced on the Illumina genome analyzer as single-end 50-, 74-, and 75-bp reads, or 75-bp paired-end reads, following the manufacturer's protocols. Image analysis and base calling was performed by Illumina pipeline versions 1.3 and 1.4 with default parameters.

Sequence reads were mapped to the reference genome (hg18) using the Maq program (12). Reads outside the targeted sequences were discarded and statistics on coverage were collected from the remaining reads using perl scripts. For indel detection, BWA was used to allow gapped alignment to the reference genome (29). SAMtools was used to call targeted bases, and any base call that deviates from reference base was regarded as a potential variation (13). Additional filters were applied as described in the SI text. Identified variants were annotated based on novelty, impact on the encoded protein, conservation, and expression using an automated pipeline (SI Text).

Samples were genotyped on the Illumina HumanCNV370-Duo BeadChips. Sample processing and labeling were performed following the manufacturer's protocols. Plink v. 1.06 was used to identify homozygous intervals from subject GIT 264–1. A sliding window of 50 SNPs was used on the tag SNPs and included no more than one possible heterozygous genotype. Resulting intervals have to meet the limit of at least 100 SNPs and 1 Mb in size.

To predict possible deletion or duplication events in the genome, LogR and B allele frequency data were extracted from SNP array and normalized to the larger sample pool to remove sample-specific noise. Based on the empirical distribution of LogR and B allele frequency values from the validated deletion and duplication SNPs, for a given SNP, the likelihoods of being 0, 1, 2, and 3 or more copies were calculated, respectively. If more than one SNP supporting deletion or duplication arose consecutively and the likelihood ratio of the SNPs exceeded threshold, they were called a deletion or duplication.

The ExonPrimer script (http://ihg2.helmholtz-muenchen.de/ihg/ExonPrimer.html) was used to generate primers for amplification of SLC26A3 coding exons using as a template genomic DNA of patients or controls (Table S5). Products were analyzed via gel electrophoresis, and amplicons sequenced using the forward or reverse primers. Disease-causing mutations were confirmed by at least 2 independent sequences from different primers.

Full-length orthologous and paralogous protein sequences from vertebrate and invertebrate species and human paralogous protein sequences were extracted from GenBank. Orthologues were confirmed based on database identity of annotation. If an orthologue could not be identified, the closest paralogue (i.e., top “hit” of a BLAST search of the respective proteome) was studied. Protein sequences were aligned using the ClustalW algorithm. GenBank accession numbers are listed in the SI Text.