Critical assessment of variant prioritization methods for rare disease diagnosis within the Rare Genomes Project

Sequencing, variant calling, and analysis by the RGP team

Genomic data were obtained by sequencing DNA purified from blood. Sequencing was performed by the Broad Institute Genomics Platform on an Illumina sequencer to 30x depth on average. Raw sequence reads were reassembled against the GRCh38 reference genome. Variants were called with GATK version 4.1.8.0 (12) in the form of single nucleotide variants (SNVs) and small insertions/deletions (indels). All data were analyzed by expert RGP variant analysts using seqr, an open-source, web-based genomic analysis tool for family-based monogenic disease analysis (seqr.broadinstitute.org/) (13). Structural variants (SVs) were not included in this challenge, but have been analyzed by the RGP team independent from the CAGI challenge.

Challenge datasets

Two datasets were provided for the CAGI6-RGP challenge, a training set and a test set. For each, a joint variant call file (VCF) was provided to the CAGI6 organizers for use in the challenge. In addition to the genomic data, clinical phenotype descriptions from patient provided information and review of medical records by a genetic counselor or medical geneticist were provided in HPO nomenclature. The diversity of phenotypes represented the range of clinical presentations routinely seen in patients referred for genetic testing. The family structure and affected status of each sequenced individual were provided, identifying the proband, sibling, mother, and father, as applicable.

