Rapid Whole-Genome Sequencing for Genetic Disease Diagnosis in Neonatal Intensive Care Units

Rapid WGS

To assess our ability to recapitulate known results, we performed rapid WGS retrospectively on DNA samples from two infants with molecular diagnoses that had previously been identified by clinical testing. Then, to assess the potential diagnostic use of rapid WGS, we prospectively performed WGS in five undiagnosed newborns with clinical presentations that strongly suggested a genetic disorder as well as their siblings.

Automation of the five main components of WGS as well as bioinformatics-based gene-variant characterization and clinical interpretation, all in an integrated workflow, made possible ~50-hour time to differential molecular diagnosis of genetic disorders (Fig. 1). Specifically, sample preparation for WGS was shortened from 16 to 4.5 hours, while a physician simultaneously entered into SSAGA clinical terms that described the neonates’ illnesses (fig. S1). For each sample, rapid WGS [2 × 100 base pair (bp) reads, including on-board cluster generation and paired-end sequencing] was performed in a single run on the Illumina HiSeq 2500 and took ~26 hours. Base calling, genomic sequence alignment, and gene variant calling took ~15 hours. The HiSeq 2500 runs yielded 121 to 139 gigabases (GB) of aligned sequences (34-to 41-fold aligned genome coverage; Table 1). Eighty-eight to 91% of bases had >99.9% likelihood of being correct (quality score >30, using Illumina software equivalent to Phred) (31, 32). We detected 4.00 ± 0.20 million nucleotides that differed from the reference genome sequence (variants) (mean ± SD) in nine samples, one from each of nine infants (Table 1).

Analytical metrics

In three samples, genome variants identified by 50-hour WGS were compared with those identified by deep targeted sequencing of either exons and 20 intron-exon boundary nucleotides of a panel of 525 recessive disease genes [Children’s Mercy Hospital Diagnostic panel 1 (CMH-Dx1)] or the exome (Table 2). CMH-Dx1 comprised 8813 exonic and intronic targets, totaling 2.1 million nucleotides (table S1) (2, 33). The exome and CMH-Dx1 methods, which used Illumina TruSeq enrichment and HiSeq 2000 sequencing, took ~19 days. In contrast, rapid WGS did not use target enrichment, was performed with the HiSeq 2500 instrument, and took ~50 hours. Samples CMH064, UDT002, and UDT173 were sequenced using these three methods, and variants were detected with a single alignment method [the Genomic Short-read Nucleotide Alignment Program (GSNAP)] (34) and variant caller [the Genome Analysis Tool Kit (GATK)] (35). Rapid WGS detected ~96% of the variants identified by a target enrichment method and ~99.5% of the variants identified by both methods had identical genotypes (Table 2), indicating that rapid WGS is highly concordant with established clinical sequencing methods (33). In contrast, analysis of the rapid WGS data set from sample CMH064 with three different alignment and variant detection methods [GSNAP/GATK, the Illumina CASAVA alignment tool, and the Burrows-Wheeler Alignment (BWA) tool] revealed surprising differences between the variants detected. Only ~80% of the variants detected using GATK/GSNAP or BWA were also detected with CASAVA (Table 2 and table S2) (36–41). This suggests that additional studies will be needed to define optimal alignment methods for clinical sequencing. Nevertheless, there was good concordance between the genotypes of variants detected by rapid WGS (using the HiSeq 2500 and CASAVA) and targeted sequencing (using exome enrichment, the HiSeq 2000, and GATK/GSNAP)—99.5% (UDT002), 99.9% (UDT173), and 99.7% (CMH064) (Table 2)—further indicating that rapid WGS is highly concordant with an established genotyping method (33). In subsequent studies, the rapid WGS technique used CASAVA for alignment and variant detection.

Genomic variants were characterized with respect to functional consequence and zygosity with a new software pipeline [Rapid Understanding of Nucleotide variant Effect Software (RUNES), fig. S2] that analyzed each sample in 2.5 hours. Samples contained a mean of 4.00 ± 0.20 million (SD) genomic variants, of which a mean of 1.87 ± 0.09 million (SD) were associated with protein-encoding genes (Table 1). Less than 1% of these variants (mean, 10,848 ± 523 SD) were also of a functional class that could potentially be disease causative (Table 1) (25–27). Of these, ~14% (mean, 1530 ± 518 SD) had an allele frequency that was sufficiently low to be a candidate for being causative in an uncommon disease (<1% allele frequency in 836 individuals sequenced at CMH) (42). Last, of these, ~71% (mean, 1083 ± 240 SD) were also of a functional class that was likely to be disease causative [American College of Medical Genetics (ACMG) categories 1 to 3] (Table 1). This set of variants was evaluated for disease causality in each patient, with priority given to variants within the candidate genes that had been nominated by an individual patient presentation.

Retrospective analyses

