Exome sequencing of 18,994 ethnically diverse patients with suspected rare Mendelian disorders

Introduction

The prevalence of rare disorders is considerable, with an estimated ~10,000 disorders having been identified with an average global prevalence of 1 in 2500 people¹. A substantial portion of these disorders have genetic etiologies and are difficult to manage clinically. The financial burden imposed on the patients, their families, and society is significant as well². One of the most significant factors that can contribute to reducing the mortality and the financial burden for rare disorder patients is receiving an accurate molecular diagnosis in a timely manner as it could allow for more precise clinical management³. Nevertheless, even after undergoing extensive testing, many patients fail to receive a definitive diagnosis, leading to what is often referred to as a ‘diagnostic odyssey’. The advent of next-generation sequencing (NGS) two decades ago has ushered in a new era of rapid and accurate diagnosis of these disorders, as evidenced by large cohort analyses^4,5,6.

Exome sequencing (ES) has emerged as a promising approach for the diagnosis of rare Mendelian disorders, providing diagnostic rates of around 30–50%^4,5,6,7. This is higher than the diagnostic rate of approximately 10% achieved through conventional array-based methods, which is also considered a genomic test⁷. By analyzing protein-coding regions of the genome, ES has proven to be a more comprehensive and cost-effective test compared to other conventional tests, especially for genetically heterogeneous disorders^8,9. Consequently, this advantage has led to the increasing adoption of ES worldwide as a routine clinical test for diagnosing rare diseases. However, the pace of adoption varies by country. Here, we report on the diagnostic performance of ES in an ethnically diverse cohort of 18,994 patients suspected of rare Mendelian disorders from 50 countries referred to a single reference laboratory in South Korea between 2020 and 2022.

Results

Patient characteristics

A total of 18,599 single probands, 374 family trios, and 21 family duos with suspected rare genetic disorders were subject to ES (Table 1). Fifty-six percent of the patients were male, and the median age of all patients at the time of sample accession was 15.1 years (range 0–107.6). The onset of symptoms during the neonatal or infancy period accounted for 45.4% of the cases, while 30.2% had their symptoms first present or noticed during adulthood. The neurodevelopmental delay (NDD) group was the largest with 22.2% of the patients and neurological (15.9%) and cardiovascular (13.3%) disorders were the next two largest groups. Almost all patients with NDD, skeletal disorder, and dermatological disorder developed symptoms from neonatal to childhood, while half of those with cardiovascular and neurological disorders developed symptoms in adulthood (Supplementary Table 1).

Table 1 Demographic characteristics of the 18,994 patients

Full size table

Ethnicity distribution

The ethnicity list the physician could select from during the order was as follows: East Asian, South Asian, Latino/Admixed American, African/African American, Non-finnish European, Ashkenazi Jewish, Amish, Middle Eastern, and Others. The predicted genetic ancestry group list was the same except for lacking Amish and Middle Eastern. As described in the methods, the genetic ancestry group for each sample was estimated based on the average probabilities that their genetic variations originated from a specific genetic ancestry group. In 38.9% of the cases (4257/10,952), the predicted genetic ancestry differed from the ordered ethnicity. Notably, 76.2% (421/552) of the patients ordered as African and 62.3% (527/845) of the patients ordered as non-Finnish European were predicted as ‘unassigned’. Additionally, 50.3% (1750/3476) of the patients ordered as ‘unassigned’ were predicted to be either East Asian or South Asian. The predicted genetic ancestry group matched the geographical definition of each country well (Fig. 1) with 99.6% from South Korea (6234/6261), 97.6% from Pakistan (991/1015), and 82.9% from South Africa (189/228) being predicted as East Asian, South Asian, and African, respectively. Based on the prediction, 60% were Asian, with 51.4% and 8.6% predicted to be East Asian and South Asian, respectively (Supplementary Table 2). There were still 25.9% of the patients predicted as unassigned and these patients were mostly from Egypt (38.9%, 1915/4926), Israel (9.6%, 472/4926), and Turkey (7.3%, 361/4926) (Fig. 1, Supplementary Table 2).