For training and contextual purposes, genome sequencing and HPO data from 35 solved RGP families were provided along with the causal variant(s) identified by the RGP team. Ancestry was not provided but was imputed for the probands using the principal component analysis and random forest model used for the Genome Aggregation Database (gnomAD) (3). Overall, the training set consisted of six proband-only families, three duos (proband and one biological parent), and 26 trios (proband and both biological parents). The majority of diagnoses were de novo (n=21), or recessive (n=8) and the majority of responsible variants had been reported in the ClinVar database as P and/or LP (14) (https://www.ncbi.nlm.nih.gov/clinvar/) at the time the challenge was announced (May 3, 2021) (Supplemental Table 1).

For test purposes, the RGP team selected 30 families for inclusion in the challenge. Fourteen were solved and 16 were unsolved. The solved families in the test set were selected more stringently than for the training set, according to the following criteria: i) the responsible gene has an established Mendelian disease-association as per the Online Mendelian Inheritance in Man database (OMIM, https://www.omim.org/) and/or published literature at the time the challenge was announced, ii) the responsible variant(s) must not have been reported as P/LP in the ClinVar database or listed in/reported as a disease mutation (DM) in the HGMD Professional database (15) (https://www.hgmd.cf.ac.uk/ac/index.php) at the time the challenge was announced (May 3, 2021), and iii) the variant(s) were classified as P, LP, or variant of uncertain significance (VUS) with evidence that is close to LP according to the ACMG/AMP guidelines (6). The causal variants in all 14 solved families had been discussed by the RGP multi-disciplinary team of physicians, genetic counselors, analysts, and molecular geneticists, and had been returned to the family via a local clinician following confirmation in a CLIA certified laboratory. The local clinicians concurred that the variants are diagnostic. The submission of these variants from RGP participants to ClinVar was intentionally delayed for the duration of the challenge. Supplemental Table 2 displays the answer key for the 30 families in the test set. Overall, the test set consisted of two proband-only families, three duos, 23 trios (proband and both biological parents), and two quads (proband, affected biological sibling, and both biological parents). All of the 16 unsolved families were trios, with the exception of one quad.

A summary of the core features of the families and diagnostic variants in the training and test sets is depicted in Figure 1.

• Download figure
• Open in new tab

Figure 1. CAGI6 RGP challenge overview of selected families.

Summary of the 35 training set families (all solved) and 30 test set families (14 solved, 16 unsolved). Imputed population ancestry, the amount of familial sequencing data provided (proband-only, duo, trio, or quad), diagnostic status, and mode of inheritance of the causal variant(s) is displayed by family. For all returnable diagnostic variants in the solved families in each set, the functional consequence according to the Variant Effect Predictor (VEP), ClinVar and HGMD reporting status at the time of announcement of the challenge (May 3, 2021), and ACMG/AMP classification are displayed by variant. NFE, Non-Finnish European; AFR, African/African American; AMR, Admixed American; ASJ, Ashkenazi Jewish; SAS, South Asian; AD, autosomal dominant; XLR, X-linked recessive; AR, autosomal recessive; P, pathogenic; LP, likely pathogenic; VUS, variant of uncertain significance; DM, disease mutation.

Challenge format

The challenge was publicly released on May 3, 2021, and open for submissions on June 8, 2021. The submission deadline was October 11, 2021. Participating teams were tasked to provide a genetic diagnosis to as many probands from the 30 families in the test set as possible by submitting predictions for each proband’s causal variant(s). The 14 solved families were included in the challenge to evaluate the performance of each model in prioritizing the causal variants (true positives). The unsolved families were included with the goal of the participating teams identifying novel, potentially causal, variants for further clinical and experimental assessment followed, where possible, by return to the families. The number of solved and unsolved families was not disclosed in the challenge description to allow the participating teams to perform the task in a manner that reflects analysis in the clinical setting. Teams were able to submit up to 100 variant predictions per proband, ranked by causal likelihood, from a maximum of six different models. The submission format, a tab-delimited text file, accepted both single (one variant per line) and proposed biallelic (two variants per line), compound heterozygous predictions. For each variant, teams were required to provide an EPCR value for the submitted variant(s) being causal on a scale of 0 to 1, with 1 indicating highest certainty. An example submission file and a validation script were provided. Predictors were informed that assessors will review how often these were the top variant(s) returned (e.g., in the top 5, 10, 20, 50, or 100 variants) but not further informed about details of the assessment metrics.

Assessment of model performance across solved families

Formatting errors in the submission files were corrected, and redundant, duplicate, and incomplete submissions were removed. Causal variant predictions for each solved proband were assessed by an independent assessor (author S.L.S). The assessor was blinded to the identity and methods of the participating teams throughout assessment, and the identities of the participating teams were only revealed once the analysis was completed. The following two numeric metrics were considered:

• (i) Mean rank points: The mean of a weighted point allocation system based on the rank position of the true positive causal variant(s) in the solved probands within the top five (100 points), top 10 (50 points), top 20 (25 points), top 50 (10 points), or top 100 (5 points) variant predictions per proband. Model performance was subsequently ranked by the mean points awarded per proband.

• (ii) F-max: The F-measure, a harmonic mean between the precision and recall for causal variant prediction in the solved probands, was calculated for all unique EPCR values for each model. The maximum F-measure (F-max) (16), corresponding EPCR threshold, and mean number of predictions submitted per proband at and above this EPCR threshold were subsequently defined for each model and model performance was ranked by the resultant F-max value.

For both numeric metrics, a bootstrapped standard error (SE) (17) was calculated over 1,000 bootstrapped samples from the probands of the 14 solved families in the test set only.

As described, the causal variants in the answer key had been formally classified as P, LP, or VUS leaning towards LP according to the ACMG/AMP guidelines; however, for the purpose of matching the teams’ predictions to the answer key, all variants were treated equivalently. In the case that a correct causal variant was submitted in combination with a second non-causal variant in a proposed biallelic, recessive prediction, the prediction was considered incorrect. For P27, a proband from a family where both the proband and the affected sibling had inherited two paternal variants in cis, where it is unknown if both or only one of the variants is required and both variants were considered equally likely to be causal by the RGP team (Supplemental Table 2), the highest-ranked variant prediction for either one of the two variants by the respective model was retained and the other was removed from the analysis.

Assessment of novel putative causal variants across solved and unsolved families

Following assessment of model performance, predictions from top performing models that i) deviated from the answer key in the solved probands and ii) were submitted for the unsolved probands, were critically evaluated in the rare disease genomics web-based analysis tool seqr (seqr.broadinstitute.org/) (13). Putative causal variants were discussed by the RGP team and, where possible, were pursued by: i) functional validation of two cases by RNA sequencing, ii) SV analysis in a separate call set generated by the GATK-SV pipeline (18) and manually reviewed in the Integrative Genome Viewer (IGV) (19) to search for a second biallelic, compound heterozygous variant in the case of recessive disease genes, and iii) submission to the Matchmaker Exchange (matchmakerexchange.org/) via seqr in the case of candidate novel disease-genes.