Patient UDT002 was a male who presented at 13 months of age with hypotonia, developmental regression. Brain magnetic resonance imaging (MRI) showed diffuse white matter changes suggesting leukodystrophy. Three hundred fifty-two disease genes were nominated by one of the three clinical terms hypotonia, developmental regression, or leukodystrophy; 150 disease genes were nominated by two terms; and 9 disease genes were nominated by all three terms (table S3). Only 16 known pathogenic variants had allele frequencies in dbSNP and the CMH cumulative database that were consistent with uncommon disease mutations. Of these, only two variants mapped to the nine candidate genes; the variants were both compound heterozygous (verified by parental testing) substitution mutations in the gene that encodes the αsubunit of the lysosomal enzyme hexosaminidase A [HEXA Chr 15:72,641,417T>C (gene symbol, chromosome number, chromosome coordinate, reference nucleotide > variant nucleotide), c.986+3A>G (transcript coordinate, reference nucleotide, variant nucleotide), and Chr15:72,640,388C>T, c.1073+1G>A]. The c.986+3A>G alters a 5′ exon–flanking nucleotide and is a known mutation that causes Tay-Sachs disease (TSD), a debilitating lysosomal storage disorder [Online Mendelian Inheritance in Man (OMIM) number 272800]. The variant had not previously been observed in our database of 651 individuals or dbSNP, which is relevant because mutation databases are contaminated with some common polymorphisms, and these can be distinguished from true mutations on the basis of allele frequency (33). The c.1073+1G>A variant is a known TSD mutation that affects an exonic splice donor site (dbSNP rs76173977). The variant has been observed only once before in our database of 414 samples, which is consistent with an allele frequency of a causative mutation in an orphan genetic disease. Thus, the known diagnosis of TSD was confirmed in patient UDT002 by rapid WGS.

Patient UDT173 was a male who presented at 5 months of age with developmental regression, hypotonia, and seizures. Brain MRI showed dysmyelination, hair shaft analysis revealed pili torti (kinky hair), and serum copper and ceruloplasmin were low. On the basis of this clinical presentation, 276 disease genes matched one of these clinical terms and 3 matched three terms (table S4). There were no previously reported disease-causing variants in these 276 genes. However, five of the candidate genes contained either variants of a type that is expected to be disease-causing based on their predicted functional consequence or missense variants of unknown significance (VUS). One of these variants was in a gene that matched all three clinical terms and was a hemizygous substitution mutation in the gene that encodes the α poly-peptide of copper-transporting adenosine triphosphatase (ATP7A Chr X:77,271,307C>T, c.2555C>T, p.P852L), aberrant forms of which are known to cause Menkes disease, a copper-transport disorder. This variant—new to our database and dbSNP—specified a nonconservative substitution in an amino acid that was highly conserved across species and had deleterious SIFT (Sorts Intolerant From Tolerant substitutions), PolyPhen2 (Polymorphism Phenotyping), and BLOSUM (BLOcks SUbstitution Matrix) scores. The known diagnosis of Menkes disease (OMIM number 309400) was recapitulated. As a further assessment of the reliability of variant detection of rapid WGS, samples UDT002 and UDT173 were aligned to the reference genome with three different alignment methods. The causative variants were recovered with each method.

Prospective analyses

Mutations in 35 genes can cause generalized, erosive dermatitis of the type found in CMH064 (table S5). The severe phenotype, negative family history, and absence of consanguinity suggested dominant de novo or recessive inheritance. No known pathogenic mutations were identified in the candidate genes that had low allele frequencies in the CMH cumulative genome and exome sequence database and similar public databases. Average coverage of the genomic regions corresponding to the candidate genes was 38.9-fold, and 98.4% of candidate gene nucleotides had >16× high-quality coverage (sufficient to rule out a heterozygous variant; table S6). Five candidate genes had 100% nucleotides with >16-fold high-quality coverage and, thus, lacked a known pathogenic mutation in an exon or within 20 nucleotides of the intron-exon boundaries. Eighteen of the candidate genes had >99% nucleotides with >16-fold high-quality coverage, and 31 had >95% nucleotides with at least this level of coverage. Furthermore, while 26 of the candidate genes had pseudogenes, paralogs, and/or repeat segments (table S6) that could potentially result in misalignment and variant miscalls, only 0.03% of target nucleotides had poor alignment quality scores.

Among the 35 candidate genes nominated by the phenotype, two rare heterozygous VUS were detected in CMH064; however, dideoxy sequencing of both healthy parents excluded one, in the keratin 14 gene, as a de novo mutation. The exomes of both parents were subsequently sequenced, and variants were examined in the trio at length. Three likely de novo mutations with excellent sequence coverage were identified in disease-causing genes. Of these, one was a candidate gene for CMH064. It was an in-frame deletion of three nucleotides in GJB2, (NM_004004), which encodes the connexin 26 protein. The variant, c.85_87del, p.Phe29del, removes a highly conserved amino acid within the first transmembrane helix (43). Dideoxy sequencing confirmed it to be a de novo mutation. Dominant, de novo GJB2 mutations have been associated with severe neonatal lethal disorders of the skin, such as keratitis-ichthyosis-deafness syndrome (KIDS), that involve the suprabasilar layers of the epidermis (OMIM number 148210) (44). The phenotype of CMH064 was atypical for KIDS, and functional studies are in progress to determine causality definitively.

