Pathogenic neurofibromatosis type 1 (NF1) RNA splicing resolved by targeted RNAseq

Neurofibromatosis type 1 (NF1) is caused by loss-of-function variants in the NF1 gene. Approximately 10% of these variants affect RNA splicing and are either missed by conventional DNA diagnostics or are misinterpreted by in silico splicing predictions. Therefore, a targeted RNAseq-based approach was designed to detect pathogenic RNA splicing and associated pathogenic DNA variants. For this method RNA was extracted from lymphocytes, followed by targeted RNAseq. Next, an in-house developed tool (QURNAs) was used to calculate the enrichment score (ERS) for each splicing event. This method was thoroughly tested using two different patient cohorts with known pathogenic splice-variants in NF1. In both cohorts all 56 normal reference transcript exon splice junctions, 24 previously described and 45 novel non-reference splicing events were detected. Additionally, all expected pathogenic splice-variants were detected. Eleven patients with NF1 symptoms were subsequently tested, three of which have a known NF1 DNA variant with a putative effect on RNA splicing. This effect could be confirmed for all 3. The other eight patients were previously without any molecular confirmation of their NF1-diagnosis. A deep-intronic pathogenic splice variant could now be identified for two of them (25%). These results suggest that targeted RNAseq can be successfully used to detect pathogenic RNA splicing variants in NF1.


INTRODUCTION
Neurofibromatosis type 1 (NF1; MIM 162200) is one of the most common autosomal dominant tumor-predisposition disorders affecting~1 in 3000 individuals 1,2 . Patients typically present with the following clinical manifestations: café-au-lait spots, dermal neurofibromas, iris hamartomas (Lisch nodules), axillary and/or inguinal freckling, and subcutaneous or plexiform neurofibromas. Other manifestations include learning disabilities, optic glioma, skeletal abnormalities and an increased risk of specific malignancies [1][2][3] . Malignancies are an important cause of the 8-15-year reduction in average life expectancy [1][2][3] . Most patients are clinically diagnosed in childhood, according to NIH consensus criteria 4 . However, genetic testing to confirm a clinical diagnosis is still warranted, because of the clinical overlap with Legius syndrome (MIM 611431; SPRED1) 1,2,5 . In addition, an accurate genetic diagnosis facilitates appropriate screening and follow-up, reproductive options (prenatal testing and pre-implantation genetic diagnosis) and access to clinical studies and potential future trials and treatments. NF1 is caused by loss-of-function variants in the tumorsuppressor gene NF1 (MIM 613113), located on chromosome 17q11.2 and consisting of~350 kb of genomic DNA. The NF1 gene produces a major 12 kb transcript NM_000267.3 that contains 57 exons and encodes for a Ras-guanosine triphosphatase (GTPase) activating protein: neurofibromin 1,6-8 . Pathogenic variants in NF1 result in loss of function of neurofibromin causing an increase in Ras signaling, affecting cell proliferation and differentiation 9 .
The diagnostics laboratory community has adopted the ACMG standards and guidelines for the interpretation of sequence variants 18 to determine their clinical relevance, including recommendations for splice variants. For variants adjacent to exon boundaries and not further than two nucleotides from the exon border, in silico tools often adequately predict changes in RNA splicing 19 . In genes were loss-of-function is a known mechanism of disease, a null variant such as a variant within the canonical splice site could get a pathogenic very strong 1 (PSV1) score. (See Richards et al 18 for the full explanation). However splice prediction tools for variants outside the consensus splice acceptor/donor sites perform poorly 20 . And although these tools are able to predict cryptic or novel splice sites, the in vivo outcome of competing splice sites for the splicing complex is difficult to predict 21 . Moreover, the specificity of predictions for the effect of genetic variants on splice enhancer or silencer sites, potential branch points 22 or even change in the branch point-exon distance 23 is limited. As a consequence, the effect of DNA variants on RNA splicing ultimately needs to be determined experimentally. A variant that shows an effect with a well-established functional study; especially a robust, reproducible and validated study in a clinical diagnostic laboratory will get a pathogenic strong 3 (PS3) score. As such a potential splice variant that shows a loss-of-function effect could get a PS3 score and in some cases a PVS1 score.
The functional study most often employed for splice variant validation is conventional Reverse Transcription (RT)-PCR, which has its drawbacks. First, it is laborious, time consuming and prone to preferential amplification of transcripts depending on primer choice. Another possible pitfall is the use of primers that amplify only a small part of the region surrounding the variant, while the genetic effect may involve skipping of multiple exons 21 . Furthermore, data are only semi-quantitative and, finally, almost every new splice isoform needs its own RT-PCR primer design. Targeted RNAseq on the other hand gives comprehensive information about RNA expression and splicing for genes/transcripts and can circumvent the aforementioned draw-backs 24 .
A diagnostic procedure was designed that uses an in-house developed tool QURNAs (Quantitative enrichment of aberrant splicing events in targeted RNAseq) to facilitate the detection and quantification of normal and/or pathogenic RNA splicing events in targeted RNA sequencing data.