Ethical considerations

The challenge data were derived from patients with rare, suspected monogenic, conditions and their close biological relatives, and included families who are medically underserved (20). Identification of putative causal variants, i.e., causal with respect to the clinical phenotype under investigation, may, if confirmed, be important for tailoring clinical interventions and obtaining social services. We did not actively search for variants unrelated to the rare condition in the family but the consent allows us to optionally provide clinical confirmation of secondary findings if they are incidentally discovered. For the purpose of this challenge, participating teams were told that pathogenic variants unrelated to the proband’s phenotype, such as might be identified as secondary or incidental variants in this challenge (21), should not be returned. All RGP participants have a consent video or phone call with a trained research coordinator to review the study protocol which includes provisions for sharing de-identified data and provide signed informed consent (Mass General Brigham IRB protocol 2016P001422). An institutionally signed (Broad-Northeastern) data transfer agreement was executed. We applied a registered access model (22) where all CAGI6 challenge predictors were required to sign and adhere to the CAGI Data Use Agreement (genomeinterpretation.org/data-use-agreement.html) but institutional signatures were not required.

Summary of submissions

Sixteen teams participated in the challenge, submitting predictions from a total of 52 models (median three models per team, range 1-6). Five teams elected to remain anonymous in the reporting, including one team (Team 6) that discovered a bug in their code during assessment and subsequently withdrew from the challenge. Between 0-100 variant predictions (single or proposed biallelic) were submitted per proband (range 0-100, median 100, mean 65). EPCR values ranged from 0-1 (median 0.32, mean 0.38) (Supplemental Figure 1). Ninety percent of predictions were single variants and 10% were proposed biallelic, compound heterozygous variants, including entries from five teams that only predicted single variants. Over half (53%) of all variant predictions were in established OMIM disease-associated genes, including entries from six teams limiting analysis to a subset of genes already implicated in disease, though the gene sets differed between groups. Eighty-four percent of predictions were in the coding sequence or direct splice region, as defined by VEP (i.e., within 1-3 bases of the exon, 3-8 bases of the intron, or in the splice polypyrimidine tract). Concordance between models for the top five ranked predictions per proband across all 30 families in the test set ranged from 0-1 (mean 0.09, standard deviation [SD] 0.15) and was significant only between different models from the same team (Supplemental Figure 2).

Summary of numeric assessment of model performance and methodology

Overall, model performance was highly variable (Figure 2A). All causal variants in the answer key were predicted within the first five rank positions by at least one model (Table 1). Our selected numeric assessment metrics for each submitted model are displayed in Table 2 and are depicted in Figure 2B. One of the top performing models from Team 9 (Invitae Moon) was able to call 13 of the 14 causal variants within the top five rank positions, followed by Team 12 (Lichtarge) with 12, Team 11 (enGenome) and Team 14 (TCS) tied with 10, and Team 5 (Exomiser) with 9 (Table 2). Here, we provide a summary of the numeric assessment of model performance and methodology for the five top performing teams. More detailed methods descriptions are provided for in total 11 of the 16 participating teams in the Supplemental Material.

• Download figure
• Open in new tab

Figure 2. Results of assessment using the 14 solved families (true positives).

A. Number of true positive diagnoses (y-axis) Identified per model (x-axis) colored by the rank position of the causal variants in the 14 solved probands. Models are ordered by their performance according to the mean rank points metric (Table 2). Team names are provided except for teams that elected to remain anonymous. B. Results of the mean rank points and F-max value numeric assessment metrics by team and model. Model 1, the primary model, for each team is indicated by the grey fill. C, Performance of models, according to the mean rank points awarded, comparing families with proband-only or duo data (i.e., an incomplete trio/quad) versus trio or quad data (i.e., a complete trio/quad).