Diagnoses suggested by the presentation in CMH076 were mitochondrial disorders, organic acidemia, or pyruvate carboxylase deficiency. Together, 75 nuclear genes and the mitochondrial genome cause these diseases (table S7). A negative family history suggested recessive inheritance that resulted from compound heterozygous or hemi-zygous variants or a heterozygous de novo dominant variant. Rapid WGS excluded known pathogenic mutations in the candidate genes. One novel heterozygous VUS was found. However, de novo occurrence of this variant was ruled out by exome sequencing of his healthy parents. No homozygous or compound heterozygous VUS with suitably low allele frequencies were identified in the known disease genes. Potential novel candidates included 929 nuclear genes that encode mitochondrial proteins but have not yet been associated with a genetic disease (45). Only one of these had a homozygous or compound heterozygous VUS with an allele frequency in dbSNP and the CMH database that was sufficiently low to be a candidate for causality in an uncommon inherited disease. Deep exome sequencing of both parents excluded this variant and did not disclose any further potentially causal variants.

A total of 174 genes are known to cause epilepsy of the type found in CMH172 (table S8). A positive family history of neonatal epilepsy and evidence of shared parental ancestry strongly suggested recessive inheritance. No known disease-causing variants or homozygous/compound heterozygous VUS with low allele frequencies were identified in these genes, which largely excluded them as causative in this patient. A genome-wide search of homozygous, likely pathogenic VUS that were novel in the CMH database and dbSNP disclosed a frame-shifting insertion in the BRCA1-associated protein required for ATM activation–1 (BRAT1, Chr 7:2,583,573-2,583,574insATCTTCTC,c.453_454insATCTTCTC, p.Leu152IlefsX70). A literature search yielded a very recent study of BRAT1 mutations in two infants with lethal, multifocal seizures, hypertonia, microcephaly, apnea, and bradycardia (OMIM number 614498) (46). Dideoxy sequencing confirmed the variant to be homozygous in CMH172 and heterozygous in both parents.

Rapid WGS was performed simultaneously on proband CMH184 (male), affected sibling (brother) CMH185, and their healthy parents, CMH186 and CMH202. Twelve genes have been associated with the clinical features of the brothers (heterotaxy and congenital heart disease; table S9). Co-occurrence in two siblings strongly suggested recessive inheritance. No known disease-causing variants or homozygous/compound heterozygous VUS with low allele frequencies were identified in these genes. A genome-wide search of novel, homozygous/compound heterozygous, likely pathogenic VUS that were common to the affected brothers and heterozygous in their parents yielded two nonsynonymous variants in the B cell CLL/lymphoma 9–like gene (BCL9L, Chr 11:118,772,350G>A,c.2102G>A, p.Gly701Asp and Chr 11:118,774,140G>A, c.554C>T, p.Ala185Val). Evidence supporting the candidacy of BCL9L for heterotaxy and congenital heart disease is presented below.

Automating medicine

Rapid WGS was made possible by two innovations. First, a widely used WGS platform has been modified to generate up to 140 GB of sequence in less than 30 hours (HiSeq 2500): Sample preparation took 4.5 hours, and 2 × 100 bp genome sequencing took 25.5 hours (Fig. 1). The total “hands-on” time for technical staff was 5 hours. Modifications included a new flowcell design and faster imaging and chemistry. Previously, NGS has either lacked sufficient sequence quantity, quality, or read lengths for clinical use of WGS or been too slow for use in acute patient care. Rapid WGS generated ~40-fold aligned genome coverage. The sequence quality was very similar to that obtained with its predecessor (HiSeq 2000), as determined by quality scores and alignment rates (48). Genotypes of nucleotide variants were >99.5% concordant with those of very deeply sequenced, partial exomes (33). The accuracy of the latter has been extensively benchmarked and is >99.9% (33).

Second, we automated much of the onerous characterization of genome variation and facilitated interpretation by restricting and prioritizing variants with respect to allele frequency (42), likelihood of a functional consequence (25), and relevance to the prompting illness. Thus, rapid WGS, as described herein, was designed for prompt disease diagnosis rather than carrier testing or newborn screening. SSAGA mapped the clinical features in ill neonates and children to disease genes. Thereby, analysis was limited only to the parts of the genome relevant to an individual patient’s presentation, in accord with guidelines for genetic testing in children (25–28). This greatly decreased the number of variants to be interpreted. In particular, SSAGA caused most incidental (secondary) findings to be masked. In the setting of acute care in the NICU, secondary findings are anticipated to impede facile interpretation, reporting, and communication with physicians and patients greatly (9, 49, 50). SSAGA also assisted in test ordering, permitting a broad selection of genes to be nominated for testing based on entry of the patients’ clinical features with easy-to-use pull-down menus. The version used herein contains ~600 recessive and mitochondrial diseases and has a diagnostic sensitivity of 99.3% for those disorders. SSAGA is likely to be particularly useful in disorders that feature clinical or genetic heterogeneity or early manifestation of partial phenotypes because it maps features to a superset of genetic disorders. SSAGA needs to be expanded to encompass dominant disorders and to the full complement of genetic diseases that meet ACMG guidelines for testing rare disorders (such as having been reported in at least two unrelated families) (26). Although neonatal disease presentations are often incomplete, only one feature is needed to match a disease gene to a presentation. In cases for which SSAGA-delimited genome analysis was negative, such as CMH064 and CMH076, a comprehensive secondary analysis was performed with limitation of variants solely to those with acceptable allele frequencies (42) and likelihood of a functional consequence (25). Nevertheless, secondary analysis was relatively facile, yielding about 1000 variants per sample.