Exome sequencing of 18,994 ethnically diverse patients with suspected rare Mendelian disorders — **Fig. 1: Ethnicities of patients referred from medical institutions located worldwide.**

Diagnostic rate

The overall diagnostic rate with positive and inconclusive results was 31.8% (6034/18,994, 95% CI: 31.1–32.4%) and 12.6% (2402/18,994, 95% CI: 12.2–13.1%) respectively. The average turnaround time was 32.9 ± 14.4 days from accessioning to reporting.

Patients with early-onset (up to childhood) symptoms were associated with higher diagnostic rate of 38.2% (4658/12,179, odds ratio: 2.45, 95% CI: 2.28–2.62) than those with later-onset (adolescence ~ adulthood) symptoms with 20.2% (1376/6815, odds ratio: 0.41, 95% CI: 0.38–0.44) (Fig. 2). Patients with regions of homozygosity (ROH) of more than 1.6%, a theoretical percentage for offsprings between fifth-degree relatives¹⁰, had a diagnostic rate of 45.6% (1593/3490, 95% CI: 44.0–47.3), significantly higher than that of 28.6% (4441/15,504, 95% CI: 27.9–29.4) in those with ROH lower than 1.6% (odds ratio: 2.09, 95% CI: 2.26–1.94). As expected, 81% of the diagnosed patients with >1.6% ROH were reported with homozygous Pathogenic (P) or Likely pathogenic (LP) variants. The inclusion of the samples from either or both parents resulted in a higher diagnostic yield (proband-only: 31.6%, duo/trio: 41.3%, odds ratio: 1.52, 95% CI: 1.25–1.87).

**Fig. 2: Factors affecting diagnostic yield.**

Lastly, the diagnostic rate was highly variable by the disease category, ranging from 3.3% (1/30, 95% CI: 0–9.8%) for psychiatric disorders to 50.4% (185/367, 95% CI: 45.3–55.5%) for dermatological disorders (Fig. 2). The NDD group with the largest number of patients had a diagnostic yield of 41.3% (1740/4208, 95% CI: 39.9–42.8%). The disease categories were correlated with early age-of-onset (Chi-square test p-value < 0.001) and the presence of ROH (Chi-square test p-value < 0.001), but not with inclusion of parents (Chi-square test p-value > 0.05).

Characteristics of reported variants, genes, and disorders

The majority of the patients were diagnosed with SNV/INDEL only (5358/6034, 88.8%), while the rest were diagnosed with CNV (641/6034, 10.6%), repeat expansions or LINE-1 insertions (Table 2, Supplementary Tables 3 and 4). Chromosomal abnormalities and other mosaic or complex abnormalities were also found (n = 20). Of the 690 recurrent variants, 23 were diagnosed at least ten times, with the highest being 63 times. At the time of reporting, 23.0% of the variants were found to be novel, never observed in public databases or literature.

Table 2 Characteristics of the variants in positive cases

Full size table

A total of 1435 genes for 1704 disorders were reported for patients with positive results with SNV/INDEL (Supplementary Tables 5 and 6). Of 1704 disorders reported, autosomal dominant (AD), autosomal recessive, X-linked, digenic recessive, and mitochondrial disorders were reported in 49.1%, 35.7%, 8.4%, 0.05%, and 0.6%, respectively. As expected, a higher AD to AR disease ratio was observed in the cardiovascular (2.4%) and tumor syndrome (10.9%) groups. Hypertrophic cardiomyopathies were the most frequently diagnosed, which was not surprising considering the high proportion of cardiovascular patients. Concordantly, variants in MYBPC3 and MYH7 genes were reported frequently along with neurofibromatosis (NF1), Moyamoya disease (RNF213), and Cardiomyopathy (TNNI3). Dual diagnoses were given to 120 patients, and 4 were given triple diagnoses. Secondary findings were reported in 2.3% of the patients (463/18,994) with 50% of the reported genes (26/50) being cardiovascular-related, followed by cancer (38%, 19/50).