Reanalysis of solved families

Given the high performance of these models, we reanalyzed the solved families in which models ranked variants higher than the causal variants identified by the RGP team in the answer key, to determine if they may contribute to disease or represent a more likely causal diagnosis; however, no compelling variants were found. To illustrate this, a detailed review of the variants prioritized by one of the top performing teams, Team 9 (Invitae Moon) in four solved probands (P2, P6, P7, and P11) is provided in the Supplemental Material.

Review of “difficult to predict” diagnoses

In genomics-driven diagnostics, failure to recognize causal variants and falsely prioritizing non-causal variants are recognized complications (5,32). We therefore reanalyzed families in the answer key for which predictors consistently failed to prioritize the causal variant(s). Several of these are described below.

The most poorly predicted diagnosis was KCND2 (c.1207C>G, p.Pro403Ala, ENST00000331113) in P6, a patient presenting with infantile-onset bilateral sensorineural hearing impairment, blindness, retinal dystrophy, hypotonia, chorea, profound global developmental delay, intellectual disability, and dystonia. Across all models, the causal variant was never reported at rank position one, was ranked at position 2-5 by just three models, and was only listed by 26 of 52 models (50%) across all variant predictions. This heterozygous ACMG/AMP LP missense variant in KCND2 explains the patient’s phenotype (33), is predicted to be deleterious by in silico prediction (REVEL 0.84 - PP3 Moderate) (8,30), and is absent from large population databases (gnomAD and TOPMed) (3,34). However, only duo sequencing was available for this family, from the proband and unaffected father; therefore, the de novo status of the variant remains unconfirmed. This may hinder models in prioritizing the variant. Indeed, calculating the mean rank points metric separately for families with proband-only or duo data versus those with trio or quad data, demonstrates a significant improvement in model performance with trio or quad data (paired Student’s T-Test p-value 0.00086) (Figure 2C). KCND2 is also not yet reported in the OMIM database as Mendelian-disease associated (last accessed April 2023). Thereby, models limiting their assessment to variants within OMIM reported Mendelian-disease associated genes would fail to prioritize this causal variant. This highlights the importance of OMIM and similar databases to the medical genomics community and the need to be able to represent novel gene-disease associations more rapidly. One such option for laboratories reporting novel Mendelian gene-disease relationships is to deposit them in the Gene Curation Coalition (GenCC) Database (https://thegencc.org/) allowing more rapid dissemination of findings to the community as well as the aggregation of many public and private gene-disease databases (35).

The second most poorly predicted diagnosis was PI4KA in P5, a patient presenting with global developmental delay, poor coordination, hypotonia, and spasticity, with an MRI-brain demonstrating cerebral hypomyelination and a dysplastic corpus callosum. Across all models, the two causal variants in this recessive gene were found at position 1-5 in nine models and were only listed by 18 of 52 models (35%) across all submitted variants. The first variant is a P nonsense variant (c.1852C>T, p.Arg618Ter, ENST00000255882; ACMG/AMP criteria applied: PVS1, PM2, PP1, PP3, and PP4). The second is a LP missense variant (c.4990G>A, p.Asp1664Asn, ENST00000255882; ACMG/AMP criteria applied: PP1, PP3, PP4, PM1 Supporting, PM2, PM3). One plausible explanation for the low prediction rate is that this is the only proband in the test set with a recessive diagnosis, requiring the models to jointly prioritize a pair of heterozygous variants. Notably, only 11 of 16 teams included proposed biallelic variants in their submissions, though several cases of recessive inheritance were included in the training set. Moreover, as with the KCND2 family above, this family was only sequenced as a duo (proband and unaffected father), demonstrating paternal inheritance of the nonsense variant, and requiring the assumption that the missense variant is maternally inherited or de novo, on the maternal haplotype, to constitute a recessive diagnosis.

The third most poorly predicted diagnosis was a splice acceptor variant in CLTC (c.1534del, p.Val512LeufsTer11, ENST00000621829), a gene associated with intellectual disability in OMIM (MIM: 617854). The proband (P19) presented with global developmental delay, hearing impairment, severely reduced visual acuity, constipation, hyperbilirubinemia, pulmonary arterial hypertension, and intracranial hemorrhage. This variant was ranked at position 1-5 by 13 models and was only listed by 18 of 52 models (35%) across all submitted variants. This de novo heterozygous LP splice acceptor variant (ACMG/AMP criteria applied: PS2, PM2, PVS1 Moderate) is predicted to cause a frameshift leading to a premature stop codon 11 amino acids downstream (in exon 10 of 31) in a highly loss-of-function constrained gene and is absent from large population databases. Moreover, since the challenge, the CLTC variant has been reported as LP in ClinVar by an independent submitter in association with intellectual disability (ClinVar variation ID: 811442). This variant arises at an acceptor splice site in the gene, thereby outside of the protein-coding region, which may explain the lower level of detection.

Finally, the fourth most poorly predicted diagnosis was TUBB8 in P27 (c.1039A>G, p.Asn347Asp and c.1033C>T, p.Leu345Phe, ENST00000568584), a female proband sequenced as a quad with her affected female sibling and both unaffected parents. In this family, two causal variants in TUBB8 were identified, inherited in cis from the unaffected father. Carriage of the causal variants by the unaffected father is explained by sex-limited expression of the oocyte maturation defect disease phenotype in females (MIM: 616780). To not exclude these variants as causal, the model would need to take sex-limited expression into consideration or allow for incomplete penetrance. Notably, amongst the 11 models ranking at least one of these causal variants at position one was Exomiser model 4, a model specifically allowing for incomplete penetrance, and all six models from DITTO, that did not take familial segregation into account.

Patterns in model performance

Following assessment of model performance, the prediction assessor was unblinded to the identity and methods of the participating teams. The wide variability in methodology, spanning stepwise filtering approaches to machine learning and artificial intelligence, did not allow for a comprehensive analysis nor use of statistical tests. However, a qualitative review of the methods demonstrated a pattern of decreasing performance when one or more of the following features were not considered by the model: i) variant call quality; e.g., depth, allele balance, and genotype quality (inclusion of sequence artifacts into submissions), ii) variant allele frequency; e.g., rare in large scale population databases such as gnomAD and TOPMed, iii) variant deleteriousness prediction; e.g., use of in silico tools and/or training on reported variants in clinical databases such as ClinVar and HGMD, iv) familial segregation within the provided dataset and inheritance mode of the respective gene, and v) relevance of the putative causal variant(s) to the proband’s phenotype. In some cases, all or most features were considered, yet the model did not identify many diagnostic variants, presumably due to the specific methodology used, information sources, and thresholds selected.