RUNES performed many laborious steps involved in variant characterization, annotation, and conversion to HGVS (Human Genome Variation Society) nomenclature in ~2 hours. RUNES unified these in an automated report that contained nearly all of the information desirable for variant interpretation, together with a cumulative variant allele frequency and a composite ACMG categorization of variant pathogenicity (fig. S2). ACMG categorization is a particularly useful standard for prioritization of the likelihood of variants being causal (26). In particular, more than 75% of coding variants were of ACMG category 4 (very unlikely to be pathogenic). Removal of such variants allowed rapid interpretation of high-likelihood pathogenic variants in relevant genes. The hands-on time for starting pipeline components and interpretation of known disease genes was, on average, less than 1 hour. Because genomic knowledge is currently limited to 1 to 2% of physicians (physician scientists, medical geneticists, and molecular pathologists), variant characterization, interpretation, and clinical guidance tools are greatly needed, as is training of medical geneticists and genetic counselors in their use.

Return of results

In blinded, retrospective analyses of two patients, rapid WGS correctly recapitulated known diagnoses. In child UDT002, two heterozygous, known mutations were identified in a gene that matched all clinical features. In male UDT173, a hemizygous (X-linked) VUS was identified in the single candidate gene matching all clinical features. The variant, a nonsynonymous nucleotide substitution, was predicted to be damaging. Rapid WGS also provided a definitive diagnosis in one of four infants enrolled prospectively. In CMH172, with refractory epilepsy, rapid WGS disclosed a novel, homozygous frame-shifting insertion in a single candidate gene (BRAT1). BRAT1 mutations were very recently reported in two unrelated Amish infants who suffered lethal, multifocal seizures (46). A molecular diagnosis was reached within 1 hour of WGS data inspection in CMH172, even though extant reference databases [Human Gene Mutation Database (HGMD) and OMIM] had not yet been updated with a BRAT1 disease association. The diagnosis was made clinically reportable by resequencing the patient and her parents. Had this diagnosis been obtained in real time, it may have expedited the decision to reduce or withdraw support. The latter decision was made in the absence of a molecular diagnosis after 5 weeks of ventilatory support, testing, and unsuccessful interventions to control seizures. Given high rates of NICU bed occupancy, accelerated diagnosis by rapid WGS has the potential to reduce the number of neonates who are turned away. The molecular diagnosis was also useful for genetic counseling of the infant’s parents to share the information with other family members at risk for carrying of this mutation. As suggested by recent guidelines (9), this case demonstrates the use of WGS for diagnostic testing when a genetic test for a specific gene of interest is not available.

In four of five affected individuals, prospective, rapid WGS provided a definitive or likely molecular diagnosis in ~50 hours. These cases demonstrated the use of WGS for diagnostic testing when a high degree of genetic heterogeneity exists, as suggested by recent guidelines (9). Confirmatory resequencing, which is necessary for return of results until rapid WGS is compliant with Clinical Laboratory Improvement Amendments (CLIA), took at least an additional 4 days. Until compliance has been established, we suggest preliminary verbal disclosure of molecular diagnoses to the neonatologist of record, followed by formal reporting upon performance of CLIA-conforming resequencing. Staged return of results of broad or complex screening tests, together with considered, expert interpretation and targeted quantification and confirmation, is likely to be acceptable in intensive care. Precedents for rapid return of interim, potentially actionable results include preliminary reporting of histopathology, radiographic, and imaging studies and interim antibiotic selection based on Gram stains pending culture and sensitivity results.

Disease gene sleuthing

Because at least 3700 monogenic disease genes remain to be identified (10), WGS will often rule out known molecular diagnoses and suggest novel candidate disease genes (23, 51). Indeed, in another prospectively enrolled family, WGS resulted in the identification of a novel candidate disease gene, providing a likely molecular diagnosis. The proband was the second affected child of healthy parents. Accurate genetic counseling regarding risk of recurrence had not been possible because the first affected child lacked a molecular diagnosis. We undertook rapid WGS of the quartet simultaneously, allowing us to further limit incidental variants by requiring recessive inheritance. Rapid WGS ruled out 14 genes known to be associated with visceral heterotaxy and congenital heart disease (HTX). Among genes that had not been associated with HTX, rapid WGS of the quartet narrowed the likely pathogenic variants to two in the BCL9L gene. BCL9L had not previously been associated with a human phenotype but is an excellent candidate gene for HTX based on its role in the Wingless (Wnt) signaling pathway, which controls numerous developmental processes, including early embryonic patterning, epithelial-mesenchymal interactions, and stem cell maintenance (51, 52).

Recently, the Wnt pathway was implicated in the left-right asymmetric development of vertebrate embryos, with a role in the regulation of ciliated organ formation and function (53–57). The key effector of Wnt signaling is β-catenin, which functions either to promote cell adhesion by linking cadherin to the actin cytoskeleton via α-catenin or to bind transcriptional coactivators in the nucleus to activate the expression of specific genes (58–60). The protein that controls the switch between these two processes is encoded by BCL9L (also known as BCL9-2) and serves as a docking protein to link β-catenin with other transcription coactivators. BCL9L and α-catenin share competitive overlapping binding sites on β-catenin; phosphorylation of β-catenin determines which pathway is activated. The p.Gly701Asp mutation found in our patients lies within the BCL9L nuclear localization signal, which is essential for β-catenin to perform transcriptional regulatory functions in the nucleus (61).