RESULTS
QURNAs detects normal and known pathogenic NF1 splicing events To validate our diagnostic procedure, targeted RNAseq was performed using blood samples of nine patients with a known pathogenic NF1 splice-variant, and two control samples. (Table 1).
First, the presence of all normal exon-exon splice junctions belonging to NF1 reference transcript NM_000267 was checked and found detectable with a median enrichment score (ERS)~1 and >5000 reads per splice junction across all samples (Fig. 1a, Supplementary Data 1 and 2). Furthermore, with QURNAs an extensive list of (naturally occurring) splicing events for NF1 was obtained (Supplementary Data 2) and compared to previously published naturally occurring events (Supplementary Data 3). From the 25 previously reported normal (rare) events, 24 were detected. Approximately half of these previously reported events are major splicing events, detected by QURNAs in almost all of the patient derived cells (median ERS~1 & reads >1000). The other half are minor events (median ERS < 1 & reads <700) and detected in a fraction of the patient derived cells (Supplementary Data 2). The only missed previously reported normal rare event is from a transcript 'ins 9a/9br' (legacy exon nomenclature, corresponding to insertion of 10 amino acids between exons 11-12 of NF1; c.1260 + 1617_1260 + 1646) that is supposed to be specific for neuronal tissue, and not detectable in blood [25][26][27] . Additionally 45 events not published before were detected in multiple samples (Supplementary Data 2).
For the nine samples with a known pathogenic splice-variant all expected (major) pathogenic splicing events were detected using QURNAs. The corresponding causal DNA variant could also be observed in the RNAseq reads (except for one -data not shown). The highest ERS per sample and concomitant major drop in ERS of the normal exon-exon splice-junctions due to alternative splicing are listed ( and 11 had additional (minor) splicing events that are related to their major pathogenic splicing event albeit with lower ERS. Events with a maximum ERS of~7 were found for the two control samples. However, these events had relatively low reads per event (< 150) and do not appear to be true splicing events based on in silico predictions ( Supplementary Fig. 2). Also, the concomitant drop in ERS of the normal exon-exon splice junctions in close proximity was absent, indicating that wild-type splicing was not significantly altered (Fig. 1a). This indicates that different data sources have to be taken into account for correct interpretation of an increased ERS, to be able to discriminate artefacts/non-splicing related events from true alternative splicing events. This is further illustrated by samples 6 and 7. For sample 6 of the validation cohort with variant c.586 + 5 G > A, previously shown to cause skipping of exon 5 28-30 , QURNAs indicated an ERS of 6.1 for exon 5 skipping (covered by 6509 reads). This is also reflected in Fig. 1a by a drop in ERS at the exon 4-5 junction (ERS~0.61) and 5-6 junction (ERS~0.63) compared to the other samples in the same run (median ERS~1 for exon junction 4-5 and median ERS~1 for 5-6 across all samples-see Fig. 1a, Supplementary Fig. 3). The ERS of 6.1 for exon 5 skipping is smaller than expected, most likely because exon 5 skipping is a (previously published) naturally occurring event (Δ5 or 'NF1-ΔE4b') 31 that is detected with a median ERS of 0.52 and median 728 reads across all samples ( Supplementary Fig. 3). Sample 7 of the validation cohort harbors variant c.1466 A > G (p.(Tyr489Cys)) in NF1. This specific variant has been reported many times as a disease-causing splice variant, since RNA studies have shown that it creates a new donor site leading to the deletion of the last 62nt of exon 13 ( Supplementary Fig. 4) 15,28,32 . This is confirmed using RNAseq (high ERS of 22.4 and covered by 1792 reads). The event corresponds to a drop in ERS to 0.59 for sample 7 at the junction of exon 13-14 compared to median ERS~1 across all samples for this event in the same run. The observed drop in ERS is due to partial loss of normal exon 13 to exon 14 splicing. However, one event with a higher ERS score was also observed for this sample (ERS 23.4; 228 reads - Table 1, Supplementary Data 1). Again, this event corresponds to an artefact based on in silico predictions ( Supplementary Fig. 4), the relatively low read number and the concomitant drop in ERS of the normal exon-exon splice junctions in close proximity is lacking (Fig. 1a).
RNAseq replication of RT-PCR based molecular diagnosis Ten RNA samples from an external diagnostic laboratory that uses RT-PCR to detect pathogenic NF1 splicing, were subsequently retested using the targeted RNAseq assay and processed with QURNAs. Initially, the RT-PCR results were blinded to us to test the performance of the targeted RNAseq assay in a simulated diagnostic setting. The RT-PCR results could be fully replicated with some minor differences and could even be complemented with splicing events that are more difficult to detect or extract from conventional RT-PCR and Sanger sequencing data (see Table  2; further discussed below). Using the targeted RNAseq diagnostic procedure all pathogenic splicing events were detected with high ERS for the aberrant event (Table 2) and the concomitant drop in ERS of the normal exon-exon splice junctions in close proximity (Fig. 1b). For these samples, as for the validation set, we could detect all normal reference exon splice junctions (median ERS~1, median reads > 800) and 24 published naturally occurring splicing events that are detectable in blood derived RNA (Supplementary Data 2, 4). Besides detection of the causal splicing events using QURNAs, the corresponding causal DNA variant could be observed in the RNAseq reads using the BAM-file (data not shown).
A minor discrepancy between the RT-PCR and RNAseq data was observed for replication sample 6. RT-PCR determined skipping of exon 18, while RNAseq in combination with QURNAs showed deletion of the first 5nt in combination with the intronic insertion of 24nt (Δ5nt 5'exon 18 + ins24nt intron 17; ERS 4.3 with 124 reads) and deletion of the first 5nt only (Δ5nt 5'exon 18; ERS 2.3 with 53 reads). Skipping of exon 18 is only observed with ERS 0.5 and 16 reads. After thoroughly reinvestigating the RT-PCR Sanger traces, the forementioned events could be observed in the background. An explanation for this difference may be preferential amplification of the smaller Δexon 18 fragment with RT-PCR. For sample 2, RNAseq identified besides the expected with RT-PCR detected deletion of the first 10nt of exon 18 (ERS 13.8, reads 160), additional exon 18 skipping (ERS 7.3, reads 75) and deletion of the first 5nt of exon 18 (ERS 3.0, reads 38, see Table 2 & Supplementary Data 3). After reevaluating the RT-PCR Sanger traces, these events could also be observed as minor events in the background.
Targeted enrichment platform comparison Part of the replication cohort samples was used to test the performance of QURNAs in combination with a different library preparation method to enrich for NF1 transcripts. PCR based targeting (NuGEN SPET) to enrich for NF1 was switched to a hybridization based enrichment (Agilent SureSelect). Comparable results were obtained with QURNAs using the different enrichment protocols. All splicing events were present using both methods with some variations in ERS and reads (Table 2 & Supplementary Fig. 5).
Testing molecular undiagnosed cases with a (suspected) clinical diagnosis of NF1 Eleven patients with a (suspected) clinical diagnosis of NF1, who lacked a definitive molecular diagnosis based on DNA diagnostics, were tested using targeted RNAseq (Sureselect, Agilent) and QURNAs data-analysis. All patients were checked for the highest ERS event(s) ( Table 3 and Supplementary Data 4 and 5) and a possible drop in the relative expression of the normal exon junctions ( Fig. 1b (patient Und-01), 2a (other patients)). When possible, heterozygous single nucleotide variants present in the coding DNA sequence were compared with the RNA sequence reads, to detect potential allelic imbalance (Table 3). For six patients no enriched pathogenic splicing event in NF1 nor SPRED1 (Supplementary Fig. 6 -only ERS for normal exon-junctions of SPRED1) could be detected based on high ERS, nor a prominent drop in ERS of the normal exonexon junction(s). Also, comparing the DNA and RNA sequencing reads, no evident allelic imbalance was observed for these patients. Interestingly, for two patients we were able to identify a deep intronic pathogenic splice variant or and for three patients resolve the exact effect on RNA splicing for variants (VUS, likely pathogenic) previously detected using DNA diagnostics (Table 3).
For undiagnosed patient 1 (in the same run as the replication samples), two major enriched splicing events between exons 38 and 39 were detected. These events introduce a cassette exon of 177 bp in intron 38 (acceptor site: ERS 51.3 and 1170 reads and donor site ERS 20.4 and 273 reads). In addition, a minor alternative acceptor site was also seen in the QURNAs output, resulting in a larger cassette exon of 208 bp at the 5'end but having the same 3'end (ERS 1.4 and 31 reads). The inclusion of these novel cassette exons, results in a drop in ERS of the normal 38-39 exon-exon junction (Figs. 1b, 2b). These aberrant splicing events are caused by a deep intronic c.5749 + 332 A > G change, that is visible in the RNAseq reads ( Supplementary Fig. 7). The presence of this variant was confirmed with Sanger DNA sequencing as a de novo pathogenic splice variant by testing the patient and parents. This variant introduces a strong splice donor site at this position and consequently the activation of the above mentioned cryptic splice acceptor sites (Fig. 2b). Both aberrant transcripts contain a predicted premature stop codon and are consequently expected to result in loss-of-function through either nonsense-mediated mRNA decay (NMD), or an abnormal protein product (Fig. 2b). Comparing various heterozygous single nucleotide variants shows modest allelic imbalance in the RNA sequence versus DNA, which indicates some NMD ( Supplementary Fig. 7). Of interest, the DNA variant is described multiple times and splicing was characterized by sequence analyses of cDNA clones with identical effect (both major and minor event) on RNA-splicing 33 , or similar effect on RNA-splicing (only major event) using direct cDNA sequencing 13,28 .  Patient 2 was included with a likely pathogenic (class 4) variant c.6580-2 A > G, which is predicted to show loss of the acceptor splice site that could lead to exon skipping. The splicing predictions however also indicate multiple cryptic acceptor sites upstream ( Supplementary Fig. 8) that could be alternatively used. Experimental analysis showed exon 43 skipping (ERS 8.5, 136 reads) and insertion of 17nt intron sequence through use of the first upstream cryptic acceptor (ERS 7.1, 97 reads) was observed ( Supplementary Fig. 8), leading to a drop in ERS of the normal exon-exon junction (Fig. 2a). The variant is visible in the RNAseq reads ( Supplementary Fig. 8). No allelic imbalance was observed (data not shown) despite absence of an NMD inhibitor in the culture, suggesting that the truncated protein is synthesized and not degraded at the mRNA level by NMD. Based on these findings the c.6580-2 A > G variant is now considered pathogenic (class 5).
For patient 5, a cassette exon including 42 bp of intron 11 was detected (acceptor site ERS 18.1 and 1497 reads, donor site ERS 28 and 2278 reads), which introduces a premature stop codon. (Supplementary Fig. 9). This event was found based on high ERS (Table 3, Supplementary Data 5) and a drop in ERS of the normal exon-11-12 junction (Fig. 2a). The aberrant transcript is caused by a c.1260 + 1604A > G change that is visible in the RNAseq reads ( Supplementary Fig. 9) and confirmed with Sanger DNA sequence analysis. Comparing various single nucleotide variants shows modest allelic imbalance in the RNA sequence versus DNA sequencing reads, which indicates some NMD ( Supplementary  Fig. 9). Interestingly, the pathogenic splicing variant c.1260 + 1604A > G is reported multiple times in the literature and found to have an identical effect on splicing 14,34,35 . Patient 6 has a clinical diagnosis of NF1 and a variant of unknown significance (VUS) c.2252 G > T (p.(Gly751Val)) in exon 19 of NF1. Her child also has the same variant and a clinical NF1 diagnosis. The glycine is conserved between species, the variant is not present in gnomAD and has, as annotated by Alamut strong disease causing in silico predictions. These predictions indicate partial loss of the (strong) acceptor site and introduction of a new donor site, potentially resulting in out-of-frame exon 19 skipping ( Supplementary Fig. 10). For this patient, RNA from short term cultured fibroblast was used, without NMD inhibition. QURNAs indicated exon 19 skipping, which is a naturally occurring event, hence the ERS is only 4.2 (covered by 4155 reads). The drop in ERS of the normal exon 18-19 junction is prominent though (Fig. 2a,  Supplementary Fig. 10). Of interest, comparing various single nucleotide variants, including c.2252 G > T, showed almost a complete loss of one allele in the RNA sequence versus DNA sequencing reads, suggesting NMD of the variant allele (Supplementary Fig. 10). This variant was previously also reported in another NF1-patient 36 . Together, these data support the classification of this variant as pathogenic.
For patient 8, RNA was directly isolated from peripheral blood lymphocytes, in contrast to the other samples in this run. This patient has a clinical diagnosis of NF1 and a VUS (class 3 variant) c.556 G > T (p.(Asp186Tyr)) in exon 5 of NF1. The aspartic acid is conserved between species and the variant is not present in gnomAD, as well as having strong disease causing in silico predictions (as annotated by Alamut). Splicing predictions show no effect on splicing through the introduction of a novel donor or acceptor site. Further in silico assessment showed effect on exonic splicing enhancer elements, and a slightly higher chance of exon skipping than the wild type (WT) allele were predicted ( Supplementary Fig. 11). Interestingly, literature suggests that the variant c.557 A > T immediately next to the one identified in this patient (c.556 G > T) is disrupting an exonic splicing enhancer element, thereby causing exon skipping 37 . Therefore this sample was included in this RNAseq study and found indeed to cause skipping of exon 5 only (ERS 3.6; 620 reads), though predominantly in combination with the inclusion of a cryptic intron between exons 4 and 5 (acceptor ERS 33.9 and 1153 reads and donor ERS 18.8 and 2739 reads; see Fig. 2b). Both events are naturally occurring albeit that the drop in ERS of exon boundaries 4-5 and 5-6 for this patient is more prominent compared to the others in the same run (Fig. 2a). The inclusion of the cryptic exon in combination with exon 5 skipping is most likely the consequence of this being an uncultured sample, since culturing of samples results in much lower levels of the cryptic exon (see Supplementary Data 2). Of interest, comparing the RNA versus the DNA sequencing reads showed an incomplete loss of the alternative T allele of the c.556 G > T variant (19% T allele present), caused by partial skipping of this exon. The other three informative single nucleotide variants do not indicate allelic imbalance (Supplementary Fig. 11), which suggests that the premature stop does not subject the transcript to strong NMD. In addition, residual presence of the missense variant p.(Asp186Tyr) could potentially have a functional effect. Together, these data support the classification of this variant as likely pathogenic. Further segregation and functional analysis should be performed to support pathogenicity.

DISCUSSION
Alternative splicing is an important mechanism wherein different mRNAs are generated from the same gene, thus increasing the coding capacity. It is often regulated in a tissue-or developmentalspecific manner. Deregulation of this mechanism can be caused by DNA variants in regulatory regions of RNA splicing, resulting in aberrant splicing. Skipping of exons and/or the introduction of a premature stop codon can lead to loss-of-function or even complete loss of the encoded protein from the variant allele. Several naturally occurring alternative transcripts have been described for the NF1 gene as well as aberrant splicing events due to DNA variants in this gene. In fact, a high percentage of up to~30% of pathogenic NF1 DNA variants have an effect on RNA splicing. Targeted RNA sequencing is currently more cost efficient compared to whole transcriptome for in depth analysis of the RNA splicing of a small set of genes. Previously, we used QURNAs to successfully detect BRCA1 and BRCA2 aberrant and normal splicing events in a small RNAseq data set 24 . Here we show that targeted RNAseq, in combination with QURNAs data-analysis (Supplementary Fig. 1), is a powerful approach to detect and distinguish normal and aberrant, pathogenic NF1 splicing events.
Not only could targeted RNAseq confirm previously reported and/or predicted changes in RNA splicing for all validation and replication samples, additional splicing events caused by the same DNA variant, but not (easily) resolved using RT-PCR, could be detected (Tables 1, 2). In these samples, with RNA extracted from lymphocytes, all normal exon-exon junctions and rare normal events previously reported in blood derived RNA could also be detected, except the neuronal tissue specific inclusion of a cassette exon of 30nt intron 11 ('ins 9a/9br'). Additionally, the use of two different enrichment methods gave similar results (Table 2, Supplementary Fig. 5).
After this thorough validation of our targeted RNAseq approach a total of 11 patients lacking a (clear) molecular diagnosis of NF1 was tested. Three patients had a clinical NF1 diagnosis and a NF1 DNA variant suspected to have an effect on NF1 splicing. The other eight patients showed some or clear NF1 symptoms, but lacked a molecular diagnosis through conventional DNA sequencing of NF1 (and SPRED1). (Likely) pathogenic NF1 splicing events could be detected in five out of these 11 patients using RNAseq (45%). For three of them this was a confirmation of a (likely) pathogenic splicing event caused by a DNA variant suspected to have an effect on splicing. For two patients (of eight without a prior molecular diagnosis) the pathogenic splicing event was caused by a deep-intronic DNA variant (25%), that was previously not detected using standard DNA diagnostics.
For six patients, no aberrant, pathogenic NF1, nor SPRED1 RNA splicing or expression could be detected. In most cases, bi-allelic gene expression could be confirmed, excluding the presence of a possible variant in the promoter region or transcription enhancer site that could result in strongly reduced or absent gene expression in the variant allele. To exclude other potential (rare) causes of NF1, the 5'UTR (up to c.−383) and the neuronal specific cassette exon of 30nt in intron 11 ('ins 9a/9br') of NF1 were also sequenced complementary, since these regions are not covered by RNAseq, nor by our standard DNA diagnostics procedure. Variants in the 5'UTR of NF1 creating or removing an upstream start codon (uAUG) and thus influencing translation can cause NF1 14,38 . However, no relevant variants were detected in these 2 regions.
In addition, retrotransposon insertions, that cause altered NF1 splicing, account for~0.4% of all pathogenic NF1 variants 17 . These are not targeted using conventional DNA diagnostics and not detected using RNAseq due to the exclusion of chimeric reads in in the RNAseq data analysis pipeline. However, these insertions can be picked up by Sanger based RT-PCR RNA sequencing, even though precise characterization of them requires additional customized efforts 17 . Another residual risk of missing variants in blood based DNA and RNA diagnostics is the presence of a mosaic postzygotic NF1 pathogenic variant, in specific affected tissue, but absent in blood. This might be the case for patient 4 ( Table 3) who presented with a segmental/mosaic NF1 phenotype restricted to the skin. Unfortunately, this affected tissue was not available for further DNA or RNA analysis. Finally, an NF1 (like) clinical phenotype, might also be explained by molecular defects in other genes, causing syndromes with overlapping phenotypes. These are, for example, Noonan syndrome (MIM PS163950; multiple genes/loci), Constitutional mismatch repair deficiency (MIM 276300; mismatch repair genes) and Proteus syndrome (MIM 176920; AKT1) 1,2,5 .
The QURNAs algorithm facilitated the identification of significant changes in RNA splicing. It uses an unbiased approach to identify and characterize aberrant and novel splice junctions without prior assumptions about transcript isoforms. The pipeline identifies all splice events present in the data and calculates their enrichment score and statistical significance compared to other samples within the same run. Together with a generic targeted RNAseq procedure, the QURNAs method can be used in routine diagnostics to discover, confirm or exclude effect of genetic variants on RNA splicing in NF1 or other genes of interest. RNAseq can potentially replace DNA based NF1 diagnostics as a primary test to detect clinically relevant variant. However, this requires the availability of reliable SNV calling on RNAseq data in the laboratory. Also, in patients for whom SNV calling may be compromised by alternative splicing events, e.g. exon skipping resulting in lower variant allele coverage, or degradation of the variant allele by NMD, DNA based diagnostics should still be run complementary to detect the causative variant.
In conclusion, the presented targeted RNAseq approach in combination with the computational QURNAs pipeline successfully detects and quantifies pathogenic NF1 splicing events driven by (deep) intronic or missense variants.

Ethics approval and consent to participate
All patients provided written or oral consent for either coded use of their material to study the clinical utility of RNAseq according to their diagnostic question or anonymous use of their material for improving diagnostic testing in general. The institutional Medical Ethics Committee (medischethische toetsingscommissie or METC) of the Maastricht University Medical Center confirmed that the Medical Research Involving Human Subjects Act (WMO) does not apply to our study (METC 2021-2666) and that an official approval of this study by the committee is not required.

Patient derived cells
White blood cells were isolated using ficoll (Sigma-Aldrich) from fresh whole blood (collected in EDTA tubes) and used either fresh or frozen in FCS with 10% of DMSO for subsequent short term-culture in a complete medium consisting of RPMI 1640 supplemented with 12.5% FCS, 1x L-glutamine (Gibco), 0.8 mM sodiumpyruvate (Gibco), 17 mM Hepes buffer (Gibco), 4.2 × 10-2 mM 2-mercaptoethanol (Gibco), 42 units/ml penicillin-streptomycin and 0.21 g/ml amphotericin B solution (Sigma-Aldrich). Lymphocyte growth was stimulated with 50 µl/ml PHA (Gibco) and 10 units/ml of IL-2 (Roche). At day 7, 4-6 h before harvesting the cells, each culture was split evenly and one part was treated with 200 µg/ml of puromycin (Sigma-Aldrich) as an inhibitor of nonsense mediated RNA decay (NMD) 39 .
For the initial validation of the targeted RNAseq procedure, RNA was isolated from uncultured peripheral blood lymphocytes of nine patients from the pre-implantation genetic testing (PGT) program at our institution, who had a pathogenic NF1 variant that leads to aberrant splicing (Table 1). Pheriphal blood lymphocytes were isolated and stored in DMSO in liquid nitrogen within 24 h after collecting the EDTA blood. Two additional samples from anonymous patients without a clinical diagnosis of NF1 were included as wild type controls. These validation samples were enriched for NF1 RNA using NuGEN SPET for RNA (see RNA isolation and Target enrichment/Library preparation).
The diagnostic performance of targeted RNAseq was further assessed using RNA samples shipped from the Centre for Medical Genetics, Ghent University Hospital, Ghent, Belgium (Table 2). RNA was isolated from blood lymphocytes (collected in EDTA tubes) that were short-term cultured. The molecular results of these samples were initially blinded, to test the ability to detect pathogenic NF1 splicing events without prior knowledge, comparable to the diagnostic setting. These replication cohort samples were enriched for NF1 using both Nugen SPET for RNA and for NF1/SPRED1 using Agilent SureSelect RNA Target Enrichment (see RNA isolation and Target enrichment/Library preparation). This allowed a direct performance comparison of these two different commercial platforms.
To further validate the targeted RNAseq approach and show its added value for the detection of pathogenic splicing events complementary to DNA diagnostics, an additional set of twelve patients was selected (Table  3). Eight patients were included because of having (or suspected of having) a clinical NF1 diagnosis without molecular confirmation. While two patients had a NF1 DNA variant with unclear effect on RNA splicing (variant of unknown significance -VUS [class 3]) and one patient had an NF1 DNA variant likely to affect splicing (likely pathogenic -LP [class 4]). One additional sample from an anonymous patient without a clinical diagnosis of NF1 was included. Available material from these patients, i.e. either lymphocytes directly isolated from blood or short-term cultured cells with or without NMD inhibition, was used to isolate RNA. Blood was collected in EDTA tubes and shipped and processed within a maximum of 48 h in our lab, except und8 that came into the lab > 48 h after blood draw. These samples were enriched for NF1/SPRED1 using SureSelect RNA Target Enrichment (see RNA isolation and Target enrichment/Library preparation).

RNA isolation and target enrichment/library preparation
Total RNA was purified and DNAseI treated using RNeasy Plus micro columns (Qiagen). Purified RNA was used as input for the SureSelect RNA Target Enrichment (below) and for the Ovation cDNA module followed by custom Ovation target enrichment, based on single primer enrichment technology (SPET) for RNA, according to the protocol of the manufacturer (NuGEN). Primer design for NF1 (NM_000267.3) was performed by the manufacturer for the coding regions and 5'/3'UTR with a probe spacing of 150 bp. Library preparation was according to the manufacturer instructions.
Additionally, the SureSelect RNA Target Enrichment for Illumina Paired-End Multiplexed Sequencing kit (Agilent; protocol version 2.2.1) was used to target the NF1 (NM_000267.3) and SPRED1 (NM_152592.2) coding regions and 5'/3'UTR. Probe design was performed using the Agilent eArray online tool, setting the region of interest to minus 100nt of the normal 5'exon boundary and plus 150nt of the normal 3'exon boundary and 6x tiling of probes. Library preparation was according to the manufacturer's instructions. Primers (SPET) and/or probes (Sureselect) design is available upon request.

Read alignment
RNAseq reads were mapped with the STAR mapper (Version 2.4.1d) and an index created from HS.GRCh37 with 92 bp overhang 40 . We performed unique mapping onto the genome without trimming and with max. 10 mismatches per read while no chimeric reads were allowed (default STAR settings). To elucidate the PCR amplification bias, de-duplication was performed either using a build in Unique Molecular Identifier (UMI) in the probe design (NuGen SPET) or by removing identical reads (Agilent Sureselect). Duplicate read pair detection was carried out with Picard tools (https://github.com/broadinstitute/picard). Start and end positions from STAR output refer to the first nucleotide in the intron (AG | gu) and last nucleotide of the intron (ag | G), respectively.

Quantitative enrichment of aberrant splicing events in targeted RNAseq-QURNAs
The QURNAs computational pipeline was developed to identify samplespecific splicing events in targeted RNAseq data, including novel and possibly aberrant transcripts or increased abundance of normal isoforms (see Supplementary fig. 1). To minimize the number of possible artefacts or non-relevant events in the QURNAs data analysis, it is important to perform QURNAs analysis on a group of at least 4 samples that are equally processed. This way, events that are either induced by sample specific treatment, and/or are more or less equally present in all samples from the group, get a neutral enrichment score of~1 in the QURNAs output after normalization and intersample comparison (see below for methodology). QURNAs is available at https://dataverse.nl/dataset.xhtml?persistentId=hdl:10411/LY8ZQ4. The mapping of short sequences to the reference genome without prior assumptions about transcript architecture, enables statistical evaluation of significance of observed RNA isoforms.

De novo prediction of splice junctions from RNAseq data
QURNAs first creates a collection of splice sites based on the reads mapped onto the genome (BAM file). For these splice sites, reads in BAM files that overlap with the splice site are counted. While iterating over all reads in the region, two additional features are counted to estimate the fraction of the reads in the region congruent with the splice junction. First, the number of reads is calculated that overlap with and contain the last nucleotide of the exon at the donor site. This includes, next to reads with a different acceptor site, "unbroken" reads corresponding to intron retention events. The second count consists of reads that "span" the splice site, thus start upstream of the nucleotide and end downstream, but do not necessarily contain the nucleotide next to the donor site in their sequence (e.g., exon skipping events or splice events that use an alternative upstream donor site).

Enrichment score (ERS) calculation
Enrichment for a splicing event in a sample is calculated as follows. For each splice site, a relative read count l spl is calculated that measures the fraction of aligned reads that are congruent with the splice site among all the reads that span the splice site. Formally, we count reads (r spl ) that are congruent with the splice spl=(start,end) i.e., the genomic alignment of the read contains a gap between exact (start, end) positions, allowing for other insertions/deletions in the read. Reads for each of the splice r spl are considered relatively to reads that span the splice site S spl (spanning reads). These start their genomic alignment upstream of the donor site and end downstream of this site but do not necessarily share the same donor or acceptor sites. To minimize the impact of spurious reads and sequencing noise, pseudocounts of at least 10 reads and up to 1% of spanning reads (max(10,S spl 1%)) are added to both r spl and s spl . The relative read count l spl for each sample is then calculated using Eq. (1). In sample i, after calculating the relative splice read count, enrichment of reads m spl i is calculated by comparing the relative read count to the average read counts from other samples using Eq. (2). To prevent samples with low sequencing coverage having inflated relative read counts due to pseudocounts (e.g. l spl i ¼ 1 for the extreme case of absent coverage, i.e. r spl i ¼ s spl i ¼ 0), the enrichment m spl i is calculated using l0 spl i without pseudocounts.

Statistical significance of enriched isoforms
To test whether splice isoforms are expressed at significantly higher levels in a sample, a two-tailed proportion test that tests the null hypothesis that proportions in several groups are the same is used. For each splice junction a group of trials as number of spanning reads (s spl ) that includes reads that support the junction r spl (positive outcome of the trial) is defined. If significant (p(spl) < 0.01), the proportion test rejects the null hypothesis that the expression of the splice junction is the same in all samples. The test assumes independence between trials in the group and this is not the case for RNAseq due to read redundancy and overlap between 150-bp reads. The marginal information contributed by additional reads, drops with the number of reads mapped to a splice junction. To reflect lower marginal information, we use square root of RNAseq read numbers using Eq. (3) where i corresponds to subsequent samples. To determine the significance of a splice event in a specific sample i, the reads for the sample are compared to the background distribution of reads in remaining samples, testing the enrichment hypothesis (one-tailed test expecting enrichment) using Eq. (4). Note that in this statistical framework it is possible that p(spl) < 0.01 but p(spl i ) > 0.01 for all samples i (the distribution of reads is different between samples, but none are significant). This is, for example the case for WT splice isoforms in mutated samples where the WT isoform becomes depleted (p(spl i ) tests only for enrichment).

Interpretation QURNAs output
The QURNAs output for all splicing events was reduced with an additional GAP filter of >3 ((pos end_STARpos start_STAR) + 1), to exclude artefacts that do not represent splicing since the intron size would be smaller than 4 nucleotides. The absolute read counts of all normal exon splice-junctions in the reference transcript were inspected to confirm good sequence coverage throughout the gene for all samples. The ERS values of these reference exon splice junctions were subsequently checked for a possible drop (normal reference exon splice junctions are expected to have a median ERS~1). Highly expressed alternative splicing events usually result in lower expression of the normal exon splice junctions in the same region, reflected by such a drop (ERS <0.8 or the lowest value of the total group of samples). Subsequently, all splicing events for a specific sample were sorted high to low based on the ERS. All events with ERS > 5 were further evaluated to discriminate possible artefacts from true alternative splicing events. Importantly, relevant ERS values can be lower than 5 when; 1) the splice variant gives rise to multiple different splicing events, instead of one or two major events; 2) the splicing events concerns an upregulated normal splicing event which is also present in other samples in the run at lower expression levels or 3) when 2 samples are present in the same run with the same effect on splicing. Alamut v2.15 (Interactive Biosoftware) was used to identify true splice-donor and acceptor sites based on in silico prediction. Regions with alternative splicing were further inspected in the BAM file in a genome browser (IGV and Alamut) to check for possible causal sequence variants in close proximity. Allelic imbalance was assessed by comparing heterozygous variants in the DNA vs the RNA sequences in IGV. A minimum of at least 2 heterozygous variants in the coding regions and/or UTRs was considered informative to assess allelic imbalance. Consistent variant allele ratio's in the sequence from 60/40% to 70/30% were considered as modest imbalance, 70/30% to 80/20% as intermediate and >80/20% as strong. But these boundaries are set arbitrarily for now and need further validation with more samples, and are also under the assumption that the imbalance is for the same allele. Since short-read sequencing was performed, the SNP phase is not known.

Nomenclature
The description of genetic variants follows the Human Genetic Variation Society (HGVS) approved guidelines 41 , where c.1 (and r.1) is the A of the ATG translation initiation codon. Alternative splicing events are those incorporating splice junctions not present in the reference transcript NM_000267. Exon numbering used by QURNAs and in the tables is systematically 1-57 using NM_000267, for conversion according to the NF1 best practice meeting 42 this numbering is +1 for all exons after exon 30 (LRG_214).

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.