Summary of variant predictions in unsolved probands

Through reanalysis of the 16 unsolved families, directed by the submitted variant predictions from the top 10 teams, two additional families (12.5%) received a genetic diagnosis. The first, by the detection of a de novo splice region variant in TCF4 (c.1228+3G>T, ENST00000398339), prioritized by eight models in total, submitted by Team 9 (Invitae Moon, model 1 at rank 1), Team 5 (Exomiser, model 1-2 at rank 1 and model 3 at rank 2), and Team 11 (enGenome, model 2 and 4 at rank 1, and model 1 and 3 at rank 2). The second, by the detection of compound heterozygous frameshift (c.706del, p.Arg236GlyfsTer8, ENST00000175506) and deep intronic (c.1137+200_1137+205del, ENST00000175506) variants in ASNS, submitted as a biallelic prediction by Team 11 (enGenome, model 1, 2, and 4 at rank 1, and model 3 at rank 2) only. Notably, four additional models from Team 9 (Invitae Moon, model 1, rank 7) and Team 2 (AIBI, model 1, 5, and 6 at rank 83-91) prioritized the ASNS frameshift variant only. In both probands, the variant(s) impact on the transcript were functionally validated by RNA sequencing and were returned to the families following confirmation in a CLIA certified laboratory (Supplemental Table 3).