BCL9L is one of two human homologs of Drosophila legless (lgs), a segment polarity gene required for Wnt signaling during development. lgs-deficient flies die as pharate adults with Wnt-related defects, including absent legs, and antennae and occasional wing defects (62). Fly embryos lacking the maternal lgs contribution display a lethal segment polarity defect. BCL9L-deficient zebrafish exhibit patterning defects of the ventrolateral mesoderm, including severe defects of trunk and tail development (60). Furthermore, inhibition of zebrafish β-catenin results in defective organ laterality (54). Overexpression of constitutively active β-catenin in medaka fish causes cardiac laterality defects (63). β-Catenin–deficient mice have defective development of heart, intestine, liver, pancreas, and stomach, including inverted cell types in the esophagus and posteriorization of the gut (64). Down-regulation of Wnt signaling in mouse and zebrafish causes randomized organ laterality and randomized side-specific gene expression. These likely reflect aberrant Wnt activity on midline formation and function of Kupffer’s vesicle, a ciliated organ of asymmetry in the zebrafish embryo that initiates left-right development of the brain, heart, and gut (56, 65). The second human homolog of lgs, BCL9, has been implicated in complex congenital heart disease in humans, of the type found in our patients (66–68). BCL9 was originally identified in precursor B cell acute lymphoblastic leukemia with a t(1:14)(q21;q32) translocation (69), linking the Wnt pathway and certain B cell leukemias or lymphomas (62). Finally, it was recently demonstrated that the Wnt/β-catenin signaling pathway regulates the ciliogenic transcription factor foxj1a expression in zebrafish (57). Decreased Wnt signal leads to disruption of left-right patterning, shorter/fewer cilia, loss of ciliary motility, and decreased foxj1a expression. Foxj1a is a member of the forkhead gene family and regulates transcriptional control of production of motile cilia (70). On the basis of this collected evidence, the symbol HTX6 has been reserved for BCL9L-associated autosomal recessive visceral heterotaxy. Additional studies are in progress to show causality definitively. These findings support clinical WGS as being valuable for research in reverse-translation studies (bedside to bench) that reveal new genetically amenable disease models.

Addressing limitations

In one remaining prospective patient, rapid WGS failed to yield a potential or definitive molecular diagnosis. Currently, WGS cannot survey every nucleotide in the genome (71). At 50× aligned coverage of the genome, WGS genotyped at least 95% of the reference genome with greater than 99.95% accuracy, using methods very similar to those used in this study (72). It has been suggested that this level of completeness is applicable for analyzing personal genomes in a clinical setting (72). In particular, GC-rich first exons of genes tend to be underrepresented (33). More complete clinical use of WGS will require higher sequencing depth, multi-platform sequencing and/or alignment methodologies, complementation by exome sequencing, or all three (73). Combined alignments with two methods of sequencing identified ~9% more nucleotide variants than one alone. However, these additions raise the cost of WGS, increase the time to clinical interpretation, and shift the cost-benefit balance.

For genetic disease diagnosis, the genomic regions that harbor known or likely disease mutations—the Mendelianome (2, 33)—must be genotyped accurately. In addition to exons and exon-intron boundaries, the Mendelianome includes some regions in the vicinity of genes that have structural variations or rearrangements. NGS of genome regions that contain pseudogenes, paralogs (genes related by genomic duplication), or repetitive motifs can be problematic. CMH064 had fulminant EB. Most EB-associated genes encode large cytoskeletal proteins with regions of constrained amino acid usage, which equate with low nucleotide complexity. In addition, several EB-associated genes have closely related paralogs or pseudogenes. These features impede unambiguous alignment of short reads, which can complicate attribution of variants by NGS. This limitation can prevent definitive exclusion of candidate genes. For example, 4.5% of nucleotides in KRT14, an EB-associated gene, had <16-fold high-quality coverage and, thus, may have failed to disclose a heterozygous variant. In CMH064, however, this possibility was excluded by targeted sequencing of the regions of KRT14 known to contain mutations that cause EB.

Furthermore, WGS is not yet effective for clinical-grade detection of all mutation types. Copy number variations and large deletions require clinical validation of research methods (33). Long, simple sequence–repeat expansions and complex rearrangements are problematic. Nevertheless, with CLIA-type adherence to standard operational processes, the component of the Mendelianome for which WGS is effective is extremely reproducible (33). Thus, the specific diseases, genes, exons, and mutation classes that are qualified for analysis, interpretation, and clinical reporting with WGS can be precisely predicted. This is of critical importance for reporting of differential diagnoses in the genetic disease arena. Thus, although insufficient alone, rapid WGS may still be a cost-effective initial screening tool for differential diagnosis of EB. In our study, all EB-associated genes had >95% nucleotides with high-quality coverage sufficient to exclude heterozygous and homozygous nucleotide variants (>16-fold); 19 of these genes had >99% nucleotides with this coverage. Hence, for rigorous testing of all EB-associated genes and mutation types, additional studies remain necessary, such as immunohistochemistry, targeted sequencing of uncallable nucleotides, and cytogenetic studies. Of 531 disease genes examined, 52 had pseudogenes, paralogs, repetitive motifs, or mutation types that may complicate WGS for comprehensive mutation detection. The comprehensiveness of WGS will be enhanced by longer reads, improved alignment methods, and validated algorithms for detecting large or complex variants (2, 4).