Internally classified likely benign variants

A total of 67,966 variants were identified in our patient cohort at a higher frequency than in gnomAD v4.0. These were the candidate variants that could be classified likely benign solely based on their frequency. Of these, 2096 were listed in ClinVar as P/LP/VUS. A total of 507 variants (2 LP, 3 Conflicts of classifications of pathogenicity, 1 not provided, and 501 VUSs) were deemed suitable for reclassification to LB, as they were either found to be too common in our samples with inconsistent phenotypes or found in patients diagnosed with different variants with unrelated phenotypes (Fig. 3a). In total, 97% (494/507) of the variants were found in at least one patient with a diagnostic variant in a different gene. More variants (410/507) were found in AR disorder genes, and 55 and 42 variants were in AD and XL disorder genes, respectively. Most of the reclassified variants were found in samples predicted to be of Asian ethnicity (60.3%) that is underrepresented in gnomAD (Fig. 3b, Supplementary Table 7).

**Fig. 3: Reclassification of variants more commonly found in the 18,994 patients than in gnomAD.**

Discussion

Here, we report insights gained from exome sequencing 18,994 consecutive patients with suspected rare disorders. As the largest and the most ethnically diverse cohort study to date that examined the diagnostic outcomes using exome sequencing^{9,10,11,12,13}, the study replicates the effectiveness of ES in diagnosing patients with a wide spectrum of rare genetic disorders. The predominant disease categories in our cohort were neurodevelopmental disorders (N = 4208, 22.2%), neurological disorders (N = 3027, 15.9%), and cardiovascular disorders (N = 2531, 13.3%). The factors affecting the diagnostic yields were also comparable to similar large studies¹⁴.

The overall diagnostic rate of 31.8% is comparable to other published exome or genome sequencing studies^{11,12,13,14,15}. However, when considering that our cohort is mostly composed of singletons with only 2.1% being duo/trios, the diagnostic rate of 31.8% for singletons is higher than the 20~25% diagnostic rate reported for singletons in other studies^13,16,17. This is mostly because of the high proportion of patients from consanguineous families (20%), which is associated with higher diagnostic rate. Proband-only patients without ROH in our cohort had a diagnostic yield of 28.4% which is more comparable to the published data. Another factor contributing to the higher diagnostic rate may be that more patients in our cohorts could have received ES as a first-tier test than in other studies that could be speculated to have included more patients who received thorough genetic workup before receiving exome sequencing test^17,18. Prior genetic testing information was limited and we could not perform analysis on how the history of prior genetic tests influences the diagnostic rate. However, the reason prior genetic testing information being scarce could be because there was no prior test done. We also investigated if different ethnic group distribution could have contributed to the higher diagnostic rate because of the four studies with over 2000 patients presenting diverse phenotypes, European Caucasians were the predominant ethnicity^13,16,19 in three studies that showed the distribution of the patients’ ethnic backgrounds. However, we did not observe a significant difference in diagnostic rates across the genetic ancestry groups among our patients with neurodevelopmental disorders (Chi-square test p-value = 0.49), the largest disease category (n = 4208, 22.2%).

Sixty-eight percent of the patients were not diagnosed, receiving either inconclusive or negative results. Inconclusive results were given to 12.6% (2402/18,994), of which 2.4% were duo/trios, not significantly different from the overall proportion of duo/trios. Some of these would have been reported positive if found to be de novo or compound heterozygous. Although we cannot predict what portion, we observed that 15 out of 374 trio cases would have been inconclusive if parents were not sequenced simultaneously. 1.6% of the inconclusive cases (39/2402) were reported with genes that were discovered recently and not annotated in OMIM yet (Supplementary Table 8). The potential diagnostic value of these genes was later demonstrated to be significant, as the genes from 23.1% of these cases (9/39) were subsequently curated with a disease in OMIM within a period of nine months following the conclusion of our study. The projected diagnostic yield including these 9 additional cases was nominal at 31.8% (difference = 0.06%). We can expect that most, if not all 9, of these cases will become positive eventually.