In a further six unsolved families, variants in putative novel disease genes were prioritized (Supplemental Table 3). For four of the six, a submission had already been made by the RGP team to Matchmaker Exchange (TPPP in P9, KCNH8 in P14, KLHL13 in P15, and THAP12 in P18). For the remaining two, new submissions were made (MRPL54 in P25 and FRY in P26). To date, Matchmaker Exchange matches warranting further consideration of these candidate genes have not been received, however, functional studies are underway for some candidates through the GREGoR consortium (https://gregorconsortium.org/). Across the remaining unsolved families, no variants identified were deemed of comparably high interest by the RGP team to pursue by functional studies or submission to Matchmaker Exchange.

Overall, there was more limited concordance in the variant predictions submitted between the top performing models in the unsolved families, compared to the solved families (Figure 3); and the vast majority of prioritized variants in the unsolved families did not merit further evaluation after review.

• Download figure
• Open in new tab

Figure 3. Concordance in the variant predictions submitted by top five performing teams in the solved and unsolved families.

Venn diagrams demonstrating the overlap in the variant predictions submitted across all probands in the solved families (left) compared to the unsolved families (right) between top performing teams.

The variants that did not merit further review in the unsolved families mostly fell into one or more of the following categories: i) heterozygous variants in dominant disease genes (according to the reported mode of inheritance in OMIM) inherited from an unaffected biological parent, and where incomplete penetrance is not expected based on current understanding, ii) heterozygous variants in dominant disease genes present in large population databases at an allele frequency higher than consistent with the prevalence of disease, where incomplete penetrance is not expected, and iii) single heterozygous variants in recessive disease genes that are unable to constitute a diagnosis without a second biallelic, compound heterozygous variant. For families with a single recessive variant and at least partial phenotype overlap with the reported phenotype, an SV call set generated by the GATK-SV pipeline was analyzed and the gene was manually reviewed in IGV with the aim to identify an SV in trans. This analysis did not, however, result in the detection of any additional variants of interest.

To provide one example, a heterozygous, maternally inherited, missense variant in GRIK2 was prioritized at rank position one by Team 9 (Invitae Moon) in P15. The variant (c.1066G>A, p.Gly356Arg, ENST00000421544) is predicted to be deleterious by in silico predictions (REVEL 0.95 – PP3 Strong) (36) and is absent in large population databases. GRIK2 is associated with dominant neurodevelopmental delay, impaired language, and ataxia (MIM: 619580) and with recessive intellectual disability (MIM: 611092). The dominant form of disease results from de novo gain-of-function variants clustering in a specific domain of ionotropic glutamate receptors, proven to affect channel kinetics and function (37,38). As the GRIK2 variant prioritized by Team 9 is inherited from the unaffected mother and falls far outside of this functional domain, it is inconsistent with being the cause of dominant disease under the assumption of complete penetrance, whereby every individual who has the variant with show signs and symptoms of the disease. The recessive form of disease results from biallelic loss-of-function variants (39,40). As the proband is lacking a second biallelic variant, the variant can also be deprioritized as a cause of recessive disease.

Returnable diagnoses identified in two unsolved families