Finally, in singleton (sporadic) cases, such as CMH064, family history is often unrevealing in distinguishing the pattern of inheritance. For example, inheritance of EB can be dominant or recessive. Of two plausible heterozygous VUS detected in candidate genes in CMH064, one was a de novo mutation in connexin 26, which is associated with KIDS, that can be fatal in neonates (43, 44). The phenotype of CMH064 was not typical for KIDS, and functional studies are in progress. For evaluation of dominantly inherited diseases, WGS requires that the parents be concomitantly tested either by rapid WGS, by exome sequencing, or by resequencing of candidate de novo variants.

Rapid WGS failed to yield a definitive molecular diagnosis for CMH076. No known mutations were found in 89 disease-associated nuclear genes or the mitochondrial genome. This was an important negative finding because a molecular diagnosis of several of these genes is “actionable.” That is, specific treatments are indicated (such as pyruvate carboxylase deficiency, thiamine responsive congenital acidosis, biotinidase deficiency, fructose 1,6-bisphosphatase deficiency, and coenzyme Q10 deficiency). Likewise, exclusion of actionable diagnoses can prevent empiric institution of inappropriate treatments. Exclusion of known genetic diseases from a differential diagnosis is also of psychosocial benefit to family members and assists in guiding physicians regarding additional testing. There were no VUS with suitable inheritance patterns, in CMH076 or in either of the healthy parents, in known disease genes or in the remaining 929 nuclear-encoded mitochondrial genes (45).

In contrast to the rapidly declining cost of WGS, the computational cost of genome analysis is largely governed by Moore’s law (74). Sequence alignment, variant calling, and genotyping took 16 hours. Extremely rapid WGS is of practical use in clinical guidance only when married to equally rapid, cost-effective, deployable, and facile interpretation and analysis (2, 4). We are continuing to improve the speed of sequence base calling, alignment, and variant calling. It is likely that this interval can be halved and that HiSeq 2500–based rapid WGS can be performed in fewer than 36 hours by the end of 2012. Clinical validation of rapid WGS, however, will take some time.

Consent

This study was approved by the Institutional Review Board of CMH. Informed written consent was obtained from adult subjects and parents of living children.

Case selection

CMH is a nonprofit children’s hospital with 314 beds, including 64 level 4 NICU beds. It provides 48% of neonatal intensive care in the Kansas City metropolitan region. In 2011, the NICU had 86% bed occupancy. Retrospective samples, UDT002 and UDT173, were selected from a validation set of 384 samples with known molecular diagnoses for one or more genetic diseases. Seven prospective samples were selected from families with probands that presented in infancy, among 143 individuals without molecular diagnoses who were enrolled between 22 November 2011 and 4 April 2012 for exome or genome sequencing.

Clinicopathological correlation and interpretation

The features of the patients’ diseases were mapped to likely candidate genes. In part, this was performed manually by a board-certified pediatrician and medical geneticist. In part, it was performed automatically by entry of terms describing the patients’ presentations into a new clinicopathological correlation tool, SSAGA (2). It was designed to enable physicians to delimit WGS analyses to genes of causal relevance to individual clinical presentations, in accord with published guidelines for genetic testing in children and with NGS (9, 25, 28). SSAGA has a menu of 227 clinical terms, which are arranged in nine categories (fig. S1). SNOMED CT (Systematized Nomenclature of Medicine—Clinical Terms) (29) map to 591 well-established recessive diseases with known causal genes (table S1). Phenotype-to-disease-to-gene mapping was informed by Gene Reviews (30), OMIM Clinical Synopsis (10), MitoCarta (45), and expert physician reviewers.

Upon entry of the features of an individual patient, SSAGA nominates the corresponding superset of relevant diseases and genes, rank ordered by number of matching terms (fig. S1). It also contains a free-form text box that allows physicians to enter findings for which no SNOMED term exists, clinical term qualifiers, relevant family history, and specific genes of interest. The diagnostic sensitivity of SSAGA improves with use, by manual updating of mappings in cases where nominations failed to include the causal gene. SSAGA is extensible to additional diseases, genes, and clinical terms. Interpretation of results was manual on the basis of ranking of variant reports yielded by RUNES on SSAGA-prioritized candidate genes, supplemented with expert gene nominations (fig. S2). In some pedigrees, the presumed pattern(s) of inheritance allowed additional variant ranking based on obligatory genotypes in affected and unaffected individuals. Aligned sequences containing variants of interest were inspected for veracity in pedigrees with the Integrative Genomics Viewer (32).

Genome and exome sequencing

Isolated genomic DNA was prepared for rapid WGS with a modification of the Illumina TruSeq sample preparation (Illumina). Briefly, 500 ng of DNA was sheared with a Covaris S2 Biodisruptor, end-repaired, A-tailed, and adaptor-ligated. Polymerase chain reaction (PCR) was omitted. Libraries were purified with SPRI beads (Beckman Coulter). Quantitation was carried out by real-time PCR. Libraries were denatured with 0.1 M NaOH and diluted to 2.8 pM in hybridization buffer.