With the samples being referred from 308 customers from 199 different medical institutions spread across 50 countries, the patient population was diverse in terms of disease categories, age groups, and ethnicity. An interesting finding was that disorders such as phenylketonuria, typically not found by exome study due to newborn screening (NBS), were diagnosed in some patients. This could have been attributable to the limited availability of NBS in certain countries or to physicians seeking exact molecular basis of phenylketonuria. We identified 35 patients with phenylketonuria. More than half (23/35) of the cases were ordered from Egypt and Romania.

The low concordance (61.1%) between the ordered ethnicity and predicted genetic ancestry group is not surprising as the ordered ethnicity is most likely self-reported or physician-determined and therefore error-prone. The prediction being limited to the genetic ancestry groups in gnomAD v2.1, which are also predicted results from PCA and random forest model, is probably contributing to the low concordance as well. The high percentage (25%) of patients being predicted as unassigned is most likely due to these patients coming from admixed populations and/or from countries underrepresented in gnomAD. Collaborations with Asian institutions enriched our dataset with over 11,000 Asian patients (60% of our dataset), exceeding Asian representation such as Asian, Asian British, and Chinese in UK biobank (2.3%, n = 10,978)²⁰. In gnomAD v4.0, East and South Asians make up 45,546 individuals. However, still, we identified novel variants or variants that were more common in our dataset deemed likely benign that were not present in gnomAD v4.0. Notably, 507 variants classified as P/LP/VUS in ClinVar could be downgraded to LB. Therefore, we expect to find more ethnicity-specific disease-causing variants and common variants that can be used to reclassify many VUS in the public databases to LP or LB.

This study has several limitations. It is possible that patients who did not have genetic conditions were included in the study, which could have negatively affected the overall diagnostic yield. Even though we cannot determine or predict what percentage of our patients did not have a genetic condition, one reason for ordering a genetic test is to rule out a genetic etiology. For these cases, when the genetic test results are negative, the physician will most likely determine non-genetic etiology for the patient. Therefore, inclusion of these cases are reasonable. Clinical details were scarce for the majority of the samples, as access to the patient charts was not available. Almost 50% of patients had only 1 (n = 5731) or 2 (n = 3351) symptoms submitted in the order. However, physicians would frequently contact the lab to provide additional clinical information when the results were negative or when genetic findings did not correlate well with clinical information and vice versa, the lab would contact the physician for more clinical details in case there are strong P/LP variants that do not fit the phenotype with suspected lack of clinical information. Therefore, the possibility of missing diagnosis because of lack of clinical information is expected to be low. Technically, variants outside the exome or in so-called ‘NGS-dead zones’ with high sequence homology are challenging to detect²¹. Mosaic variants with low variant allele fractions can be missed and a VUS can be neglected due to lack of knowledge. ES also has lower resolution for CNV compared to genome sequencing, and cannot detect copy-number neutral structural variants with breakpoints outside the exonic regions. Finally, mitochondrial genome analysis was not included in samples received earlier in the study (n = 4438). Therefore, it is possible that a diagnosis was missed in these samples. Also, since the mitochondrial genome was not specifically targeted and only off-target reads were used, coverage across the mitochondrial genome varied across samples and was in general lower (Mean: 110.1, SD: 64.2) than the autosome, which may have resulted in further missing the diagnostic variants, especially lower-level heteroplasmic variants. An appropriate disclaimer was provided in the report. However, all these limitations should be similar for any large cohort diagnostic ES studies^13,17,22.