For P1, Team 11 (enGenome) prioritized compound heterozygous putative loss-of-function variants in ASNS at rank position 1-3 across four submitted models; a maternally inherited frameshift variant (c.706del, p.Arg236GlyfsTer8, ENST00000175506) and a paternally inherited deep intronic 6 base pair deletion (c.1137+200_1137+205del, ENST00000175506). ASNS is a disease gene associated with asparagine synthetase deficiency (MIM: 615574) and is a phenotype match for the proband, who presented with Lennox-Gastaut syndrome, infantile spasms, microcephaly, hypotonia, nystagmus, optic nerve hypoplasia, partial agenesis of the corpus callosum, and delayed myelination. Loss-of-function of ASNS is an established disease mechanism in autosomal recessive asparagine synthetase deficiency (41,42). The frameshift variant is rare in large population databases (absent in gnomAD, reported in 1/264,690 alleles in TOPMed) and has recently (Feb, 2022) been reported as P in ClinVar (ClinVar variation ID: 1411238). The variant leads to a premature stop codon in the middle of the gene, in exon six of 13, and is expected to result in a truncated protein. The variant is classified as LP according to ACMG/AMP guidelines (criteria applied: PVS1 and PM2 Supporting). The deep intronic indel between exons 10 and 11 (200bp away from the exon), is absent from large population databases and has a moderate SpliceAI score (0.2) (43) predicting acceptor gain. RNA sequencing analysis performed on blood from the proband demonstrated evidence of complex splice disruption, including intron retention and novel exon creation, leading to a premature stop codon in the middle of the gene (Figure 4A). In light of this evidence, the variant was classified as LP according to ACMG/AMP guidelines (criteria applied: PS3, PM3, and PM2 Supporting). ASNS was deemed a clinical fit by the family’s local physician. A cerebrospinal fluid (CSF) asparagine level was measured in the proband and was found to be within normal range. Though low CSF asparagine level would further support the diagnosis, normal levels have previously been reported in patients with ASNS defects, due to limitations in the sensitivity of the assay (41,44). The family is now pursuing oral asparagine therapy.

• Download figure
• Open in new tab

Figure 4.

Confirmatory RNA sequencing in P1 and P3. For both A and B, in the top panel, paired end reads from the RNA sequencing BAM file are displayed for the proband. In the lower panels, the RNA sequencing read pileup tract is displayed with the novel (orange) and known (blue) junctions annotated in the proband and in aggregated data from GTEx controls, respectively. Beneath, the gene transcript isoforms are displayed. A, RNA sequencing analysis performed on blood in P1 compared to normalized GTEx blood samples (n=755) (48). The results for ASNS (displaying exon 9 and 10) demonstrate evidence of splice disruption due to a deep intronic indel (indicated by the red box in the proband) with cryptic exon creation and intron 9 read-through. B, RNA sequencing analysis performed on an EBV-transformed lymphoblastoid cell line (LCL) in P3 compared to normalized GTEx lymphocyte samples (n=174). The results for TCF4 (displaying exon 10 to 13) demonstrate evidence of splice disruption due to a near-splice variant (indicated by the red line in the proband) with skipping of exon 11 in approximately 20% of reads. E, exon.

In P3, three top performing teams, Team 9 (Invitae Moon), Team 5 (Exomiser), and Team 11 (enGenome), prioritized a de novo variant in TCF4 (c.1228+3G>T, ENST00000398339), a disease gene associated with dominant Pitt-Hopkins syndrome (PHS, MIM 610954). This splice region variant has a moderate SpliceAI score (0.72) predicting donor loss and is absent from large population databases. Moreover, it is a putative loss-of-function variant in a highly loss-of-function constrained gene (pLI score 1, LOEUF 0.22, gnomAD) for which loss-of-function is an established disease mechanism (45). This TCF4 variant was flagged during analysis in seqr by the RGP team. However, at the time, it was considered non-compelling due to the absence of classical PHS features in the proband, such as dysmorphism, including a large beaked nose, wide mouth, fleshy lips, and clubbed fingertips, and abnormal breathing patterns, presenting as hyperventilation episodes. The phenotypic spectrum of TCF4 has, however, since been expanded to include neurodevelopmental delay in the absence of classical PHS (46). Moreover, upon re-contacting the family for additional clinical information and to request photographs, abnormal breathing patterns and mild dysmorphic features supporting PHS were confirmed. The variant has recently (Aug, 2021) been independently reported in ClinVar as LP (variation ID: 1204043), and has been reported in a study generating patient-specific induced pluripotent stem cells to model PHS (47). RNA sequencing analysis performed on cultured lymphoblasts from the proband demonstrated evidence of splice disruption with exon skipping in the middle of the gene, in exon 11 of 20 (Figure 4B). The variant was thereby classified as LP according to ACMG/AMP guidelines (criteria applied: PVS1 and PM2 Supporting).