Samples for rapid WGS were each loaded onto two flowcells, followed by sequencing on Illumina HiSeq 2500 instruments that were set to rapid run mode. Cluster generation, followed by 2 × 100 cycle sequencing reads, separated by paired-end turnaround, were performed automatically on the instrument.

Isolated genomic DNA was also prepared for Illumina TruSeq exome or custom gene panel sequencing with standard Illumina TruSeq protocols. Enrichment for the custom gene panel was performed twice by Illumina hybrid selection with 20,477 eighty-nucleotide probes for 8366 genomic regions, representing exons and 20 intron-exon boundary nucleotides. It encompassed 2,158,661 bp, 525 genes, and 591 recessive diseases (2, 33) (table S1). The probes were designed to target 350 nucleotide genomic targets, with an average density of 2.4 probes per target (range, 2 to 56). Custom gene panel–enriched samples were sequenced on HiSeq 2000 instruments with TruSeq v3 reagents to a depth of >3 GB of singleton 100-bp reads in samples UDT173 and UDT002, respectively; 32.9 and 38.3% of base pairs were on target defined with a 0-bp extension, representing 469- and 501-fold enrichment in samples UDT173 and UDT002, respectively. Exome-enriched samples were enriched twice with standard Illumina hybrid selection and were sequenced on HiSeq 2000 instruments with TruSeq v3 reagents to a depth of >8 GB of singleton 100-bp reads per sample.

Genome and exome sequencing were performed as research, not in a manner that complies with routine diagnostic tests as defined by the CLIA guidelines.

Sequence analysis

CASAVA 1.8.2 (Illumina) performed gapped ELAND alignment of HiSeq 2500 sequences to the reference nuclear and mitochondrial genome sequences [Hg19 and GRCH37 (NC_012920.1), respectively] as well as variant identification. HiSeq 2000 sequences were aligned to the reference nuclear and mitochondrial genome sequences with GSNAP, and variants were identified and genotyped with the GATK (2, 34–36). Sequence analysis used base-call files, FASTQ files that contain sequences and base-call quality scores, the compressed binary version of the Sequence Alignment/Map format (a representation of nucleotide sequence alignments), and Variant Call Format (a format for nucleotide variants). Nucleotide variants were annotated with RUNES (2), our variant characterization pipeline, which incorporated VEP (Variant Effect Predictor) (37), comparisons to NCBI dbSNP, known disease mutations from the HGMD (38), and additional in silico prediction of variant consequences with ENSEMBL and UCSC gene annotations (39, 40) (fig. S2). RUNES assigned each variant an ACMG pathogenicity category (25–27) and an allele frequency on the basis of 722 patients sequenced since October 2011. Rapid WGS in CMH064 and exome sequences of his parents were also analyzed by Clinical Sequence Miner (deCODE Genetics).

Patient 1

CMH064 was a male born at 33 weeks gestation with erosive dermatoses. He was delivered vaginally following induction for preeclampsia. Desquamation and erythroderma from the hairline to occiput were present at birth (Fig. 2A). Denuded, hyperpigmented, and partially scarred lesions were noted above the upper lip, over the mentum, and in place of eyebrows. He had a truncated foreskin. His nails were dystrophic and yellowed. There were no vesicles, pustules, blisters, or mucosal lesions. Family history was positive for psoriasis. His mother had a healthy daughter from a prior union; there was no history of fetal loss. His father was healthy.

Cultures and herpesvirus PCR were negative. He developed severe neutropenia by day 3. Skin sloughing worsened. Rigid bronchoscopy and intubation were necessary because of fibrinous oropharyngeal exudate.

Skin biopsy histology revealed acantholysis, loss of cohesion between keratinocytes, and empty lacunae (Fig. 3A). There was focal dermal infiltration with neutrophils and lymphocytes and complete sloughing of the epidermal layer with focal clefting at the suprabasal layer (Fig. 3B). Immunofluorescence staining was negative for IgA (immunoglobulin A), IgM, and IgG except for linear staining for C3. Additional skin immunofluorescence studies revealed slightly reduced plakoglobin and desmoplakin and normal laminin 332; collagen types 4, 7, and 17; and plakophilin-1. Electron microscopy confirmed the absence of dermoepidermal junction (DEJ) separation and showed focally widened spaces between keratinocytes and cell vacuolization from the DEJ to the stratum corneum. Hemidesmosomes were normal. Some keratinocytes had large solitary vacuoles, abnormal condensation of keratin filaments, and perinuclear pallor. Some desmosomes had ragged edges. There were no intracellular inclusions. Negative laboratory studies included karyotype, Ro, La, Smith, ribonucleoprotein, and Scl-70 autoantibodies. Igs were unremarkable apart from an elevated serum IgA.

Sloughing of the skin, mucosal surfaces, and cornea continued to worsen, and by day 30, his activity level had markedly decreased (Fig. 2B). His fingers were edematous and discolored and had retained only three nails. On day 39, he developed purulent drainage from facial lesions. Skin cultures were positive for Escherichia coli and Enterococcus faecalis, and blood cultures for E. coli. Antibiotics were administered. He was thrombocytopenic and anemic, necessitating numerous transfusions. On day 47, ultrasound revealed nonocclusive portal vein and left brachiocephalic vein thrombi. By day 54, he developed metabolic acidosis, bloody stools, and persistent tachycardia. Medical interventions were withdrawn and he died on day 54. At autopsy, suprabasal acantholysis was present in the skin and the esophageal mucosa. Dideoxy sequencing of candidate genes KRT5, DSP, JUP, TP63, and KRT14 exons 1, 4, and 6 (the regions harboring most KRT14 mutations) was negative.