The overall diagnostic rate and the factors influencing the diagnostic yield were comparable to similar large studies. A notable finding was the presence of numerous novel or rare variants that appeared more frequently in this dataset compared to public datasets, indicating that certain populations remain underrepresented in public datasets. Since allele frequency is a key factor for determining variant pathogenicity, it is essential to have a sufficiently large and well-matched control dataset that can provide an accurate estimation of variant allele frequency. As we accumulate more data, we expect to continue improving the diagnostic performance and contribute to the rare disease community by providing unique and valuable genomic data. Also, since it is possible that there are subsets of patients who have disease-causing variants in the same gene not yet associated with human disease, future work includes querying for these patients as a first step of novel gene discoveries.

Methods

Patient inclusion

The study includes 18,994 patients referred to a single reference laboratory between 2020 and 2022 from 199 medical institutions across 50 countries. No exclusions were made based on age of onset or disease categories. Patients or patients’ legal guardians consented and had options to be informed of medically actionable secondary findings as defined by the American College of Medical Genetics (ACMG) guidelines^23,24. The study was conducted in a diagnostic setting and as all the samples and data were de-identified throughout, Institutional Review Boards (IRB) approval was not required.

Exome sequencing and bioinformatics analysis

Whole blood, buccal swabs, or extracted genomic DNA samples were collected from each patient. Genomic DNA was extracted from blood and buccal swab samples using QIAamp DNA Blood Mini Kit (Qiagen) and AccuBuccal DNA Preparation Kit (Accugene), respectively. Exome sequencing was performed using IDT xGen Exome Research Panel v2 (Integrated DNA Technologies, Coralville, Iowa, USA) for exome capture and NovaSeq 6000 system (Illumina, San Diego, CA, USA) for sequencing as 150 bp paired-end read. The base call (BCL) sequence files were converted and demultiplexed to FASTQ files using bcl2fastq v2.20.0.422 (Illumina, San Diego, CA, USA). Sequence reads were aligned to the Genome Reference Consortium Human Build 37 (GRCh37) and Revised Cambridge Reference Sequence (rCRS) of the mitochondrial genome using BWA-mem 0.7.17 to generate BAM files²⁵. Aligned BAM files were sorted using samtools (v.1.9)²⁶. Potential PCR duplicates were marked using Picard (v.2.20.8) (http://broadinstitute.github.io/picard/). GATK (v3.8) was used for Base Quality Score Recalibration (BQSR) and haplotype calling^27,28. Multiple programs were used to detect different types of variants (Supplementary Fig. 1). Single nucleotide variants (SNV) and small insertions/deletions (INDEL) within the targeted regions +/-100 bp flanking regions were called using GATK (v.3.8)^27,28. Copy number variants were called based on the depth-of-coverage (DOC) information using CoNIFER v0.2.2²⁹ and 3bCNV, an internally developed tool. Aneuploidy was predicted based on the DOC information of each chromosome. Repeat expansions were called with ExpansionHunter v5.0.0 using RepeatCatalogs-v1.0.0³⁰. Mobile element insertions (MEI) were called using Mobile Element Locator Tool (MELT) v2.2.2³¹. AutoMap v1.2 was used for detecting regions of homozygosity (ROH) for all patients³². The mean DOC was 154.6x per exome with a minimum 98% of the targeted region covered at 20x DOC. Various quality control metrics were checked to fall within the acceptable ranges (Supplementary Table 9).

Variant analysis by EVIDENCE

Variants were annotated, filtered, and prioritized using EVIDENCE, an internally developed system that incorporates a daily updated database module, a variant classification module, and a symptom similarity scoring module³³. Daily updated databases include public databases, in-house variant databases, and manually curated literature databases (Supplementary Table 10). Once the variants were annotated by variant effect predictor (VEP)³⁴ using the latest databases, variants with allele frequency >5% in gnomAD v2.1 (https://gnomad.broadinstitute.org/) were removed except for those previously reported as pathogenic (P) or likely pathogenic (LP) at least once in ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/). Variants were then classified as P, LP, variant of uncertain significance (VUS), likely benign (LB), or benign (B) according to the customized ACMG and the American Molecular Pathology (AMP) guidelines^35,36,37. Finally, the symptom similarity score was calculated for each variant by measuring how well matched the patient’s symptoms are to the known symptoms associated with the disease the variant occurred in. This was done by modifying the method from Kohler et al.³⁸, using the depth of HPO (Human Phenotype Ontology) terms instead of the information content of terms described in the study. For repeat expansions, the variant was considered pathogenic if the repeat number was equal to or above the normal repeat number range reported in the literature and OMIM database. MEI were considered as potentially disease-causing if the insertion occurred within a coding exon or at a known disease-causing region outside the exon. The final list of variants to be manually reviewed by medical geneticists to select reportable variants consisted of variants predicted to change amino acids within all known disease genes. Synonymous or intronic variants predicted to alter splicing by spliceAI³⁹ score >0.2 were also included. Genes only reported in literature as potential disease genes but not yet registered in OMIM were also included. Variants were prioritized by how pathogenic the variants were and how similar the patient’s symptoms were to the symptoms caused by the gene-disease the variant is found in. All reportable variants were examined using the Integrative Genomics Viewer (IGV, version 2.9.4)⁴⁰ for their technical validity and confirmed by Sanger sequencing if the variant quality was low^41,42.

Variant reporting

A positive report consisted of one or more variants that were determined to be disease-causing: for an autosomal dominant disease or an X-linked disease, one heterozygous or hemizygous P/LP variant and for an autosomal recessive disease, one homozygous P/LP variant, two P/LP (potential) compound heterozygous variants in a known disease gene that would fit the phenotype. An inconclusive report included one heterozygous or hemizygous VUS in a known autosomal dominant or X-linked disease gene that would fit the phenotype well, or potential compound heterozygous variants with one P/LP variant and one VUS, homozygous VUS variant, or only one P/LP variant in a known autosomal recessive disease gene. If P/LP variants were found in a gene that had not been previously associated with a disease defined in OMIM but reported in the literature, they were reported as inconclusive. Finally, if no clinically significant variant was found, a negative report was generated.

Genetic ancestry group prediction

The genetic ancestry group of each sample was predicted using an internally developed genetic ancestry group prediction model based on the allele frequency information of each population represented in gnomAD v2.1 data. To avoid variant sampling bias, only variants with an allele frequency of 0.05~0.95, allele count higher than 100, and allele number higher than 2000 across all populations were used. Linkage disequilibrium (LD) pruning was not applied. The genetic ancestry group score for each sample was defined by the geometric mean of conditional probabilities that variants would have occurred from each ethnic group. A genetic ancestry group with the highest score was assigned to each sample.

$${genetic; ancestry; group; score}={prod }_{x=1}^{n}Pr {({V}_{x}left|Eright.)}^{frac{l}{n}}={prod }_{x=1}^{n}{Allele; frequency; in},{E}^{frac{l}{n}},$$

where E means each genetic ancestry group.

Selection of likely benign variants

Variants that appeared in our dataset too frequently to be causal for rare diseases were selected and classified as likely benign. These variants are absent in gnomAD 4.0 as heterozygous for autosomal dominant (AD) genes, homozygous for autosomal recessive (AR) genes, and hemizygous for X-linked genes. As different genes/diseases have different penetrance, expressivity, and age-of-onset, a different threshold of allele frequency (AF) was used to select variants in the internal database: the gnomAD allele frequency of the most common ClinVar P/LP variant of each gene with 9 or more ClinVar submissions was used as the threshold and only the variants with AF above this threshold were further investigated for the possibility of being benign. We manually curated these variants as follows: for AD genes, if the variant was found as homozygous in at least one patient with no related symptoms or found in at least one patient diagnosed with variant(s) in a different gene while not having symptoms related to the gene harboring the variant in question, the variant was classified as likely benign. For AR genes, if the homozygous variant was found in at least one patient diagnosed with variant(s) in a different gene while not having symptoms related to the gene harboring the variant in question, the variant was classified as likely benign.