Patient 2

CMH076 was a male born at term with lactic acidosis, cardiomyopathy, and corneal clouding. He was born to a primigravid mother whose pregnancy was notable for decreased movements at 35 weeks gestation. His mother and father were healthy. Variable decelerations in heart rate were noted on the day before delivery. Labor was complicated by prolonged rupture of membranes, and delivery was by vacuum extraction for meconium staining. Apgar scores were 2, 3, and 5 at 1, 5, and 10 min, respectively. He had poor respiratory effort and hypotonia and required intubation. Upon transfer to CMH on day 2, he had lactic acidosis (lactate, 12 mmol/dl), coagulopathy, and cloudy corneas. Multiple cultures were negative. Echocardiogram showed chamber enlargement, reduction in biventricular function, noncompaction cardiomyopathy, mild tricuspid insufficiency, and mild aortic insufficiency. Urine testing revealed normal amino acids and elevated 3-methyglutaconic acid, 3-methylglutaric acid, and 2-ethyl-3-hydroxy-propionic acid. Long-chain fatty acids, acyl-carnitines, lysosomal hydrolases, β-galactosidase, β-glucuronidase, sphingomyelinase, glucocerebrosidase, α-L-iduronidase, and α-glucosaminidase were normal. Pressors were required for hypotension, and acidosis increased. He was diagnosed with hypoxic ischemic encephalopathy. On day 3, lactate was 28.2 mmol/dl. On day 5, respiratory distress worsened, accompanied by bloody endotracheal secretions; arterial pH was 7.04 and lactate was 22.0 mmol/dl. Medical interventions were withdrawn at the family’s request, and he expired on day 5. Postmortem testing by array comparative genomic hybridization (aCGH) and sequencing for mitochondrial transfer RNAs and TAZ, associated with Barth syndrome, were normal.

Patient 3

CMH172 was a female with intractable epilepsy. She was delivered at 39 weeks gestation by Cesarean section after an uncomplicated pregnancy. No exposure in utero to drugs, alcohol, or medications was reported. Birth weight was normal, length was 46 cm (<3%), and head circumference was 33 cm (<3%). Amniotic fluid was meconium-stained. Apgar scores were 6, 7, and 8 at 1, 5, and 10 min, respectively. Family history was positive for a female cousin with profound intellectual disability and infrequent seizures, and two cousins by a consanguineous marriage who died at 6 and 8 weeks of age of intractable epilepsy; all were from the same small Mexican town as the proband. Seizures started 1 hour after delivery. Antibiotics were given empirically until cultures and cerebrospinal fluid (CSF) herpesvirus PCR returned negative. Seizures continued despite multiple antiepileptic medications. CSF (including glycine level and CSF/plasma ratio) and brain MRI were normal. Electroencephalogram (EEG) showed focal epileptiform and sharp wave activity. Blood ammonia, electrolytes, pH, and glucose were normal. Oral feeding was poor. She was intubated and required increasing respiratory support for low SaO₂ and bradycardia. Ophthalmologic examination and radiologic skeletal survey were normal. An echocardiogram revealed a patent foramen ovale, tricuspid regurgitation, and peripheral pulmonary stenosis. Her karyotype was normal. aCGH was not diagnostic, but multiple tracts of homozygosity suggested shared parental ancestry. A repeat brain MRI at age 3 weeks was normal. Upon transfer to CMH at 5 weeks of age, she was small but symmetric, with bitemporal narrowing, micrognathia, flat nasal bridge, upslanted palpebral fissures, uplifted ear lobes, redundant helices, and fifth finger clinodactyly. She had hypertonia, persistence of cortical thumbs, hyperreflexia, clonus, and facial twitching. B6 challenge improved her EEG transiently, followed by return of multifocal sharp waves. Serum amino acids and urine organic acids were normal. Recurrent seizures continued both clinically and by EEG. After lengthy discussion, the parents requested withdrawal of support.

Patient 4

CMH184 was a male with visceral heterotaxy and congenital heart disease (dextro-transposition of the great arteries, total anomalous pulmonary venous return with pulmonary veins connecting to the right atrium, a large ventricular septal defect, pulmonary valve and main pulmonary artery atresia, mildly hypoplastic branch pulmonary arteries, patent ductus arteriosus with ductal-dependant left-to-right flow, and large atrial septal defect with obligate right-to-left flow). There was situs inversus of the spleen, liver, and stomach, with the aorta on the right of the spine and inferior vena cava on the left. Family history was positive for a 6-year-old brother (CMH185) with the same findings (dextrocardia, ventricular inversion, double outlet right ventricle, pulmonary stenosis, small pulmonary arteries, interrupted inferior vena cava with azygous continuation, and situs inversus of the liver and spleen). His parents (mother, CMH186, and father, CMH202) and two other siblings (one male and one female) were healthy. Testing of ZIC3, associated with X-linked recessive heterotaxy 1, was normal. Patient 4 remains in the NICU and is undergoing cardiac surgery.