Benchmarking RNA-seq Tools for Real-World Diagnostic Applications

Benchmarking RNA-seq Tools for Real-World Diagnostic Applications | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Benchmarking RNA-seq Tools for Real-World Diagnostic Applications Sarah Silverstein, Kaushik Ganapathy, Sandra Donkervoort, Veronique Bolduc, and 11 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-8622051/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Background Pediatric neuromuscular diseases are genetically and clinically heterogeneous. A substantial proportion remain without a definitive genetic diagnosis despite available clinical molecular testing. RNA-sequencing (RNA-seq) can be used to complement genome or exome sequencing to elucidate or to identify the functional impact of variants of uncertain significance, but when manually analyzed is limited to candidate DNA variants or phenotype-driven gene lists. Open-source computational tools have been developed to systematically and unbiasedly analyze RNA-seq data for aberrant splicing, expression, or allelic imbalance. However, best use practice of these tools is yet to be established. Methods To assess the performance of selected tools, we collected RNA-seq from 97 previously diagnosed samples to establish a truth set for benchmarking. Pathogenic variants were categorized as: true positives with confirmed aberrant RNA events and true negatives with no transcriptomic effect. We assessed performance of eight commonly used tools for splicing, expression and allelic imbalance analysis. We then applied the optimal strategy to 74 undiagnosed RNA-seq samples to identify new candidate diagnoses. Results Across 68 diagnosed probands with aberrant RNA events, tools correctly identified 28 diagnoses. Splicing analysis tools provided most of the findings, but allelic imbalance tools uniquely identified 4, underscoring their value. Conversely, the false positive rate was highest for the splice tools and lowest for expression analysis. Application of tools led to identification of candidate variants for only 9 out of 74 undiagnosed patients. Conclusions Inclusion of RNA-seq tools can expedite variant prioritization, characterization and interpretation in the diagnostic pipeline but remain complementary to manual analysis of loci where candidate variants were identified by DNA sequencing. RNA-seq transcriptomics rare disease diagnostics pediatric neuromuscular disease computational tools variant interpretation Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Background Pediatric neuromuscular disorders are characterized by broad genetic and clinical heterogeneity, with approximately 30–60% of cases remaining molecularly undiagnosed following exome or genome sequencing[ 1 – 4 ]. Recent reports suggest that incorporating RNA-seq data derived from clinically accessible tissues such as whole blood, peripheral blood mononuclear cells, fibroblasts, urine, or muscle can increase diagnostic yield by an additional 2.7–32%[ 5 – 17 ]. In a typical clinical-diagnostic bioinformatics workflow of DNA sequencing, variants are prioritized based on phenotypically-driven analysis, possibly identifying candidate variants with effect on splicing or gene/protein expression. If RNA-seq data is obtained, splicing and expression of the genes with these variants are then manually inspected in the RNA-seq BAM file data. This process is time-consuming and, despite the availability of published guidelines, remains heavily reliant on the expertise of individual analysts and the availability of appropriate control and tissue specific datasets[ 18 ]. Recently, open-source computational tools have been developed to automate the analysis of transcriptomic data in the context of rare disease[ 19 – 24 ]. These tools are designed to identify aberrant RNA events—such as abnormal splicing, aberrant gene expression, or allelic imbalance—in a systematic and agnostic manner. Unlike manual methods, these tools employ an RNA-first strategy, enabling transcriptome-wide analysis independent of prior DNA variant prioritization. This approach has the potential to identify pathogenic events that may be overlooked during manual review. Applying open-source strategies to RNA-seq datasets derived from clinically accessible tissues has resulted in modest diagnostic yield, with the highest yield in cases where initial DNA sequencing identified a variant of uncertain significance (VUS) predicted to alter splicing and prioritized [ 5 – 7 , 15 , 17 , 25 – 28 ]. In cases without any diagnostic leads, RNA-first strategies have resulted in limited diagnostic yield [ 5 – 7 , 15 , 17 , 25 – 27 ]. To date, comprehensive benchmarking of these tools using previously diagnosed samples harboring pathogenic variants leading to known aberrant RNA events has yet to be undertaken. Furthermore, prior studies have typically sequenced large batches of samples and analyzed them collectively to mitigate batch effects. This strategy does not reflect real-world clinical constraints, where disease-relevant tissue is often limited, requiring evaluation of open-source tools using diverse approaches. We assessed a range of RNA-seq analysis tools to identify those best suited for our diagnostic pipeline and to assess their limitations. The most effective strategy was then applied to a set of clinically well-characterized, undiagnosed cases to enable candidate variant discovery. We also evaluated the capacity of RNA-seq to function as a standalone diagnostic tool, versus being utilized as an adjunct to DNA-based approaches. Methods Patient Recruitment Patients were referred by their neurologist or geneticist. Written informed consent and age- appropriate assent for study procedures were obtained by a qualified investigator [protocol 12-N-0095 approved by the National Institute of Neurological Disorders and Stroke, National Institutes of Health IRB]. Medical history was obtained, including clinical evaluations and muscle biopsy, were performed as part of the standard diagnostic examination. Samples for research- based testing, including skin, blood and muscle biopsy samples, were obtained using standard procedures. Patients with confirmed genetic diagnosis were included as true positive or true negative benchmarking samples, and patients without candidates after evaluation were included for our undiagnosed samples analysis (Supplemental Methods). RNA-seq Alignment Quality control of FASTQ files were assessed with FastQC (v0.12.1). BAM files were generated using the GTEXv10 pipeline ( https://github.com/broadinstitute/gtex-pipeline ) using reference genome GRCh38 (Gencode, version 39). Briefly, FASTQ files were aligned using STAR (v2.7.10a) and optical duplicates marked by Picard Mark Duplicates (v3.3.0)[ 29 ]. Alignment metrics were assessed with RNA-SeQC (v2.4.2)[ 30 ]. Principal component analysis was performed using Transcript per million (TPM) values generated by RNA-SeQC and the prcomp function in R package stats (v4.4.0). Samples not clustering with other samples annotated to be derived from the same tissue source (either muscle, fibroblast or cycloheximide treated fibroblast) were removed (Fig S1 ). True Positive Definition Diagnosed samples that passed quality control and contained diagnostic aberrant RNA events were manually reviewed in the Integrative Genomics Viewer (IGV)[ 31 , 32 ]. Each sample contributed one true positive gene for evaluation. Each true positive gene could be included as a true positive across more than one analysis type ( e.g . splicing and expression) based on underlying mechanism. Splicing true positives were defined as containing ≥ 3 reads supporting a novel splice junction created by the pathogenic variant. Expression true positives were defined by a confirmed loss-of-function mechanism with normalized transcript abundance ≤ 50% or ≥ 100% of the corresponding normalized GTEx TPM. Genes were normalized to the geometric mean of HPRT1 , PPIA , YWHAZ and TBP . Allelic imbalance true positives were defined by monoallelic loss-of-function events supported by single nucleotide polymorphism (SNP) analysis demonstrating a marked allelic ratio ≥ 70:30 at the site of the rare variant. True Negative Definition Patients harboring pathogenic variants with no direct effect on splicing or expression (missense, in frame indels, terminal gene frameshifts) were considered as “true negatives” for RNA anomaly. RNA-seq data for these patients were manually reviewed in IGV to confirm that splicing was canonical and that no expression changes or allelic imbalance were present. Computational Tools All tools were run on Biowulf ( https://hpc.nih.gov ). Current RNA-seq tools for rare disease use two main strategies: outlier-based (one-vs-all) and event-based (case-vs-control) approaches [ 33 ]. Outlier-based methods assume that each rare disease is unique within a heterogeneous cohort, an assumption that holds for some muscle disorder cohorts. However, it was not a valid assumption for our dataset, which includes recurrent variants and genes. Additionally, small-batch sequencing within the cohort likely introduced variability across samples. Event-based tools that compare cases to matched controls were considered as a strategy to mitigate these issues. We attempted to evaluate both case-versus-control and one-versus-all strategies when possible; however, this was only possible for splicing analysis in practice. One vs All Splice Tools (Outlier based approaches) Three popular tools—FRASER, FRASER2.0, and LeafcutterMD [ 20 , 23 , 34 ] —were selected for evaluation. Since these tools require ≥ 30 samples, we created sample groups by tissue of origin, mixing diagnosed and undiagnosed RNA-seq samples from the same tissue: muscle (n = 71), fibroblasts (n = 37), and cycloheximide-treated fibroblasts (CHX, n = 79). Muscle samples were sequenced in multiple small batches over 14 years using different library preparations (Table S1 ), primarily from the quadriceps (Fig S2 ). Although FRASER recommends analyzing stranded and unstranded libraries separately, preliminary analysis showed higher true positive recall when both library preparation methods were combined (Fig S3 E-F). Therefore, all three tools were run with mixed library preparation. FRASER was executed via the DROP pipeline (v1.2.2), using default parameters and a false discovery rate (FDR) of 5% (Fig S3 A-C) [ 24 ]. FRASER2.0 was similarly run via the DROP pipeline (v1.4.0) with standard settings and an FDR of 5% (Fig S3 D) [ 24 ]. LeafcutterMD (v0.2.9) was applied with preprocessing of BAM files to generate intron clusters using Leafcutter, with recommended parameters (minimum 50 splice junctions per cluster and a maximum intron distance of 500 kb)[ 34 , 35 ]. Outlier detection was restricted to groups with at least 20 reads and a minimum intron cluster size of 50. Resulting junctions were annotated with Ensembl gene IDs using Bedops (v2.4.41), followed by conversion to gene symbols utilizing the annotationDBI package in R. Correction for multiple hypothesis testing was performed using a FDR of 5%. Case vs Control Splice Tools (Events-based approaches) Two commonly used tools, rMATS-Turbo and MINTIE, were selected[ 19 , 22 ]. Each patient sample was run against a set of ten healthy control samples. We used ten samples for the control cohort based on MINTIE recommendations. Control samples were selected from the GTEx database to match the tissue source when available (Table S2 ). Muscle control samples were filtered for RNA Integrity Number (RIN) > 8, ischemic time < 1000 minutes, and absence of autolysis. Fibroblast controls were selected based on tissue type "skin" and tissue type detail "cells – cultured fibroblasts," further filtered for RIN = 10 and ischemic time < 400 minutes. Due to the lack of available disease-free controls for CHX samples, case vs controls approaches were not tested. MINTIE (v0.4.1) was run using FASTQ files from each patient sample against the ten selected GTEx controls, with recommended parameters and additional settings (splice_motif_mismatch = 3, FDR = 5%)[ 22 ]. Output was filtered to include only splice-related outliers, excluding large deletions, insertions, and other structural variants. rMATS-Turbo (v4.1.2) was executed on BAM files, with each patient sample compared to ten controls[ 19 ]. Standard parameters were used, with additional flags --novelSS and --allow-clipping. The maximum exon length for novel splice isoforms was set to 700 bp (-mel 700) to accommodate detection of large pseudoexons based on the largest in our true positive set. Correction for multiple hypothesis testing was performed using an FDR of 5%. All five MATS output files were included for assessment of true positive events, regardless of whether the detected aberrant splice category matched the expected event (e.g., intronic cryptic donor activation detected as intron retention). Expression Tools OUTRIDER was run on all three cohorts as previously described in a one vs all manner via the DROP pipeline (v1.2.2) with recommended parameters which included an FDR of 1%[ 21 , 23 ]. Since OUTRIDER is designed as a one vs all approach, we performed a second run for the muscle patient samples with an addition of 100 control muscle samples downloaded from GTEx to test if the influx of healthy samples improves detection (Table S2 ). Allelic Imbalance Tools The approach to detecting allelic imbalance in transcriptome data relies on heterozygous presence of variants in DNA. Thus, we restricted our test set to muscle samples where exome or genome variant call files (VCFs) were readily available, for a subset cohort of 59 samples. Two tools were selected, MAE module in the DROP pipeline [ 24 ] which uses variant level allele specific expression (ASE) and ANEVA-h, which uses gene level ASE[ 36 ]. The MAE module in DROP (v1.4.0) was run using standard parameters except modification of the allelic balance threshold to 70 instead of 80, reflecting our true positive dataset[ 9 , 24 ]. ANEVA-h required preprocessing of VCF and BAM files as follows. Haplotype-resolved, gene-level ASE data were generated for each case to enable regulatory outlier detection using ANEVA-h. RNA-seq libraries were processed and aligned to the GRCh38 reference genome as previously described. Corresponding genotype data were obtained from matched genome sequencing (GS) and phased using EAGLE (v2.4.1), leveraging the 1000 Genomes Project Phase 3 (30x) high-coverage reference panel for population-based phasing[ 37 , 38 ]. Phased genotypes and aligned RNA-seq BAM files were used as input to phASER (v1.1.1), which was run following best practices for haplotype-resolved ASE quantification [ 39 ]. Gene-level ASE values were computed using the phaser_gene_ae.py module, with heterozygous exonic SNPs aggregated by haplotype using GENCODE v39 transcript annotations [ 39 ]. ANEVA-h was then applied using the ANEVADOT_test function from the ANEVA-DOT R package, with default parameters [ 40 ]. Analyses were performed using reference VG estimates derived from 706 skeletal muscle RNA-seq samples from GTEx v8, ensuring tissue-specific calibration ( https://zenodo.org/records/15226130 , https://zenodo.org/records/15226575 ). A minimum coverage threshold of 10 reads per gene was applied for inclusion. Multiple testing correction was performed using the Benjamini–Hochberg FDR procedure, and genes with FDR-adjusted p-values < 0.05 were designated as significant ASE outliers. True Positive and Negative Detection Tool performance was evaluated based on detection of known aberrant RNA events. For splicing, a correctly detected true positive was defined by an exact coordinate match to the validated splice event. For expression, the true positive gene had to appear among reported outliers. For allelic imbalance, MAE required inclusion of the correct rare variant, while ANEVA-H required inclusion of the corresponding gene. True negatives were assessed by the absence of reported outliers, with correct detection defined as no outliers identified in the sample. Performance Metrics Confusion matrices for True Positive (TP), True Negatives (TN), False Positives (FP) and False Negatives (FN) were established for each tool, allowing calculation of Precision, Accuracy, Recall, False Positive Rate (FPR), F1 score, and Mathews Correlation Coefficient (MCC). $$\:\text{Recall}=\frac{TP}{TP+FN}$$ $$\:\text{Precision}=\frac{TP}{TP+FP}$$ $$\:\text{Accuracy}=\frac{TP+TN}{TP+TN+FP+FN}$$ $$\:\text{FPR}=\frac{FP}{FP+TN}$$ $$\:\text{MCC}=\frac{TP\times\:TN-FP\times\:FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$$ $$\:F1=2\times\:\frac{\text{Precision}\times\:\text{Recall}}{\text{Precision}+\text{Recall}}$$ Variant Annotation Individual VCFs generated from exome sequencing and genome sequencing were annotated for population allele frequencies and splice prediction scores using the VEP_annotate tool (v113)[ 41 ]. SpliceAI scores were obtained as precomputed datasets from Biowulf[ 42 ]. Allele frequencies for each variant were sourced from gnomAD genomes v4.0 ( https://gnomad.broadinstitute.org/ ). Subsequent filtering of annotated VCFs was performed using VEP_filter (v113) to retain only variants with an allele frequency less than 0.01 or 0.0005 (a high and low threshold relevant to rare disease to assess how different cutoffs affect final outlier number) and with no homozygotes reported in gnomAD. Filtered VCFs were then converted to BED format using Bedops (version 2.4.41 megarow) with the convert2bed function. Outlier Evaluation For splicing, an ensemble approach with variant filtration was used. Chromosomal coordinates for each outlier splice event were extracted into a master bed file (see Supplemental Methods). The master bed file was sorted using bedops (v2.4.41) and intersected with the filtered annotated variant bed file using the --range − 50:50 -e set operation to identify variants overlapping splice outliers. Variants that intersected the outlier bed file were then filtered for spliceAI score ≥ 0.1, segregated in seqr and manually evaluated with the integrative genomics viewer (IGV)[ 32 , 43 ] (see Supplemental Methods). For OUTRIDER and ANEVA-h, the number of variants identified per individual was sufficiently low to allow for manual inspection of all events without additional filtration. Candidate gene variants identified for undiagnosed patients were submitted to GeneMatcher[ 44 ]. Results Truth set selection We used RNA-seq data from 97 previously diagnosed patient samples to benchmark RNA-seq analysis tools. To ensure RNA-seq samples were correct, each BAM file was confirmed to harbor the known pathogenic variant in IGV. RNA-seq samples with independently verified variants affecting RNA were manually validated to contain evidence of aberrant splicing, expression or allelic imbalance in IGV. This set was defined as our true positive set (76 samples from 68 probands, Fig 1A,C, see methods for detailed definition). RNA-seq samples with confirmed pathogenic missense variants, inframe indels and premature truncating variants without nonsense mediated decay (NMD) and not expected to cause RNA changes were manually confirmed to contain no visible transcriptome changes in the pathogenic gene. This set was defined as our true negative set (21 samples, Fig 1A, see methods for detailed definition). Recall rate of splice tools identifies ensemble strategy as most sensitive To determine the performance of splice tools, we compared five tools using both case–control and one-vs-all designs to recall true positive splice events by tissue source (68 true positive variants from muscle, fibroblast and cycloheximide treated fibroblasts, Fig. 1A-B). In muscle derived true positive samples, case–control tools were most sensitive: rMATS (56%) > MINTIE (37.5%) > FRASER (33%) > FRASER2.0 (12.5%) > LeafcutterMD (0%) (16 true positive variants, Table S3, Fig. 2A). FRASER2.0 was excluded from further analysis as it underperformed FRASER. In fibroblast derived true positives, rMATS again performed best (86%) > MINTIE (28%) > FRASER (14%) > LeafcutterMD (0%) (7 true positive variants, Table S4, Fig 2A). Cycloheximide (CHX) treatment of fibroblasts can halt NMD by blocking protein synthesis, increasing detection of degraded aberrant events. In CHX-treated fibroblast true positives, only outlier-based methods were tested due to lack of control data, FRASER detected 33% and LeafcutterMD none (45 true positive variants, Table S5, Fig 2A). To evaluate if underperforming tools contributed unique detections, we examined recall by splice event category. FRASER detected two intron retention events, contributing one variant missed by other tools, supporting its inclusion (Fig 2B). MINTIE uniquely captured multiple splice events of a BET1 variant (NM_005868.6: c.202G>C; p.Asp68His) in fibroblast and muscle tissue (Fig 2B, Tables S3-4). Although case-control methods have the highest recall rates, an ensemble approach aggregating outliers across multiple tools is most sensitive. In muscle, when tool output is combined, 68% (11/16) true positives are recalled, indicating greater overall yield with ensemble strategies (Fig 2B). To understand the limitations of tool use, we examined the true positives undetected by any tool (Tables S3-5). These reflected common limitations: low mapping quality ( POMGNT1 NM_001290129.2 c.636C>T; p.205 =), few spanning splice junction reads due to degradation ( RYR1 NM_000540.3: c.3178+587A>G; NEB NM_001164508.2: c.21417+3A>G), complex loci with age specific alternative exon usage or highly repetitive regions ( NEB NM_001164508.2: c.21417+3A>G, IGHMBP2 NM_002180.3: c.1235+450 G>A), or inaccurate splice junction splitting ( RYR1 NM_000540.3: c.3178+587A>G , COL6A1 NM_001848.3: c.190+176C>T, IGHMBP2 NM_002180.3: c.1235+450 G>A). These cases demonstrate the continuing importance of manual investigation guided by clinical and genetic hypothesis. Clinical application of splice detection tools requires a manageable number of outliers for review. We assessed output interpretability by comparing outlier counts reported on average per sample. In muscle, MINTIE and rMATS reported thousands of outliers, particularly when strand-specific data were compared to unstranded controls (Fig. 2C). In fibroblast and cycloheximide-treated fibroblast samples, outlier counts clustered by sequencing center, reflecting batch and library preparation effects (Figs. S4A–D). To reduce outlier burden, we filtered for events overlapping rare DNA variants (Fig. 2D). Among ten true positive samples with genomic data, filtering for AF ≤ 0.01, a lenient cutoff for rare disease, yielded a mean of 1,221 variants per sample; applying a stricter threshold (AF ≤ 0.0005, segregation when available, and SpliceAI ≥ 0.1) reduced this to 1–3 variants (Figs. 2D–E, Table S6, Supplemental Methods). In nine of ten cases, the causal variant remained post-filtering. The missed case ( POMT2 NM_013382.7) the first variant (c.333+1G>T) associated splice event was not included in the master bed file (see Supplemental Methods), while a second intronic variant (c.1725+78G>A) was not captured by exome sequencing, highlighting the added value of including RNA-based variant calling (Table S6). Overall, rare variant filtration is successful at overcoming the challenge presented by the high number of outlier calls. Recall rate of expression tools is low To determine the utility of expression analysis, we benchmarked OUTRIDER (15 true positive genes, Fig 1A-B). Across muscle true positive cases the recall rate was 7%, and the addition of 100 healthy control samples did not improve the recall (14 true positive genes,Table 1, Table S7). Fibroblast true positives were not detected at all, and no true positives were assessed from the CHX samples (1 true positive gene, Table 1, Table S7). Overall, the combined recall rate was 6.6% (1/15), and inclusion of expression analysis did not identify diagnoses missed by other analysis types (Fig 1C). The average number of outliers reported per sample were low (Table 1). Notably, the inclusion of control samples increased the average number of outliers while not providing additional diagnostic benefit (Table 1). To understand the limitations of OUTRIDER, we evaluated the undetected true positives. TTN underexpression recurs four times in our group of true positives and may not be considered statistically significant when using the one vs all approach. No patterns emerged for why any of the other true positives were missed (Table S7). Recall rate of allelic imbalance highlights ANEVA-H To benchmark allelic imbalance detection, we tested the MAE module from DROP and ANEVA-H (8 true positive genes, Fig. 1A-B). Both detected 5 of 8 true positives (62.5% recall; Fig. 3A, Table S8). The tools overlapped on four variants, each identifying one unique event, for a combined recall of 75% (6/8). Allelic imbalance analysis yielded 4 diagnoses missed by splice or expression outlier methods (68 total diagnoses, Fig. 1C), highlighting its complementary diagnostic value. ANEVA-H reported markedly fewer outliers per sample than MAE (mean 44 vs 167; Fig. 3B), supporting its clinical adaptability. Undetected cases involved MSTO1 and JAG2 . MSTO1 ( NM_018116.4:c.1433A>G; p.Tyr478Cys & deletion) mapping is confounded by a pseudogene, and JAG2 ( NM_002226.5: c.2515G>A; p.Gly839Arg & large deletion) was listed as homozygous in the VCF, masking monoallelic expression (Table S8). Two true positives were only detected by one tool each. ANEVA-H specifically missed RYR1 (NM_000540.3: c.14126C>T; p.Thr4709Met & c.3178+587A>G) due to absent haplotype-specific SNPs in RYR1 (required for gene level ASE generation), and MAE specifically missed TTN (NM_001267550.2: c.35794G>T; p.Glu11932Ter & c.22973C>A; p.Ser7658Ter) because of its size and repetitive regions (Table S8). Gene-level estimates from ANEVA-H may be more sensitive to subtle effects than variant-level MAE output. True negative analysis reveals high false positive rates Unnecessary evaluation of RNA outliers in patients whose final diagnoses do not involve RNA-mediated mechanisms can divert limited clinical resources. To estimate how often such outliers occur, twenty-one true negative (defined in methods) RNA-seq samples were analyzed across all computational tools for presence or absence of outlier calls, as presence of any number of outliers will require clinical evaluation without benefit (Fig 1A). Among splicing tools, rMATS, MINTIE, LeafCutterMD, and FRASER reported outliers for all true negative samples (100%), while FRASER2.0 had a lower rate (33%) (Fig. 4A, Tables S3-5, Fig. S5). For expression tools, OUTRIDER with GTEx controls showed higher false positive rates than OUTRIDER alone (83% vs. 56%), confirming that OUTRIDER alone is both more sensitive and specific (Table S7, Fig 4A). For allelic imbalance tools, MAE produced more false positives (100%) than ANEVA-H (80%), further supporting ANEVA-H as the better-performing method given its lower outlier burden (Table S8, Fig. 4A). Poor predictive performance across tools To evaluate tool accuracy and precision, true positive and true negative detections were aggregated across tissue types for each tool (Table S9). Among splicing tools, rMATS achieved the highest precision, accuracy, recall, and F1 score (Fig. 4B, Table S9). For expression analysis, OUTRIDER alone outperformed the OUTRIDER + GTEx combination across all metrics (Fig. 4C, Table S9). In allelic imbalance analysis, ANEVA-H exceeded MAE on every metric (Fig. 4D, Table S9). Notably, Mathews Correlation Coefficient (MCC) values were negative across all tools, indicating limited overall predictive power (Fig 4B-D, Table S9). Across all analyses, only 28 of 68 (41%) known diagnoses were correctly detected (Fig. 1C), consistent with the MCC scores. Application to undiagnosed samples We applied the top performing approaches integrating ensemble splice analysis with variant filtration, OUTRIDER, and ANEVA-H to 74 undiagnosed samples (Fig. 5A). In 28 muscle-derived samples, splice analysis identified 129 candidate variants (mean 4.6 per sample), of which 39 showed clear aberrant splicing on manual review. Nine variants were prioritized for further evaluation after segregation and analysis for a possible second hit if not de novo (Figs. 5B-C, Table S10). None were located in OMIM-listed disease genes. OUTRIDER analysis identified one candidate in a CHX treated sample, TRIP4 , with reduced expression, in a patient with phenotypically compatible disease (Fig. S6A–B); however, exome, genome, and optical mapping revealed no causal variants, and follow-up RNA-seq confirmed normal expression and splicing in untreated fibroblasts (Fig. S6C–D). ANEVA-H did not identify additional candidates. Together, these analyses demonstrate that splicing-based approaches contributed novel candidate findings, whereas expression and allelic imbalance analyses did not improve diagnostic yield. Discussion Over a 14-year period, RNA-seq was integrated into our research diagnostic pipeline for patients with suspected Mendelian neuromuscular disorders, improving diagnostic yield but leaving a substantial fraction unresolved. To assess how transcriptome-centric methods could further enhance yield, we benchmarked commonly used RNA analysis tools. Splicing tools consistently demonstrated the highest diagnostic sensitivity, detecting 24/68 diagnoses, in line with prior studies [ 5 – 7 , 14 , 15 ]. While outlier-based methods FRASER and FRASER2 are the most widely implemented, case–control approaches such as rMATS-turbo achieved greater sensitivity in our dataset. A prior report analyzing whole blood RNA-seq reported similar findings, while another report using muscle-derived RNA-seq found greater sensitivity with FRASER [ 5 , 7 ]. Differences in FDR thresholds applied to rMATS and FRASER may explain the discrepancy [ 5 ]. Although combining multiple tools improved recall in our dataset, it also increased outlier burden. Reported strategies to reduce outliers include restricting to clinically relevant genes or panels [ 5 , 6 , 15 , 16 ]. While this enhances interpretability, it does not consistently result in new diagnoses [ 15 – 17 ]. We chose not to restrict results this way given that our patients had undergone prior extensive phenotype based genetic analyses without success. Treatment of cells with cycloheximide has been proposed to enhance detection of splice events subject to NMD, while the diagnostic yield reported has only been modest [ 15 , 17 ]. With only one sample recurring in our treated and untreated true positives, we could not perform comprehensive analysis to ascertain the diagnostic benefit of using CHX. Expression outlier analysis assumes that loss-of-function (LoF) variants trigger NMD and reduce transcript abundance. Using OUTRIDER, with and without GTEx controls, we observed limited sensitivity and no novel diagnostic findings, consistent with previous reports [ 5 – 7 , 15 , 16 ]. Some outliers corresponded to known heterozygous pathogenic variants, while many others lacked supporting rare variants, likely reflecting pathway-level variation or common polymorphisms [ 6 , 14 ]. Tissue specificity and genome build inconsistencies may further confound interpretation [ 45 , 46 ]. Nevertheless, higher yields reported in other cohorts—such as NIH Undiagnosed Diseases Program fibroblast datasets[ 42 , 43 ] —suggest that expression analysis may perform better in phenotypically diverse or batch sequenced cohorts [ 47 , 48 ]. Allelic imbalance analysis remains underutilized in RNA-seq diagnostics, with only a few studies reporting its clinical application and modest diagnostic yield [ 5 , 14 , 16 ]. While variant level MAE and gene level ANEVA-H have similar sensitivity in our dataset, ANEVA-H produces fewer outliers and has a lower false positive rate. Genes reported as outliers by ANEVA-H frequently harbored structural variants or LoF variants, consistent with prior studies [ 49 ], yet lacked compatible inheritance or phenotypic correlation to make them candidate diagnoses. Imbalance analysis recovered four diagnoses in our true positive benchmarking missed by other tools, supporting its inclusion in RNA analysis pipelines. However, methods such as ANEVA-H remain constrained by haplotype phasing requirements, supporting ensemble approaches for now. This study has several limitations. An optimal study would be prospective, using pre-defined criteria for negative and positive cases. However, our study required a minority subset (splicing and transcription modifying variants) of a set of rare diseases. We felt that a study of retrospective data, with pre-analytically defined case categorization, had a strong potential to generate useful information. That said, our findings would benefit from prospective implementation and replication. Additionally, variants with clinically interpretable mechanisms were selected for assessment of tools, reflective of downstream evaluation of outliers in the clinic. This favors clearly defined aberrant events and likely underestimates tool sensitivity, making discernment of a single tool’s superiority challenging. True negatives were evaluated by presence or absence of outliers to capture time spent assessing outliers for patients without underlying pathogenic aberrant RNA events. Evaluating true negatives by binary outlier classification may overestimate false positive rates, as some outliers may represent genuine biological variation rather than artifacts, and additionally, underestimates predictive ability of tools by lowering MCC scores. Challenges in identifying clear diagnostic candidates in our 74 unresolved cases likely reflect both biological and technical factors. Firstly, unsolved cases may harbor non-RNA-affecting pathogenic mechanisms which would be missed by RNA analytical tools [ 42 , 50 , 51 ]. Secondly, muscle and fibroblast tissues may not be the relevant tissue for capturing pathogenic events in patients with neurogenic or neuromuscular junction phenotypes. Moreover, reliance on exome data in some cases excludes deep intronic variants, and addition of RNA variant calling or genome sequencing could provide more comprehensive detection. Whether to apply RNA-centric versus RNA-complementary approaches remain a central question for diagnostic transcriptomics. Manual investigation of genomic variants in RNA-seq files is inherently complementary. In our benchmarking, RNA-based tools only recovered 28 of 68 known diagnoses with RNA events, and application of tools to undiagnosed cases required rare variant filtration, underscoring the interdependence of transcriptomic and genomic data with current technology. Similar trends have been reported across studies, where diagnostic gains stem largely from RNA evidence supporting previously identified VUSs or genes prioritized due to one pathogenic hit and phenotypic fit [ 5 – 7 , 15 , 16 , 27 ]. RNA-centric pipelines relying solely on RNA variant calls instead may miss pathogenic splice events when the variant is absent in the resulting transcript, while hybrid genome-plus-transcriptome strategies offer the most comprehensive diagnostic coverage [ 7 , 33 , 49 , 52 ]. Given the performance of the current tools, RNA approaches remain complementary to genomic analyses, with manual investigation still a necessary analytical tool. Conclusion With currently available RNA-seq tools, we recommend using ensemble approaches in diagnostic pipelines rather than any single tool. Additionally, we recommend manual inspection of all VUS’s predicted to alter RNA even if no aberrant events are identified using tools. As RNA-seq continues to mature as a diagnostic tool, integrated genome-transcriptome pipelines, improved reference datasets, disease specific transcriptomic signatures and tissue-specific expression ranges will be essential for reducing false positives and enhancing interpretability [ 25 , 53 ]. Increased use of long-read RNA-seq will improve splice detection in difficult to map regions and genes with complex isoforms. Ultimately, such unified analytical frameworks will expedite variant prioritization and accelerate the path to molecular diagnosis. Abbreviations NMD Nonsense Mediated Decay CHX Cycloheximide ASE Allele Specific Expression IGV Integrative Genomics Viewer RNA-Seq RNA-Sequencing VUS Variant of Uncertain Significance LoF Loss of Function NIH National Institutes of Health SNP Single Nucleotide Polymorphism BAM Binary Alignment File AF Allele Frequency VCF Variant Call File MCC Mathews Correlation Coefficient TPM Transcripts Per Million FPR False Positive Rate Declarations Ethics approval and consent to participate Written informed consent and age- appropriate assent for study procedures were obtained by a qualified investigator [protocol 12-N-0095 approved by the National Institute of Neurological Disorders and Stroke, National Institutes of Health IRB]. Consent for publication All participating patients consented to publication of deidentified material Availability of data and materials The datasets generated and/or analyzed during the current study are available in the dbGAP repository phs001272.v2.p1. Competing interests The authors declare that they have no competing interests. Funding The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from dbGaP accession number phs000424.v8.p2 on 09/25/2023. This research was supported by the Intramural Research Program of the National Institutes of Health (NIH) [grant 1ZIANS003129]. The contributions of the NIH author(s) are considered Works of the United States Government. The findings and conclusions presented in this paper are those of the author(s) and do not necessarily reflect the views of the NIH or the U.S. Department of Health and Human Services. Sequencing and analysis were provided by the Broad Institute of MIT and Harvard Center for Mendelian Genomics (Broad CMG) and was funded by the National Human Genome Research Institute, the National Eye Institute, the National Heart, Lung and Blood Institute and grant UM1HG008900 to Daniel MacArthur and Heidi Rehm. This study makes use of data shared through the seqr platform with funding provided by National Institutes of Health grants R01HG009141 and UM1HG008900. PM was supported by the National Institutes of Health under award number R01GM140287 Authors' contributions SS conceptualized the study design, performed all data analysis, and prepared the manuscript. KG ran the ANEVA-H data analysis and contributed manuscript preparation. SG contributed to the study design and manuscript preparation. SD performed initial genetic analysis for each patient used in the study and contributed to clinical interpretation of new candidate diagnoses. VB collected RNA from patients with COL6-RDs. JM ran FRASER and OUTRIDER on the cycloheximide treated RNA-seq samples. YH collected RNA from patient muscle and fibroblast for sequencing. PU contributed to study design and bioinformatic support. VG and BW contributed patient diagnoses used as true positives in the study. RO contribute to clinical evaluation of patients and interpretation of new candidate diagnoses. ARF contribute to clinical evaluation of patients. PM supervised ANEVA-H data analysis. DA supervised SS and conceptualized study design and manuscript preparation. CB supervised SS and conceptualized study design and manuscript preparation. All authors read and approved the final manuscript. Acknowledgements This work utilized the computational resources of the NIH HPC Biowulf cluster ( https://hpc.nih.gov ). With great thanks to referring physicians, patients and their families. We would like to thank Christopher Mendoza, Gilberto (“Mike”) Averion and Kia Brooks for their help and support in coordinating the patient evaluations. References Zambon AA, Muntoni F. Congenital muscular dystrophies: What is new? Neuromuscular Disorders. Elsevier; 2021;31:931–42. https://doi.org/10.1016/J.NMD.2021.07.009 Zambon AA, Pini V, Bosco L, Falzone YM, Munot P, Muntoni F, et al. Early onset hereditary neuronopathies: an update on non-5q motor neuron diseases. Brain [Internet]. Oxford Academic; 2023 [cited 2024 Jan 9];146:806–22. https://doi.org/10.1093/BRAIN/AWAC452 Kang PB, Morrison L, Iannaccone ST, Graham RJ, Bönnemann CG, Rutkowski A, et al. Evidence-based guideline summary: evaluation, diagnosis, and management of congenital muscular dystrophy: Report of the Guideline Development Subcommittee of the American Academy of Neurology and the Practice Issues Review Panel of the American Association of Neuromuscular & Electrodiagnostic Medicine. Neurology [Internet]. Neurology; 2015 [cited 2023 Dec 13];84:1369–78. https://doi.org/10.1212/WNL.0000000000001416 Maggi L, Scoto M, Cirak S, Robb SA, Klein A, Lillis S, et al. Congenital myopathies--clinical features and frequency of individual subtypes diagnosed over a 5-year period in the United Kingdom. Neuromuscul Disord [Internet]. Neuromuscul Disord; 2013 [cited 2024 Jan 8];23:195–205. https://doi.org/10.1016/J.NMD.2013.01.004 Segarra-Casas A, Domínguez-González C, Natera-de Benito D, Kapetanovic S, Hernández-Laín A, Estévez-Arias B, et al. Translating Muscle RNAseq Into the Clinic for the Diagnosis of Muscle Diseases. Ann Clin Transl Neurol [Internet]. John Wiley and Sons Inc; 2025 [cited 2025 Jun 25];0:1–15. https://doi.org/10.1002/ACN3.70078;REQUESTEDJOURNAL:JOURNAL:23289503;JOURNAL:JOURNAL:23289503;WGROUP:STRING:PUBLICATION Luo X, Xiao B, Liang L, Zhang K, Xu T, Liu H, et al. Blood RNA-seq in rare disease diagnostics: a comparative study of cases with and without candidate variants. J Transl Med [Internet]. BioMed Central Ltd; 2025 [cited 2025 Jun 25];23. https://doi.org/10.1186/S12967-025-06609-W, Jaramillo Oquendo C, Wai HA, Rich WI, Bunyan DJ, Thomas NS, Hunt D, et al. Identification of diagnostic candidates in Mendelian disorders using an RNA sequencing-centric approach. Genome Med [Internet]. BioMed Central Ltd; 2024 [cited 2025 Jun 25];16:1–19. https://doi.org/10.1186/S13073-024-01381-W/FIGURES/6 Frésard L, Smail C, Ferraro NM, Teran NA, Li X, Smith KS, et al. Identification of rare-disease genes using blood transcriptome sequencing and large control cohorts. Nat Med. Nature Research; 2019;25:911–9. https://doi.org/10.1038/s41591-019-0457-8 Kremer LS, Bader DM, Mertes C, Kopajtich R, Pichler G, Iuso A, et al. Genetic diagnosis of Mendelian disorders via RNA sequencing. Nat Commun. Nature Publishing Group; 2017;8. https://doi.org/10.1038/ncomms15824 Cummings BB, Marshall JL, Tukiainen T, Lek M, Donkervoort S, Foley AR, et al. Improving genetic diagnosis in Mendelian disease with transcriptome sequencing. Sci Transl Med [Internet]. American Association for the Advancement of Science; 2017 [cited 2024 Jan 16];9. https://doi.org/10.1126/scitranslmed.aal5209 Murdock DR, Dai H, Burrage LC, Rosenfeld JA, Ketkar S, Müller MF, et al. Transcriptome-directed analysis for Mendelian disease diagnosis overcomes limitations of conventional genomic testing. J Clin Invest [Internet]. American Society for Clinical Investigation; 2021 [cited 2024 Jan 16];131. https://doi.org/10.1172/JCI141500 Gonorazky HD, Naumenko S, Ramani AK, Nelakuditi V, Mashouri P, Wang P, et al. Expanding the Boundaries of RNA Sequencing as a Diagnostic Tool for Rare Mendelian Disease. Am J Hum Genet [Internet]. Am J Hum Genet; 2019 [cited 2024 Jan 16];104:466–83. https://doi.org/10.1016/J.AJHG.2019.01.012 Lee H, Huang AY, Wang L kai, Yoon AJ, Renteria G, Eskin A, et al. Diagnostic utility of transcriptome sequencing for rare Mendelian diseases. Genetics in Medicine. Elsevier; 2020;22:490–9. https://doi.org/10.1038/S41436-019-0672-1 Yépez VA, Gusic M, Kopajtich R, Mertes C, Smith NH, Alston CL, et al. Clinical implementation of RNA sequencing for Mendelian disease diagnostics. Genome Medicine 2022 14:1 [Internet]. BioMed Central; 2022 [cited 2022 May 4];14:1–26. https://doi.org/10.1186/S13073-022-01019-9 De Cock L, D’haenens E, Vantomme L, Backers L, Beyens A, Claes KB, et al. Cracking rare disorders: a new minimally invasive RNA-seq protocol. NPJ Genom Med [Internet]. NPJ Genom Med; 2025 [cited 2025 Jun 25];10:45. https://doi.org/10.1038/S41525-025-00502-7, Riquin K, Isidor B, Mercier S, Nizon M, Colin E, Bonneau D, et al. Integrating RNA-Seq into genome sequencing workflow enhances the analysis of structural variants causing neurodevelopmental disorders. J Med Genet [Internet]. BMJ Publishing Group; 2023 [cited 2025 Jun 25];61:47–56. https://doi.org/10.1136/JMG-2023-109263, Dekker J, Schot R, Bongaerts M, de Valk WG, van Veghel-Plandsoen MM, Monfils K, et al. Web-accessible application for identifying pathogenic transcripts with RNA-seq: Increased sensitivity in diagnosis of neurodevelopmental disorders. Am J Hum Genet [Internet]. Cell Press; 2023 [cited 2025 Jun 25];110:251–72. https://doi.org/10.1016/j.ajhg.2022.12.015 Bournazos AM, Riley LG, Bommireddipalli S, Ades L, Akesson LS, Al-Shinnag M, et al. Standardized practices for RNA diagnostics using clinically accessible specimens reclassifies 75% of putative splicing variants. Genetics in Medicine. Elsevier; 2022;24:130–45. https://doi.org/10.1016/J.GIM.2021.09.001 Wang Y, Xie Z, Kutschera E, Adams JI, Kadash-Edmondson KE, Xing Y. rMATS-turbo: an efficient and flexible computational tool for alternative splicing analysis of large-scale RNA-seq data. Nat Protoc [Internet]. Springer Nature; 2024 [cited 2025 Jun 25];19:1083–104. https://doi.org/10.1038/S41596-023-00944-2, Scheller IF, Lutz K, Mertes C, Yépez VA, Gagneur J. Improved detection of aberrant splicing with FRASER 2.0 and the intron Jaccard index. Am J Hum Genet [Internet]. Cell Press; 2023 [cited 2025 Jun 25];110:2056–67. https://doi.org/10.1016/j.ajhg.2023.10.014 Brechtmann F, Mertes C, Matusevičiūtė A, Yépez VA, Avsec Ž, Herzog M, et al. OUTRIDER: A Statistical Method for Detecting Aberrantly Expressed Genes in RNA Sequencing Data. Am J Hum Genet [Internet]. Am J Hum Genet; 2018 [cited 2024 Feb 5];103:907–17. https://doi.org/10.1016/J.AJHG.2018.10.025 Cmero M, Schmidt B, Majewski IJ, Ekert PG, Oshlack A, Davidson NM. MINTIE: identifying novel structural and splice variants in transcriptomes using RNA-seq data. Genome Biol [Internet]. BioMed Central Ltd; 2021 [cited 2024 Feb 5];22:1–26. https://doi.org/10.1186/S13059-021-02507-8/TABLES/2 Mertes C, Scheller IF, Yépez VA, Çelik MH, Liang Y, Kremer LS, et al. Detection of aberrant splicing events in RNA-seq data using FRASER. Nature Communications 2021 12:1 [Internet]. Nature Publishing Group; 2021 [cited 2024 Feb 5];12:1–13. https://doi.org/10.1038/s41467-020-20573-7 Yépez VA, Mertes C, Müller MF, Klaproth-Andrade D, Wachutka L, Frésard L, et al. Detection of aberrant gene expression events in RNA sequencing data. Nature Protocols 2021 16:2 [Internet]. Nature Publishing Group; 2021 [cited 2024 Jan 24];16:1276–96. https://doi.org/10.1038/s41596-020-00462-5 Zhao S, Macakova K, Sinson JC, Dai H, Rosenfeld J, Zapata GE, et al. Clinical validation of RNA sequencing for Mendelian disorder diagnostics. Am J Hum Genet [Internet]. Cell Press; 2025 [cited 2025 Jun 30];112:779–92. https://doi.org/10.1016/j.ajhg.2025.02.006 He W Bin, Xiao WJ, Dai CL, Wang YR, Li XR, Gong F, et al. RNA splicing analysis contributes to reclassifying variants of uncertain significance and improves the diagnosis of monogenic disorders. J Med Genet [Internet]. BMJ Publishing Group; 2022 [cited 2025 Jun 29];59:1010–6. https://doi.org/10.1136/JMEDGENET-2021-108013, Stark JC, Pipko N, Liang Y, Szuto A, Ting Tsoi C, Dickson MA, et al. Clinical applications of and molecular insights from RNA sequencing in a rare disease cohort. Genome Medicine 2025 17:1 [Internet]. BioMed Central; 2025 [cited 2025 Jul 14];17:1–24. https://doi.org/10.1186/S13073-025-01494-W Lewis RG, O’Shea JM, Pizzo L, Wen T, Fulmer ML, Zhao J, et al. RNA sequencing provides functional insights and diagnostic resolution in previously unsolved rare disease cases. BMC Medical Genomics 2025 18:1 [Internet]. BioMed Central; 2025 [cited 2025 Nov 16];18:182-. https://doi.org/10.1186/S12920-025-02227-Z Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics [Internet]. Bioinformatics; 2013 [cited 2025 Jun 25];29:15–21. https://doi.org/10.1093/BIOINFORMATICS/BTS635, Deluca DS, Levin JZ, Sivachenko A, Fennell T, Nazaire MD, Williams C, et al. RNA-SeQC: RNA-seq metrics for quality control and process optimization. Bioinformatics [Internet]. Bioinformatics; 2012 [cited 2024 Feb 19];28:1530–2. https://doi.org/10.1093/BIOINFORMATICS/BTS196 Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol [Internet]. Nat Biotechnol; 2011 [cited 2025 Jun 25];29:24–6. https://doi.org/10.1038/NBT.1754, Thorvaldsdottir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013;14:178–92. https://doi.org/10.1093/bib/bbs017 Montgomery SB, Bernstein JA, Wheeler MT. Toward transcriptomics as a primary tool for rare disease investigation. Cold Spring Harb Mol Case Stud [Internet]. Cold Spring Harbor Laboratory Press; 2022 [cited 2025 Jun 26];8:a006198. https://doi.org/10.1101/MCS.A006198 Jenkinson G, Li YI, Basu S, Cousin MA, Oliver GR, Klee EW. LeafCutterMD: an algorithm for outlier splicing detection in rare diseases. Bioinformatics [Internet]. Bioinformatics; 2020 [cited 2022 Nov 3];36:4609–15. https://doi.org/10.1093/BIOINFORMATICS/BTAA259 Li YI, Knowles DA, Humphrey J, Barbeira AN, Dickinson SP, Im HK, et al. Annotation-free quantification of RNA splicing using LeafCutter. Nature Genetics 2017 50:1 [Internet]. Nature Publishing Group; 2017 [cited 2022 Aug 21];50:151–8. https://doi.org/10.1038/s41588-017-0004-9 Ganapathy KR, Broly M, Silverstein S, Mendoza M, Song E, Kotis B, et al. Improved Identification of Large-effect Rare Genetic Variants using Haplotype Aggregated Allele-specific Expression Data. medRxiv [Internet]. Cold Spring Harbor Laboratory Press; 2025 [cited 2026 Jan 4];2025.12.16.25341855. https://doi.org/10.64898/2025.12.16.25341855 Loh PR, Danecek P, Palamara PF, Fuchsberger C, Reshef YA, Finucane HK, et al. Reference-based phasing using the Haplotype Reference Consortium panel. Nat Genet [Internet]. Nature Publishing Group; 2016 [cited 2025 Jun 25];48:1443–8. https://doi.org/10.1038/NG.3679, Byrska-Bishop M, Evani US, Zhao X, Basile AO, Abel HJ, Regier AA, et al. High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios. Cell [Internet]. Elsevier B.V.; 2022 [cited 2025 Jun 25];185:3426-3440.e19. https://doi.org/10.1016/J.CELL.2022.08.004, Castel SE, Mohammadi P, Chung WK, Shen Y, Lappalainen T. Rare variant phasing and haplotypic expression from RNA sequencing with phASER. Nat Commun [Internet]. Nature Publishing Group; 2016 [cited 2025 Jun 25];7:12817. https://doi.org/10.1038/NCOMMS12817 Mohammadi P, Castel SE, Cummings BB, Einson J, Sousa C, Hoffman P, et al. Genetic regulatory variation in populations informs transcriptome analysis in rare disease. Science (1979) [Internet]. American Association for the Advancement of Science; 2019 [cited 2025 Jun 25];366:351–6. https://doi.org/10.1126/SCIENCE.AAY0256/SUPPL_FILE/AAY0256_TABLES5.XLSX McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The Ensembl Variant Effect Predictor. Genome Biol [Internet]. BioMed Central Ltd.; 2016 [cited 2025 Sep 14];17:1–14. https://doi.org/10.1186/S13059-016-0974-4/TABLES/8 Jaganathan K, Kyriazopoulou Panagiotopoulou S, McRae JF, Darbandi SF, Knowles D, Li YI, et al. Predicting Splicing from Primary Sequence with Deep Learning. Cell [Internet]. Cell Press; 2019 [cited 2025 Nov 4];176:535-548.e24. https://doi.org/10.1016/j.cell.2018.12.015 Pais LS, Snow H, Weisburd B, Zhang S, Baxter SM, DiTroia S, et al. seqr: A web-based analysis and collaboration tool for rare disease genomics. Hum Mutat [Internet]. John Wiley and Sons Inc; 2022 [cited 2025 Jun 25];43:698–707. https://doi.org/10.1002/HUMU.24366, Sobreira N, Schiettecatte F, Valle D, Hamosh A. GeneMatcher: A Matching Tool for Connecting Investigators with an Interest in the Same Gene. Hum Mutat [Internet]. John Wiley and Sons Inc.; 2015 [cited 2025 Sep 14];36:928–30. https://doi.org/10.1002/HUMU.22844, Ungar RA, Goddard PC, Jensen TD, Degalez F, Smith KS, Jin CA, et al. Impact of genome build on RNA-seq interpretation and diagnostics. Am J Hum Genet [Internet]. Cell Press; 2024 [cited 2025 Jul 1];111:1282–300. https://doi.org/10.1016/j.ajhg.2024.05.005 Karam R, Wengrod J, Gardner LB, Wilkinson MF. Regulation of nonsense-mediated mRNA decay: Implications for physiology and disease. Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms. 2013;1829:624–33. https://doi.org/10.1016/j.bbagrm.2013.03.002 Cassini T, Silverstein S, Behan M, Tifft CJ, Malicdan MC, Adams DR, et al. Mitochondrial trifunctional protein deficiency caused by a deep intronic deletion leading to aberrant splicing. JIMD Rep [Internet]. John Wiley and Sons Inc; 2025 [cited 2025 Jul 1];66. https://doi.org/10.1002/JMD2.12459, Silverstein S, Cassini T, Fu J, Pusey B, Macnamara E, Frost FG, et al. RNA sequencing driven diagnosis expands the phenotypic spectrum of NBAS deficiency. Mol Genet Metab [Internet]. Academic Press Inc.; 2025 [cited 2025 Jul 1];145. https://doi.org/10.1016/j.ymgme.2025.109105 Ferraro NM, Strober BJ, Einson J, Abell NS, Aguet F, Barbeira AN, et al. Transcriptomic signatures across human tissues identify functional rare genetic variation. Science. 2020;369. https://doi.org/10.1126/science.aaz5900 Soemedi R, Cygan KJ, Rhine CL, Wang J, Bulacan C, Yang J, et al. Pathogenic variants that alter protein code often disrupt splicing. Nat Genet [Internet]. Nat Genet; 2017 [cited 2025 Nov 4];49:848–55. https://doi.org/10.1038/NG.3837 Chong R, Insigne KD, Yao D, Burghard CP, Wang J, Hsiao YHE, et al. A Multiplexed Assay for Exon Recognition Reveals that an Unappreciated Fraction of Rare Genetic Variants Cause Large-Effect Splicing Disruptions. Mol Cell [Internet]. Cell Press; 2019 [cited 2025 Nov 4];73:183-194.e8. https://doi.org/10.1016/j.molcel.2018.10.037 Guo H, Urban AE, Wong WH. Prioritizing disease-related rare variants by integrating gene expression data. PLoS Genet [Internet]. Public Library of Science; 2024 [cited 2025 Sep 1];20:e1011412. https://doi.org/10.1371/JOURNAL.PGEN.1011412 Arriaga MT, Mendez R, Ungar RA, Bonner DE, Matalon DR, Lemire G, et al. Transcriptome-wide outlier approach identifies individuals with minor spliceopathies. medRxiv [Internet]. medRxiv; 2025 [cited 2025 Sep 15]; https://doi.org/10.1101/2025.01.02.24318941 Table Table 1: Results of true positive expression analysis across Muscle, Fibroblast and Cycloheximide (CHX)- Treated Fibroblasts sets. Muscle samples were run with and without an additional 100 healthy samples from GTEx. Cohort Muscle Muscle + GTEx Fibroblast CHX Sample size 71 171 37 79 Recall 1/14 (7%) 1/14 (7%) 0/1 (0%) NA Avg Outliers/sample 3.08 4.02 0.8 1.21 Additional Declarations No competing interests reported. Supplementary Files 120125SupplementalMethods.docx 120125SupplementalTables.xlsx Supplemental Tables Table S1: Details on sequencing and library preparation across cohorts Table S2: GTEx accession numbers for samples used in this study Table S3: Muscle splicing true positives Table S4: Fibroblast splicing true positives Table S5: CHX-treated fibroblast splicing true positives Table S6: Variant filtration of detected true positives in muscle Table S7: Expression true positives Table S8: Allelic imbalance true positives Table S9: Confusion matrices and performance metrics Table S10: Splicing candidate variants in 29 undiagnosed muscle samples 100626SupplementalFigures.pdf Supplemental Figures Figure S1: PCA analysis of cohorts Figure S2: Breakdown of biopsy locations for muscle RNA-seq Figure S3: FRASER parameters and preliminary tests Figure S4: Outliers reported for splicing tools in fibroblast and cycloheximide treated fibroblast samples Figure S5: Specificity rates per tool per tissue origin Figure S6: Spurious find in cycloheximide treated fibroblast OUTRIDER analysis Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-8622051","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":581085250,"identity":"e9e75d5f-86dd-48bd-8e63-bc14b09d5161","order_by":0,"name":"Sarah Silverstein","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABCklEQVRIie2RMUsDMRTHXwnklj92PWmpXyGQQcRiv0rKQVxEOjro0UkX4dZKv0wk0C4R15NOInRzkFtOHDSndqrRuonkB4G8P/nxHi9EkcgfhfmDrc8ClIzJbqTw96vyCsxmCq0UolR9JCF2Lm7uqxHlXZ5ms6fnut/d3X68tiM66wnztSLcoexMyIKnOptCaexNj5Wd0FwGFdLUARlwOMlIWYjFkbCg2fBqHBisWLIXUO6V26pVq1eIO/e9QqXmvgsDTy4ZQRmIEo1yOiwChiiXfL95w5NzyaAzCNcMJoxshzZWaLbAST5oM/bQqvsHAzF3svJJjweU1RLWkh9/c53890okEon8V94AlUVMIlGHok8AAAAASUVORK5CYII=","orcid":"","institution":"National Institutes of Health","correspondingAuthor":true,"prefix":"","firstName":"Sarah","middleName":"","lastName":"Silverstein","suffix":""},{"id":581085251,"identity":"b4ad7cdd-8f34-4be8-8bfb-1c43317f6f4e","order_by":1,"name":"Kaushik Ganapathy","email":"","orcid":"","institution":"Scripps Research Institute","correspondingAuthor":false,"prefix":"","firstName":"Kaushik","middleName":"","lastName":"Ganapathy","suffix":""},{"id":581085252,"identity":"9fe4f6c6-4d70-426d-a7e5-85b9ec2cd807","order_by":2,"name":"Sandra Donkervoort","email":"","orcid":"","institution":"National Institutes of Health","correspondingAuthor":false,"prefix":"","firstName":"Sandra","middleName":"","lastName":"Donkervoort","suffix":""},{"id":581085254,"identity":"345526fb-820e-4779-a314-0bbc905f921c","order_by":3,"name":"Veronique Bolduc","email":"","orcid":"","institution":"National Institutes of Health","correspondingAuthor":false,"prefix":"","firstName":"Veronique","middleName":"","lastName":"Bolduc","suffix":""},{"id":581085255,"identity":"ee252b75-c045-4b9a-9a93-446246bd2fca","order_by":4,"name":"Ying Hu","email":"","orcid":"","institution":"National Institutes of Health","correspondingAuthor":false,"prefix":"","firstName":"Ying","middleName":"","lastName":"Hu","suffix":""},{"id":581085256,"identity":"5718a12b-e250-45e1-bf00-7b3c960312cf","order_by":5,"name":"Justin Moy","email":"","orcid":"","institution":"Boston University","correspondingAuthor":false,"prefix":"","firstName":"Justin","middleName":"","lastName":"Moy","suffix":""},{"id":581085258,"identity":"f8e1e15d-489d-4a46-826b-93e466a202f2","order_by":6,"name":"Prech Uapinyoying","email":"","orcid":"","institution":"National Institutes of Health","correspondingAuthor":false,"prefix":"","firstName":"Prech","middleName":"","lastName":"Uapinyoying","suffix":""},{"id":581085261,"identity":"fa501947-7ea5-4009-ba51-1c554b6fe31d","order_by":7,"name":"Svetlana Gorokhova","email":"","orcid":"","institution":"Aix-Marseille University","correspondingAuthor":false,"prefix":"","firstName":"Svetlana","middleName":"","lastName":"Gorokhova","suffix":""},{"id":581085263,"identity":"e5c473b5-b498-4787-b4b8-2a0a088aa1b5","order_by":8,"name":"Vijay Ganesh","email":"","orcid":"","institution":"Broad Institute","correspondingAuthor":false,"prefix":"","firstName":"Vijay","middleName":"","lastName":"Ganesh","suffix":""},{"id":581085265,"identity":"23a50c4e-0a12-4ec0-b053-593f8f386094","order_by":9,"name":"Ben Weisburd","email":"","orcid":"","institution":"Broad Institute","correspondingAuthor":false,"prefix":"","firstName":"Ben","middleName":"","lastName":"Weisburd","suffix":""},{"id":581085266,"identity":"e72dee27-670b-4a57-bdc1-08778431d00e","order_by":10,"name":"Rotem OrBach","email":"","orcid":"","institution":"National Institutes of Health","correspondingAuthor":false,"prefix":"","firstName":"Rotem","middleName":"","lastName":"OrBach","suffix":""},{"id":581085267,"identity":"b6f45a8c-f829-4e45-a7fb-8046763cae7a","order_by":11,"name":"A. Reghan Foley","email":"","orcid":"","institution":"National Institutes of Health","correspondingAuthor":false,"prefix":"","firstName":"A.","middleName":"Reghan","lastName":"Foley","suffix":""},{"id":581085268,"identity":"a07e798a-31f9-4213-8a14-51e452d72d5e","order_by":12,"name":"Pejman Mohammadi","email":"","orcid":"","institution":"Seattle Children's Research Institute","correspondingAuthor":false,"prefix":"","firstName":"Pejman","middleName":"","lastName":"Mohammadi","suffix":""},{"id":581085269,"identity":"940dcd72-7581-4adc-ad7e-2961590fd2c4","order_by":13,"name":"David Adams","email":"","orcid":"","institution":"National Human Genome Research Institute","correspondingAuthor":false,"prefix":"","firstName":"David","middleName":"","lastName":"Adams","suffix":""},{"id":581085270,"identity":"03deb455-c1af-4982-b9bb-3a5a6b310bb7","order_by":14,"name":"Carsten Bonnemann","email":"","orcid":"","institution":"National Institutes of Health","correspondingAuthor":false,"prefix":"","firstName":"Carsten","middleName":"","lastName":"Bonnemann","suffix":""}],"badges":[],"createdAt":"2026-01-16 20:38:14","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-8622051/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-8622051/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":101365240,"identity":"38739f8d-0908-43c2-9be4-3c26842e626c","added_by":"auto","created_at":"2026-01-29 00:54:04","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1516690,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eTruth cohorts included for analysis, workflow design and overall results.\u003c/em\u003e A) Truth set comprises of 97 samples, split into true positive and true negative sets. True positives and negatives were further categorized by tissue of origin and by the number of variants or genes contributed to each analysis type B) Benchmarking was performed on the eight tools with diverse methodolgies covering three RNA analysis approaches: aberrant splicing, allelic imbalance and expression outlier detection (ASE = allele specific expression). C) A total of 76 true positive samples from 68 probands with pathogenic variants leading to confirmed RNA events. Probands have causative variants leading to more than one aberrant RNA event type. Each proband is reported here by genetic diagnosis along with which analysis type the event was considered a true positive for. Red dots for undetected event, green dots for correct detection and grey dots when the sample was not tested.\u003c/p\u003e","description":"","filename":"OnlineFig1.png","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/6885f35dcf2ce81d5b27a110.png"},{"id":101365244,"identity":"8a7ca448-99da-4524-bab6-ebfc697e623a","added_by":"auto","created_at":"2026-01-29 00:54:06","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":2765554,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eResults of true positive aberrant splice analysis\u003c/em\u003e. A) Recall rates of true positive splice variants for rMATS, MINTIE, LeafcutterMD, FRASER and FRASER2.0, broken down by tissue of origin. FRASER2.0 was only tested in the muscle set. For CHX set, only LeafcutterMD and FRASER were evaluated. Muscle = muscle-derived RNA-seq (16 true positive variants), Fibroblast = fibroblast derived RNA-seq (7 true positive variants), and CHX = cycloheximide-treated fibroblast derived RNA-seq (45 true positive variants) B) True positive splice variants create more than one aberrant splice event. Splice events were categorized as follows: donor event (CD = cryptic donor, ND = novel donor), acceptor event (CA = cryptic acceptor, NA= novel acceptor), exon skip (ES), pseudoexon (PE) and retained intron (RI). To evaluate any splice category preference among tools, each splice event from muscle true positive set is reported clustered by category with a dot if correctly detected by corresponding tool. An ensemble approach detects the most events. C) Violin plots showing the number of outliers reported for each sample by each tool in the muscle set, colored by library preparation differences. D) Variant based filtration strategy: outliers from all tools were merged into an ensemble file and intersected with rare DNA variants to reduce outlier burden, reporting outliers associated with rare variants only. E) Out of 16 true positive splice variants in the muscle set, 11 variants correctly detected (arising from 10 probands) were evaluated with variant filtration. Outliers were aggregated across all splice tools, excluding LeafcutterMD (average 13745.5 outliers/sample). Two allele frequencies were evaluated to assess maximum and minimum final outliers after variant based filtration. A lenient threshold representing the maximum variant allele frequency used for rare disease analysis (AF\u0026lt;0.01) with average outlier-associated variants after filtration 1220.6. Variants were then segregated, searched for autosomal dominant (AD) or autosomal recessive (AR) inheritance patterns and restricted to a stringent allele frequency threshold commonly applied in rare disease (AF\u0026lt;0.0005), reducing outliers to an average of 144.2 and 154.2, respectively. A final filtration by SpliceAI score \u0026lt; 0.1 reduced outliers to 2.9 and 1.7, respectively. Across 10 probands, 9 cases were correctly included in the final outlier list.\u003c/p\u003e","description":"","filename":"OnlineFig2.png","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/12e8e6705720936a03594009.png"},{"id":101365242,"identity":"828bf1d3-34b7-4c67-b7c3-1ad815411061","added_by":"auto","created_at":"2026-01-29 00:54:05","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":593494,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eResults of true positive allelic imbalance analysis\u003c/em\u003e. A)Recall rates of 8 true positive genes for both tested tools, ANEVA-H and MAE. True positives were only tested in muscle derived RNA-seq samples B) Number of outliers reported for each sample per tool (59 total samples run).\u003c/p\u003e","description":"","filename":"OnlineFig3.png","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/88ee111aecb8c3d4c5297881.png"},{"id":101365243,"identity":"5685e2db-785e-45d9-b1bc-fec6dc29b53b","added_by":"auto","created_at":"2026-01-29 00:54:05","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":2468849,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003ePrecision, accuracy and predictive power\u003c/em\u003e. A) True negative samples (n=21) were evaluated across all tools by presence or absence of reported outliers. Results were aggregated across tissue sources to quanity the false positive rate for each tool tested, with 100% false positive rate reported when all true negative samples had outliers. Bar plots show the false positive rate for each tool, categorized by analysis subtypes splicing (rMATS = 9, MINTIE =9, FRASER = 21, FRASER2.0 = 6, LeafCutterMD=21), expression (OUTRIDER =9, OUTRIDER+GTEx =6) and allelic imbalance (ANEVA-H =5, MAE =5). B) True positve samples were similarly aggegated across tissue sources for each tool, allowing for determination of true positives, false positives, true negatives and false negatives for each tool (Table S9). Accuracy, precision, recall, F1 score and MCC were calculated. Radar plot show metrics for splice tools C) for allelic imbalance tools and D) for expression tools.\u003c/p\u003e","description":"","filename":"OnlineFig4.png","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/a8b79ba91c124610cd4878bb.png"},{"id":101365217,"identity":"49d66695-14ff-4d03-b1e7-b8c68dc674ba","added_by":"auto","created_at":"2026-01-29 00:54:03","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":711560,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eApplication of RNA analysis to undiagnosed samples\u003c/em\u003e. A) 74 undiagnosed samples were evaluated by RNA analytical tools to identify candidate variants. Figure shows the numbers split by tissue source and list of analytical tools applied to each set of undiagnosed samples. B) For splice analysis, an ensemble approach was used along with filtration by rare DNA variants (AF \u0026lt; 0.0005). Bar plot shows subsequent steps of the filtration method in reducing the number of outliers to evaluate. C) 129 final outliers were manually evaluated across 28 patients, resulting in 9 final candidates submitted to geneMatcher.\u003c/p\u003e","description":"","filename":"OnlineFig5.png","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/e30660ea11575e64d5b8a68b.png"},{"id":102397169,"identity":"6ef2f82f-74e2-47df-9f5f-38554d2b896e","added_by":"auto","created_at":"2026-02-11 10:06:04","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3042050,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/8d91b252-1913-48d8-91f0-7eb815d3874d.pdf"},{"id":101365198,"identity":"e4842b2c-4cf2-41a4-ab47-7ece0c8339d2","added_by":"auto","created_at":"2026-01-29 00:54:00","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":29902,"visible":true,"origin":"","legend":"","description":"","filename":"120125SupplementalMethods.docx","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/5423b3a14b977cd74ec57f30.docx"},{"id":101365265,"identity":"d2f6124c-2ddc-495f-a131-875387aa9d25","added_by":"auto","created_at":"2026-01-29 00:54:08","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":41367,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eSupplemental Tables\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTable S1: Details on sequencing and library preparation across cohorts\u003c/p\u003e\n\u003cp\u003eTable S2: GTEx accession numbers for samples used in this study\u003c/p\u003e\n\u003cp\u003eTable S3: Muscle splicing true positives\u003c/p\u003e\n\u003cp\u003eTable S4: Fibroblast splicing true positives\u003c/p\u003e\n\u003cp\u003eTable S5: CHX-treated fibroblast splicing true positives\u003c/p\u003e\n\u003cp\u003eTable S6: Variant filtration of detected true positives in muscle\u003c/p\u003e\n\u003cp\u003eTable S7: Expression true positives\u003c/p\u003e\n\u003cp\u003eTable S8: Allelic imbalance true positives\u003c/p\u003e\n\u003cp\u003eTable S9: Confusion matrices and performance metrics\u003c/p\u003e\n\u003cp\u003eTable S10: Splicing candidate variants in 29 undiagnosed muscle samples\u003c/p\u003e","description":"","filename":"120125SupplementalTables.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/2a92d86aded0d260e8d0255b.xlsx"},{"id":101365266,"identity":"10e528a1-d929-4330-8423-10a318a3178f","added_by":"auto","created_at":"2026-01-29 00:54:08","extension":"pdf","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":2611676,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003eSupplemental Figures\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eFigure S1: PCA analysis of cohorts\u003c/p\u003e\n\u003cp\u003eFigure S2: Breakdown of biopsy locations for muscle RNA-seq\u003c/p\u003e\n\u003cp\u003eFigure S3: FRASER parameters and preliminary tests\u003c/p\u003e\n\u003cp\u003eFigure S4: Outliers reported for splicing tools in fibroblast and cycloheximide treated fibroblast samples\u003c/p\u003e\n\u003cp\u003eFigure S5: Specificity rates per tool per tissue origin\u003c/p\u003e\n\u003cp\u003eFigure S6: Spurious find in cycloheximide treated fibroblast OUTRIDER analysis\u003c/p\u003e","description":"","filename":"100626SupplementalFigures.pdf","url":"https://assets-eu.researchsquare.com/files/rs-8622051/v1/788931eb0bc2f602368550b3.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Benchmarking RNA-seq Tools for Real-World Diagnostic Applications","fulltext":[{"header":"Background","content":"\u003cp\u003ePediatric neuromuscular disorders are characterized by broad genetic and clinical heterogeneity, with approximately 30\u0026ndash;60% of cases remaining molecularly undiagnosed following exome or genome sequencing[\u003cspan additionalcitationids=\"CR2 CR3\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e]. Recent reports suggest that incorporating RNA-seq data derived from clinically accessible tissues such as whole blood, peripheral blood mononuclear cells, fibroblasts, urine, or muscle can increase diagnostic yield by an additional 2.7\u0026ndash;32%[\u003cspan additionalcitationids=\"CR6 CR7 CR8 CR9 CR10 CR11 CR12 CR13 CR14 CR15 CR16\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. In a typical clinical-diagnostic bioinformatics workflow of DNA sequencing, variants are prioritized based on phenotypically-driven analysis, possibly identifying candidate variants with effect on splicing or gene/protein expression. If RNA-seq data is obtained, splicing and expression of the genes with these variants are then manually inspected in the RNA-seq BAM file data. This process is time-consuming and, despite the availability of published guidelines, remains heavily reliant on the expertise of individual analysts and the availability of appropriate control and tissue specific datasets[\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eRecently, open-source computational tools have been developed to automate the analysis of transcriptomic data in the context of rare disease[\u003cspan additionalcitationids=\"CR20 CR21 CR22 CR23\" citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. These tools are designed to identify aberrant RNA events\u0026mdash;such as abnormal splicing, aberrant gene expression, or allelic imbalance\u0026mdash;in a systematic and agnostic manner. Unlike manual methods, these tools employ an RNA-first strategy, enabling transcriptome-wide analysis independent of prior DNA variant prioritization. This approach has the potential to identify pathogenic events that may be overlooked during manual review. Applying open-source strategies to RNA-seq datasets derived from clinically accessible tissues has resulted in modest diagnostic yield, with the highest yield in cases where initial DNA sequencing identified a variant of uncertain significance (VUS) predicted to alter splicing and prioritized [\u003cspan additionalcitationids=\"CR6\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e, \u003cspan additionalcitationids=\"CR26 CR27\" citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]. In cases without any diagnostic leads, RNA-first strategies have resulted in limited diagnostic yield [\u003cspan additionalcitationids=\"CR6\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e, \u003cspan additionalcitationids=\"CR26\" citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. To date, comprehensive benchmarking of these tools using previously diagnosed samples harboring pathogenic variants leading to known aberrant RNA events has yet to be undertaken. Furthermore, prior studies have typically sequenced large batches of samples and analyzed them collectively to mitigate batch effects. This strategy does not reflect real-world clinical constraints, where disease-relevant tissue is often limited, requiring evaluation of open-source tools using diverse approaches.\u003c/p\u003e \u003cp\u003eWe assessed a range of RNA-seq analysis tools to identify those best suited for our diagnostic pipeline and to assess their limitations. The most effective strategy was then applied to a set of clinically well-characterized, undiagnosed cases to enable candidate variant discovery. We also evaluated the capacity of RNA-seq to function as a standalone diagnostic tool, versus being utilized as an adjunct to DNA-based approaches.\u003c/p\u003e"},{"header":"Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003ePatient Recruitment\u003c/h2\u003e \u003cp\u003ePatients were referred by their neurologist or geneticist. Written informed consent and age- appropriate assent for study procedures were obtained by a qualified investigator [protocol 12-N-0095 approved by the National Institute of Neurological Disorders and Stroke, National Institutes of Health IRB]. Medical history was obtained, including clinical evaluations and muscle biopsy, were performed as part of the standard diagnostic examination. Samples for research- based testing, including skin, blood and muscle biopsy samples, were obtained using standard procedures. Patients with confirmed genetic diagnosis were included as true positive or true negative benchmarking samples, and patients without candidates after evaluation were included for our undiagnosed samples analysis (Supplemental Methods).\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eRNA-seq Alignment\u003c/h3\u003e\n\u003cp\u003eQuality control of FASTQ files were assessed with FastQC (v0.12.1). BAM files were generated using the GTEXv10 pipeline (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/broadinstitute/gtex-pipeline\u003c/span\u003e\u003cspan address=\"https://github.com/broadinstitute/gtex-pipeline\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) using reference genome GRCh38 (Gencode, version 39). Briefly, FASTQ files were aligned using STAR (v2.7.10a) and optical duplicates marked by Picard Mark Duplicates (v3.3.0)[\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e]. Alignment metrics were assessed with RNA-SeQC (v2.4.2)[\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e]. Principal component analysis was performed using Transcript per million (TPM) values generated by RNA-SeQC and the prcomp function in R package stats (v4.4.0). Samples not clustering with other samples annotated to be derived from the same tissue source (either muscle, fibroblast or cycloheximide treated fibroblast) were removed (Fig \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e\n\u003ch3\u003eTrue Positive Definition\u003c/h3\u003e\n\u003cp\u003eDiagnosed samples that passed quality control and contained diagnostic aberrant RNA events were manually reviewed in the Integrative Genomics Viewer (IGV)[\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e, \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e]. Each sample contributed one true positive gene for evaluation. Each true positive gene could be included as a true positive across more than one analysis type (\u003cem\u003ee.g\u003c/em\u003e. splicing and expression) based on underlying mechanism. Splicing true positives were defined as containing\u0026thinsp;\u0026ge;\u0026thinsp;3 reads supporting a novel splice junction created by the pathogenic variant. Expression true positives were defined by a confirmed loss-of-function mechanism with normalized transcript abundance\u0026thinsp;\u0026le;\u0026thinsp;50% or \u0026ge;\u0026thinsp;100% of the corresponding normalized GTEx TPM. Genes were normalized to the geometric mean of \u003cem\u003eHPRT1\u003c/em\u003e, \u003cem\u003ePPIA\u003c/em\u003e, \u003cem\u003eYWHAZ\u003c/em\u003e and \u003cem\u003eTBP\u003c/em\u003e. Allelic imbalance true positives were defined by monoallelic loss-of-function events supported by single nucleotide polymorphism (SNP) analysis demonstrating a marked allelic ratio\u0026thinsp;\u0026ge;\u0026thinsp;70:30 at the site of the rare variant.\u003c/p\u003e\n\u003ch3\u003eTrue Negative Definition\u003c/h3\u003e\n\u003cp\u003ePatients harboring pathogenic variants with no direct effect on splicing or expression (missense, in frame indels, terminal gene frameshifts) were considered as \u0026ldquo;true negatives\u0026rdquo; for RNA anomaly. RNA-seq data for these patients were manually reviewed in IGV to confirm that splicing was canonical and that no expression changes or allelic imbalance were present.\u003c/p\u003e\n\u003ch3\u003eComputational Tools\u003c/h3\u003e\n\u003cp\u003eAll tools were run on Biowulf (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://hpc.nih.gov\u003c/span\u003e\u003cspan address=\"https://hpc.nih.gov\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cem\u003e).\u003c/em\u003e Current RNA-seq tools for rare disease use two main strategies: outlier-based (one-vs-all) and event-based (case-vs-control) approaches [\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e]. Outlier-based methods assume that each rare disease is unique within a heterogeneous cohort, an assumption that holds for some muscle disorder cohorts. However, it was not a valid assumption for our dataset, which includes recurrent variants and genes. Additionally, small-batch sequencing within the cohort likely introduced variability across samples. Event-based tools that compare cases to matched controls were considered as a strategy to mitigate these issues. We attempted to evaluate both case-versus-control and one-versus-all strategies when possible; however, this was only possible for splicing analysis in practice.\u003c/p\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eOne vs All Splice Tools (Outlier based approaches)\u003c/h2\u003e \u003cp\u003eThree popular tools\u0026mdash;FRASER, FRASER2.0, and LeafcutterMD [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e, \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e, \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e] \u0026mdash;were selected for evaluation. Since these tools require\u0026thinsp;\u0026ge;\u0026thinsp;30 samples, we created sample groups by tissue of origin, mixing diagnosed and undiagnosed RNA-seq samples from the same tissue: muscle (n\u0026thinsp;=\u0026thinsp;71), fibroblasts (n\u0026thinsp;=\u0026thinsp;37), and cycloheximide-treated fibroblasts (CHX, n\u0026thinsp;=\u0026thinsp;79). Muscle samples were sequenced in multiple small batches over 14 years using different library preparations (Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e), primarily from the quadriceps (Fig \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). Although FRASER recommends analyzing stranded and unstranded libraries separately, preliminary analysis showed higher true positive recall when both library preparation methods were combined (Fig \u003cspan refid=\"MOESM3\" class=\"InternalRef\"\u003eS3\u003c/span\u003eE-F). Therefore, all three tools were run with mixed library preparation.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFRASER was executed via the DROP pipeline (v1.2.2), using default parameters and a false discovery rate (FDR) of 5% (Fig \u003cspan refid=\"MOESM3\" class=\"InternalRef\"\u003eS3\u003c/span\u003eA-C) [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. FRASER2.0 was similarly run via the DROP pipeline (v1.4.0) with standard settings and an FDR of 5% (Fig \u003cspan refid=\"MOESM3\" class=\"InternalRef\"\u003eS3\u003c/span\u003eD) [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. LeafcutterMD (v0.2.9) was applied with preprocessing of BAM files to generate intron clusters using Leafcutter, with recommended parameters (minimum 50 splice junctions per cluster and a maximum intron distance of 500 kb)[\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e, \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e]. Outlier detection was restricted to groups with at least 20 reads and a minimum intron cluster size of 50. Resulting junctions were annotated with Ensembl gene IDs using Bedops (v2.4.41), followed by conversion to gene symbols utilizing the annotationDBI package in R. Correction for multiple hypothesis testing was performed using a FDR of 5%.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eCase vs Control Splice Tools (Events-based approaches)\u003c/h3\u003e\n\u003cp\u003eTwo commonly used tools, rMATS-Turbo and MINTIE, were selected[\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. Each patient sample was run against a set of ten healthy control samples. We used ten samples for the control cohort based on MINTIE recommendations. Control samples were selected from the GTEx database to match the tissue source when available (Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). Muscle control samples were filtered for RNA Integrity Number (RIN)\u0026thinsp;\u0026gt;\u0026thinsp;8, ischemic time\u0026thinsp;\u0026lt;\u0026thinsp;1000 minutes, and absence of autolysis. Fibroblast controls were selected based on tissue type \"skin\" and tissue type detail \"cells \u0026ndash; cultured fibroblasts,\" further filtered for RIN\u0026thinsp;=\u0026thinsp;10 and ischemic time\u0026thinsp;\u0026lt;\u0026thinsp;400 minutes. Due to the lack of available disease-free controls for CHX samples, case vs controls approaches were not tested.\u003c/p\u003e \u003cp\u003eMINTIE (v0.4.1) was run using FASTQ files from each patient sample against the ten selected GTEx controls, with recommended parameters and additional settings (splice_motif_mismatch\u0026thinsp;=\u0026thinsp;3, FDR\u0026thinsp;=\u0026thinsp;5%)[\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. Output was filtered to include only splice-related outliers, excluding large deletions, insertions, and other structural variants. rMATS-Turbo (v4.1.2) was executed on BAM files, with each patient sample compared to ten controls[\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]. Standard parameters were used, with additional flags --novelSS and --allow-clipping. The maximum exon length for novel splice isoforms was set to 700 bp (-mel 700) to accommodate detection of large pseudoexons based on the largest in our true positive set. Correction for multiple hypothesis testing was performed using an FDR of 5%. All five MATS output files were included for assessment of true positive events, regardless of whether the detected aberrant splice category matched the expected event (e.g., intronic cryptic donor activation detected as intron retention).\u003c/p\u003e\n\u003ch3\u003eExpression Tools\u003c/h3\u003e\n\u003cp\u003eOUTRIDER was run on all three cohorts as previously described in a one vs all manner via the DROP pipeline (v1.2.2) with recommended parameters which included an FDR of 1%[\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e]. Since OUTRIDER is designed as a one vs all approach, we performed a second run for the muscle patient samples with an addition of 100 control muscle samples downloaded from GTEx to test if the influx of healthy samples improves detection (Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e).\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eAllelic Imbalance Tools\u003c/h2\u003e \u003cp\u003eThe approach to detecting allelic imbalance in transcriptome data relies on heterozygous presence of variants in DNA. Thus, we restricted our test set to muscle samples where exome or genome variant call files (VCFs) were readily available, for a subset cohort of 59 samples. Two tools were selected, MAE module in the DROP pipeline [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e] which uses variant level allele specific expression (ASE) and ANEVA-h, which uses gene level ASE[\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e]. The MAE module in DROP (v1.4.0) was run using standard parameters except modification of the allelic balance threshold to 70 instead of 80, reflecting our true positive dataset[\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e]. ANEVA-h required preprocessing of VCF and BAM files as follows. Haplotype-resolved, gene-level ASE data were generated for each case to enable regulatory outlier detection using ANEVA-h. RNA-seq libraries were processed and aligned to the GRCh38 reference genome as previously described. Corresponding genotype data were obtained from matched genome sequencing (GS) and phased using EAGLE (v2.4.1), leveraging the 1000 Genomes Project Phase 3 (30x) high-coverage reference panel for population-based phasing[\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e, \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e]. Phased genotypes and aligned RNA-seq BAM files were used as input to phASER (v1.1.1), which was run following best practices for haplotype-resolved ASE quantification [\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e]. Gene-level ASE values were computed using the phaser_gene_ae.py module, with heterozygous exonic SNPs aggregated by haplotype using GENCODE v39 transcript annotations [\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e]. ANEVA-h was then applied using the ANEVADOT_test function from the ANEVA-DOT R package, with default parameters [\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e]. Analyses were performed using reference VG estimates derived from 706 skeletal muscle RNA-seq samples from GTEx v8, ensuring tissue-specific calibration ( \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://zenodo.org/records/15226130\u003c/span\u003e\u003cspan address=\"https://zenodo.org/records/15226130\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://zenodo.org/records/15226575\u003c/span\u003e\u003cspan address=\"https://zenodo.org/records/15226575\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). A minimum coverage threshold of 10 reads per gene was applied for inclusion. Multiple testing correction was performed using the Benjamini\u0026ndash;Hochberg FDR procedure, and genes with FDR-adjusted p-values\u0026thinsp;\u0026lt;\u0026thinsp;0.05 were designated as significant ASE outliers.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eTrue Positive and Negative Detection\u003c/h2\u003e \u003cp\u003eTool performance was evaluated based on detection of known aberrant RNA events. For splicing, a correctly detected true positive was defined by an exact coordinate match to the validated splice event. For expression, the true positive gene had to appear among reported outliers. For allelic imbalance, MAE required inclusion of the correct rare variant, while ANEVA-H required inclusion of the corresponding gene. True negatives were assessed by the absence of reported outliers, with correct detection defined as no outliers identified in the sample.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003ePerformance Metrics\u003c/h2\u003e \u003cp\u003eConfusion matrices for True Positive (TP), True Negatives (TN), False Positives (FP) and False Negatives (FN) were established for each tool, allowing calculation of Precision, Accuracy, Recall, False Positive Rate (FPR), F1 score, and Mathews Correlation Coefficient (MCC).\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:\\text{Recall}=\\frac{TP}{TP+FN}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$$\\:\\text{Precision}=\\frac{TP}{TP+FP}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equc\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e\n$$\\:\\text{Accuracy}=\\frac{TP+TN}{TP+TN+FP+FN}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equd\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equd\" name=\"EquationSource\"\u003e\n$$\\:\\text{FPR}=\\frac{FP}{FP+TN}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Eque\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Eque\" name=\"EquationSource\"\u003e\n$$\\:\\text{MCC}=\\frac{TP\\times\\:TN-FP\\times\\:FN}{\\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$$\u003c/div\u003e\u003c/div\u003e\u003cdiv id=\"Equf\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equf\" name=\"EquationSource\"\u003e\n$$\\:F1=2\\times\\:\\frac{\\text{Precision}\\times\\:\\text{Recall}}{\\text{Precision}+\\text{Recall}}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eVariant Annotation\u003c/h2\u003e \u003cp\u003eIndividual VCFs generated from exome sequencing and genome sequencing were annotated for population allele frequencies and splice prediction scores using the VEP_annotate tool (v113)[\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e]. SpliceAI scores were obtained as precomputed datasets from Biowulf[\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e]. Allele frequencies for each variant were sourced from gnomAD genomes v4.0 (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://gnomad.broadinstitute.org/\u003c/span\u003e\u003cspan address=\"https://gnomad.broadinstitute.org/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eSubsequent filtering of annotated VCFs was performed using VEP_filter (v113) to retain only variants with an allele frequency less than 0.01 or 0.0005 (a high and low threshold relevant to rare disease to assess how different cutoffs affect final outlier number) and with no homozygotes reported in gnomAD. Filtered VCFs were then converted to BED format using Bedops (version 2.4.41 megarow) with the convert2bed function.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eOutlier Evaluation\u003c/h2\u003e \u003cp\u003eFor splicing, an ensemble approach with variant filtration was used. Chromosomal coordinates for each outlier splice event were extracted into a master bed file (see Supplemental Methods). The master bed file was sorted using bedops (v2.4.41) and intersected with the filtered annotated variant bed file using the --range \u0026minus;\u0026thinsp;50:50 -e set operation to identify variants overlapping splice outliers. Variants that intersected the outlier bed file were then filtered for spliceAI score \u0026ge; 0.1, segregated in seqr and manually evaluated with the integrative genomics viewer (IGV)[\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e] (see Supplemental Methods). For OUTRIDER and ANEVA-h, the number of variants identified per individual was sufficiently low to allow for manual inspection of all events without additional filtration.\u003c/p\u003e \u003cp\u003eCandidate gene variants identified for undiagnosed patients were submitted to GeneMatcher[\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e].\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cp\u003e\u003cem\u003eTruth set selection\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eWe used RNA-seq data from 97 previously diagnosed patient samples to benchmark RNA-seq analysis tools. To ensure RNA-seq samples were correct, each BAM file was confirmed to harbor the known pathogenic variant in IGV. RNA-seq samples with independently verified variants affecting RNA were manually validated to contain evidence of aberrant splicing, expression or allelic imbalance in IGV. This set was defined as our true positive set (76 samples from 68 probands, Fig 1A,C, see methods for detailed definition). RNA-seq samples with confirmed pathogenic missense variants, inframe indels and premature truncating variants without nonsense mediated decay (NMD) and not expected to cause RNA changes were manually confirmed to contain no visible transcriptome changes in the pathogenic gene. This set was defined as our true negative set (21 samples, Fig 1A, see methods for detailed definition). \u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eRecall rate of splice tools identifies ensemble strategy as most sensitive\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo determine the performance of splice tools, we compared five tools using both case–control and one-vs-all designs to recall true positive splice events by tissue source (68 true positive variants from muscle, fibroblast and cycloheximide treated fibroblasts, Fig. 1A-B). In muscle derived true positive samples, case–control tools were most sensitive: rMATS (56%) \u0026gt; MINTIE (37.5%) \u0026gt; FRASER (33%) \u0026gt; FRASER2.0 (12.5%) \u0026gt; LeafcutterMD (0%) (16 true positive variants, Table S3, Fig. 2A). FRASER2.0 was excluded from further analysis as it underperformed FRASER. In fibroblast derived true positives, rMATS again performed best (86%) \u0026gt; MINTIE (28%) \u0026gt; FRASER (14%) \u0026gt; LeafcutterMD (0%) (7 true positive variants, Table S4, Fig 2A). Cycloheximide (CHX) treatment of fibroblasts can halt NMD by blocking protein synthesis, increasing detection of degraded aberrant events. In CHX-treated fibroblast true positives, only outlier-based methods were tested due to lack of control data, FRASER detected 33% and LeafcutterMD none (45 true positive variants, Table S5, Fig 2A).\u003c/p\u003e\n\u003cp\u003eTo evaluate if underperforming tools contributed unique detections, we examined recall by splice event category. FRASER detected two intron retention events, contributing one variant missed by other tools, supporting its inclusion (Fig 2B). MINTIE uniquely captured multiple splice events of a \u003cem\u003eBET1\u003c/em\u003e variant (NM_005868.6: c.202G\u0026gt;C; p.Asp68His) in fibroblast and muscle tissue (Fig 2B, Tables S3-4). Although case-control methods have the highest recall rates, an ensemble approach aggregating outliers across multiple tools is most sensitive. In muscle, when tool output is combined, 68% (11/16) true positives are recalled, indicating greater overall yield with ensemble strategies (Fig 2B).\u003c/p\u003e\n\u003cp\u003eTo understand the limitations of tool use, we examined the true positives undetected by any tool (Tables S3-5). These reflected common limitations: low mapping quality (\u003cem\u003ePOMGNT1\u0026nbsp;\u003c/em\u003eNM_001290129.2 c.636C\u0026gt;T; p.205 =), few spanning splice junction reads due to degradation (\u003cem\u003eRYR1\u003c/em\u003e NM_000540.3: c.3178+587A\u0026gt;G;\u003cem\u003e\u0026nbsp;NEB\u0026nbsp;\u003c/em\u003eNM_001164508.2: c.21417+3A\u0026gt;G), complex loci with age specific alternative exon usage or highly repetitive regions (\u003cem\u003eNEB\u0026nbsp;\u003c/em\u003eNM_001164508.2: c.21417+3A\u0026gt;G, \u003cem\u003eIGHMBP2\u003c/em\u003e NM_002180.3: c.1235+450 G\u0026gt;A), or inaccurate splice junction splitting (\u003cem\u003eRYR1\u003c/em\u003e NM_000540.3: c.3178+587A\u0026gt;G\u003cem\u003e\u0026nbsp;, COL6A1\u003c/em\u003e NM_001848.3: c.190+176C\u0026gt;T, \u003cem\u003eIGHMBP2\u0026nbsp;\u003c/em\u003eNM_002180.3: c.1235+450 G\u0026gt;A). These cases demonstrate the continuing importance of manual investigation guided by clinical and genetic hypothesis.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eClinical application of splice detection tools requires a manageable number of outliers for review. We assessed output interpretability by comparing outlier counts reported on average per sample. In muscle, MINTIE and rMATS reported thousands of outliers, particularly when strand-specific data were compared to unstranded controls (Fig. 2C). In fibroblast and cycloheximide-treated fibroblast samples, outlier counts clustered by sequencing center, reflecting batch and library preparation effects (Figs. S4A–D).\u003c/p\u003e\n\u003cp\u003eTo reduce outlier burden, we filtered for events overlapping rare DNA variants (Fig. 2D). Among ten true positive samples with genomic data, filtering for AF ≤ 0.01, a lenient cutoff for rare disease, yielded a mean of 1,221 variants per sample; applying a stricter threshold (AF ≤ 0.0005, segregation when available, and SpliceAI ≥ 0.1) reduced this to 1–3 variants (Figs. 2D–E, Table S6, Supplemental Methods). In nine of ten cases, the causal variant remained post-filtering. The missed case (\u003cem\u003ePOMT2\u003c/em\u003e NM_013382.7) the first variant (c.333+1G\u0026gt;T) associated splice event was not included in the master bed file (see Supplemental Methods), while a second intronic variant (c.1725+78G\u0026gt;A) was not captured by exome sequencing, highlighting the added value of including RNA-based variant calling (Table S6). Overall, rare variant filtration is successful at overcoming the challenge presented by the high number of outlier calls.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eRecall rate of expression tools is low\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo determine the utility of expression analysis, we benchmarked OUTRIDER (15 true positive genes, Fig 1A-B). Across muscle true positive cases the recall rate was 7%, and the addition of 100 healthy control samples did not improve the recall (14 true positive genes,Table 1, Table S7). Fibroblast true positives were not detected at all, and no true positives were assessed from the CHX samples (1 true positive gene, Table 1, Table S7). Overall, the combined recall rate was 6.6% (1/15), and inclusion of expression analysis did not identify diagnoses missed by other analysis types (Fig 1C). The average number of outliers reported per sample were low (Table 1). Notably, the inclusion of control samples increased the average number of outliers while not providing additional diagnostic benefit (Table 1).\u003c/p\u003e\n\u003cp\u003eTo understand the limitations of OUTRIDER, we evaluated the undetected true positives. \u003cem\u003eTTN\u003c/em\u003e underexpression recurs four times in our group of true positives and may not be considered statistically significant when using the one vs all approach. No patterns emerged for why any of the other true positives were missed (Table S7).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eRecall rate of allelic imbalance highlights ANEVA-H\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo benchmark allelic imbalance detection, we tested the MAE module from DROP and ANEVA-H (8 true positive genes, Fig. 1A-B). Both detected 5 of 8 true positives (62.5% recall; Fig. 3A, Table S8). The tools overlapped on four variants, each identifying one unique event, for a combined recall of 75% (6/8). Allelic imbalance analysis yielded 4 diagnoses missed by splice or expression outlier methods (68 total diagnoses, Fig. 1C), highlighting its complementary diagnostic value. ANEVA-H reported markedly fewer outliers per sample than MAE (mean 44 vs 167; Fig. 3B), supporting its clinical adaptability.\u003c/p\u003e\n\u003cp\u003eUndetected cases involved \u003cem\u003eMSTO1\u003c/em\u003eand \u003cem\u003eJAG2\u003c/em\u003e. \u003cem\u003eMSTO1 (\u003c/em\u003eNM_018116.4:c.1433A\u0026gt;G; p.Tyr478Cys \u0026amp; deletion) mapping is confounded by a pseudogene, and \u003cem\u003eJAG2\u003c/em\u003e \u003cem\u003e(\u003c/em\u003eNM_002226.5: c.2515G\u0026gt;A; p.Gly839Arg \u0026amp; large deletion) was listed as homozygous in the VCF, masking monoallelic expression (Table S8). Two true positives were only detected by one tool each. ANEVA-H specifically missed \u003cem\u003eRYR1\u003c/em\u003e (NM_000540.3: c.14126C\u0026gt;T; p.Thr4709Met \u0026amp; c.3178+587A\u0026gt;G) due to absent haplotype-specific SNPs in \u003cem\u003eRYR1\u003c/em\u003e (required for gene level ASE generation), and MAE specifically missed \u003cem\u003eTTN\u003c/em\u003e (NM_001267550.2: c.35794G\u0026gt;T; p.Glu11932Ter \u0026amp; c.22973C\u0026gt;A; p.Ser7658Ter) because of its size and repetitive regions (Table S8). Gene-level estimates from ANEVA-H may be more sensitive to subtle effects than variant-level MAE output.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eTrue negative analysis reveals high false positive rates\u0026nbsp;\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eUnnecessary evaluation of RNA outliers in patients whose final diagnoses do not involve RNA-mediated mechanisms can divert limited clinical resources. To estimate how often such outliers occur, twenty-one true negative (defined in methods) RNA-seq samples were analyzed across all computational tools for presence or absence of outlier calls, as presence of any number of outliers will require clinical evaluation without benefit (Fig 1A). Among splicing tools, rMATS, MINTIE, LeafCutterMD, and FRASER reported outliers for all true negative samples (100%), while FRASER2.0 had a lower rate (33%) (Fig. 4A, Tables S3-5, Fig. S5). For expression tools, OUTRIDER with GTEx controls showed higher false positive rates than OUTRIDER alone (83% vs. 56%), confirming that OUTRIDER alone is both more sensitive and specific (Table S7, Fig 4A). For allelic imbalance tools, MAE produced more false positives (100%) than ANEVA-H (80%), further supporting ANEVA-H as the better-performing method given its lower outlier burden (Table S8, Fig. 4A).\u003c/p\u003e\n\u003cp\u003e\u003cem\u003ePoor predictive performance across tools\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eTo evaluate tool accuracy and precision, true positive and true negative detections were aggregated across tissue types for each tool (Table S9). Among splicing tools, rMATS achieved the highest precision, accuracy, recall, and F1 score (Fig. 4B, Table S9). For expression analysis, OUTRIDER alone outperformed the OUTRIDER + GTEx combination across all metrics (Fig. 4C, Table S9). In allelic imbalance analysis, ANEVA-H exceeded MAE on every metric (Fig. 4D, Table S9). Notably, Mathews Correlation Coefficient (MCC) values were negative across all tools, indicating limited overall predictive power (Fig 4B-D, Table S9). Across all analyses, only 28 of 68 (41%) known diagnoses were correctly detected (Fig. 1C), consistent with the MCC scores.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eApplication to undiagnosed samples\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eWe applied the top performing approaches integrating ensemble splice analysis with variant filtration, OUTRIDER, and ANEVA-H to 74 undiagnosed samples (Fig. 5A). In 28 muscle-derived samples, splice analysis identified 129 candidate variants (mean 4.6 per sample), of which 39 showed clear aberrant splicing on manual review. Nine variants were prioritized for further evaluation after segregation and analysis for a possible second hit if not de novo (Figs. 5B-C, Table S10). None were located in OMIM-listed disease genes. OUTRIDER analysis identified one candidate in a CHX treated sample, \u003cstrong\u003e\u003cem\u003eTRIP4\u003c/em\u003e\u003c/strong\u003e, with reduced expression, in a patient with phenotypically compatible disease (Fig. S6A–B); however, exome, genome, and optical mapping revealed no causal variants, and follow-up RNA-seq confirmed normal expression and splicing in untreated fibroblasts (Fig. S6C–D). ANEVA-H did not identify additional candidates. Together, these analyses demonstrate that splicing-based approaches contributed novel candidate findings, whereas expression and allelic imbalance analyses did not improve diagnostic yield.\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eOver a 14-year period, RNA-seq was integrated into our research diagnostic pipeline for patients with suspected Mendelian neuromuscular disorders, improving diagnostic yield but leaving a substantial fraction unresolved. To assess how transcriptome-centric methods could further enhance yield, we benchmarked commonly used RNA analysis tools.\u003c/p\u003e \u003cp\u003eSplicing tools consistently demonstrated the highest diagnostic sensitivity, detecting 24/68 diagnoses, in line with prior studies [\u003cspan additionalcitationids=\"CR6\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. While outlier-based methods FRASER and FRASER2 are the most widely implemented, case\u0026ndash;control approaches such as rMATS-turbo achieved greater sensitivity in our dataset. A prior report analyzing whole blood RNA-seq reported similar findings, while another report using muscle-derived RNA-seq found greater sensitivity with FRASER [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e, \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. Differences in FDR thresholds applied to rMATS and FRASER may explain the discrepancy [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e]. Although combining multiple tools improved recall in our dataset, it also increased outlier burden. Reported strategies to reduce outliers include restricting to clinically relevant genes or panels [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e, \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. While this enhances interpretability, it does not consistently result in new diagnoses [\u003cspan additionalcitationids=\"CR16\" citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. We chose not to restrict results this way given that our patients had undergone prior extensive phenotype based genetic analyses without success.\u003c/p\u003e \u003cp\u003eTreatment of cells with cycloheximide has been proposed to enhance detection of splice events subject to NMD, while the diagnostic yield reported has only been modest [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e]. With only one sample recurring in our treated and untreated true positives, we could not perform comprehensive analysis to ascertain the diagnostic benefit of using CHX.\u003c/p\u003e \u003cp\u003eExpression outlier analysis assumes that loss-of-function (LoF) variants trigger NMD and reduce transcript abundance. Using OUTRIDER, with and without GTEx controls, we observed limited sensitivity and no novel diagnostic findings, consistent with previous reports [\u003cspan additionalcitationids=\"CR6\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. Some outliers corresponded to known heterozygous pathogenic variants, while many others lacked supporting rare variants, likely reflecting pathway-level variation or common polymorphisms [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e]. Tissue specificity and genome build inconsistencies may further confound interpretation [\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e, \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e]. Nevertheless, higher yields reported in other cohorts\u0026mdash;such as NIH Undiagnosed Diseases Program fibroblast datasets[\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e] \u0026mdash;suggest that expression analysis may perform better in phenotypically diverse or batch sequenced cohorts [\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e, \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eAllelic imbalance analysis remains underutilized in RNA-seq diagnostics, with only a few studies reporting its clinical application and modest diagnostic yield [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e]. While variant level MAE and gene level ANEVA-H have similar sensitivity in our dataset, ANEVA-H produces fewer outliers and has a lower false positive rate. Genes reported as outliers by ANEVA-H frequently harbored structural variants or LoF variants, consistent with prior studies [\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e], yet lacked compatible inheritance or phenotypic correlation to make them candidate diagnoses. Imbalance analysis recovered four diagnoses in our true positive benchmarking missed by other tools, supporting its inclusion in RNA analysis pipelines. However, methods such as ANEVA-H remain constrained by haplotype phasing requirements, supporting ensemble approaches for now.\u003c/p\u003e \u003cp\u003eThis study has several limitations. An optimal study would be prospective, using pre-defined criteria for negative and positive cases. However, our study required a minority subset (splicing and transcription modifying variants) of a set of rare diseases. We felt that a study of retrospective data, with pre-analytically defined case categorization, had a strong potential to generate useful information. That said, our findings would benefit from prospective implementation and replication. Additionally, variants with clinically interpretable mechanisms were selected for assessment of tools, reflective of downstream evaluation of outliers in the clinic. This favors clearly defined aberrant events and likely underestimates tool sensitivity, making discernment of a single tool\u0026rsquo;s superiority challenging. True negatives were evaluated by presence or absence of outliers to capture time spent assessing outliers for patients without underlying pathogenic aberrant RNA events. Evaluating true negatives by binary outlier classification may overestimate false positive rates, as some outliers may represent genuine biological variation rather than artifacts, and additionally, underestimates predictive ability of tools by lowering MCC scores.\u003c/p\u003e \u003cp\u003eChallenges in identifying clear diagnostic candidates in our 74 unresolved cases likely reflect both biological and technical factors. Firstly, unsolved cases may harbor non-RNA-affecting pathogenic mechanisms which would be missed by RNA analytical tools [\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e, \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e, \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e]. Secondly, muscle and fibroblast tissues may not be the relevant tissue for capturing pathogenic events in patients with neurogenic or neuromuscular junction phenotypes. Moreover, reliance on exome data in some cases excludes deep intronic variants, and addition of RNA variant calling or genome sequencing could provide more comprehensive detection.\u003c/p\u003e \u003cp\u003eWhether to apply RNA-centric versus RNA-complementary approaches remain a central question for diagnostic transcriptomics. Manual investigation of genomic variants in RNA-seq files is inherently complementary. In our benchmarking, RNA-based tools only recovered 28 of 68 known diagnoses with RNA events, and application of tools to undiagnosed cases required rare variant filtration, underscoring the interdependence of transcriptomic and genomic data with current technology. Similar trends have been reported across studies, where diagnostic gains stem largely from RNA evidence supporting previously identified VUSs or genes prioritized due to one pathogenic hit and phenotypic fit [\u003cspan additionalcitationids=\"CR6\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e, \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. RNA-centric pipelines relying solely on RNA variant calls instead may miss pathogenic splice events when the variant is absent in the resulting transcript, while hybrid genome-plus-transcriptome strategies offer the most comprehensive diagnostic coverage [\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e, \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e]. Given the performance of the current tools, RNA approaches remain complementary to genomic analyses, with manual investigation still a necessary analytical tool.\u003c/p\u003e"},{"header":"Conclusion","content":"\u003cp\u003eWith currently available RNA-seq tools, we recommend using ensemble approaches in diagnostic pipelines rather than any single tool. Additionally, we recommend manual inspection of all VUS\u0026rsquo;s predicted to alter RNA even if no aberrant events are identified using tools. As RNA-seq continues to mature as a diagnostic tool, integrated genome-transcriptome pipelines, improved reference datasets, disease specific transcriptomic signatures and tissue-specific expression ranges will be essential for reducing false positives and enhancing interpretability [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e, \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e]. Increased use of long-read RNA-seq will improve splice detection in difficult to map regions and genes with complex isoforms. Ultimately, such unified analytical frameworks will expedite variant prioritization and accelerate the path to molecular diagnosis.\u003c/p\u003e"},{"header":"Abbreviations","content":"\u003cdiv class=\"DefinitionList\"\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eNMD\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eNonsense Mediated Decay\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eCHX\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eCycloheximide\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eASE\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eAllele Specific Expression\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eIGV\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eIntegrative Genomics Viewer\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eRNA-Seq\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eRNA-Sequencing\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eVUS\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eVariant of Uncertain Significance\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eLoF\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eLoss of Function\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eNIH\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eNational Institutes of Health\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eSNP\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eSingle Nucleotide Polymorphism\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eBAM\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eBinary Alignment File\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eAF\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eAllele Frequency\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eVCF\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eVariant Call File\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eMCC\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eMathews Correlation Coefficient\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eTPM\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eTranscripts Per Million\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv class=\"DefinitionListEntry\"\u003e \u003cdiv class=\"Term\"\u003eFPR\u003c/div\u003e \u003cdiv class=\"Description\"\u003e \u003cp\u003eFalse Positive Rate\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cem\u003eEthics approval and consent to participate\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eWritten informed consent and age- appropriate assent for study procedures were obtained by a qualified investigator [protocol 12-N-0095 approved by the National Institute of Neurological Disorders and Stroke, National Institutes of Health IRB].\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eConsent for publication\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eAll participating patients consented to publication of deidentified material\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eAvailability of data and materials\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe datasets generated and/or analyzed during the current study are available in the dbGAP repository phs001272.v2.p1.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eCompeting interests\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eFunding\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThe Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from dbGaP accession number phs000424.v8.p2 on 09/25/2023.\u003c/p\u003e\n\u003cp\u003eThis research was supported by the Intramural Research Program of the National Institutes of Health (NIH) [grant 1ZIANS003129]. The contributions of the NIH author(s) are considered Works of the United States Government. The findings and conclusions presented in this paper are those of the author(s) and do not necessarily reflect the views of the NIH or the U.S. Department of Health and Human Services.\u003c/p\u003e\n\u003cp\u003eSequencing and analysis were provided by the Broad Institute of MIT and Harvard Center for Mendelian Genomics (Broad CMG) and was funded by the National Human Genome Research Institute, the National Eye Institute, the National Heart, Lung and Blood Institute and grant UM1HG008900 to Daniel MacArthur and Heidi Rehm. This study makes use of data shared through the seqr platform with funding provided by National Institutes of Health grants R01HG009141 and UM1HG008900.\u003c/p\u003e\n\u003cp\u003ePM was supported by the National Institutes of Health under award number R01GM140287\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eAuthors\u0026apos; contributions\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eSS conceptualized the study design, performed all data analysis, and prepared the manuscript. KG ran the ANEVA-H data analysis and contributed manuscript preparation. SG contributed to the study design and manuscript preparation. SD performed initial genetic analysis for each patient used in the study and contributed to clinical interpretation of new candidate diagnoses. VB collected RNA from patients with COL6-RDs. JM ran FRASER and OUTRIDER on the cycloheximide treated RNA-seq samples. YH collected RNA from patient muscle and fibroblast for sequencing. PU contributed to study design and bioinformatic support. VG and BW contributed patient diagnoses used as true positives in the study. RO contribute to clinical evaluation of patients and interpretation of new candidate diagnoses. ARF contribute to clinical evaluation of patients. PM supervised ANEVA-H data analysis. DA supervised SS and conceptualized study design and manuscript preparation. CB supervised SS and conceptualized study design and manuscript preparation. All authors read and approved the final manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cem\u003eAcknowledgements\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eThis work utilized the computational resources of the NIH HPC Biowulf cluster\u003cem\u003e\u0026nbsp;(\u003c/em\u003e\u003cem\u003ehttps://hpc.nih.gov\u003c/em\u003e\u003cem\u003e).\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003eWith great thanks to referring physicians, patients and their families. We would like to thank Christopher Mendoza, Gilberto (\u0026ldquo;Mike\u0026rdquo;) Averion and Kia Brooks for their help and support in coordinating the patient evaluations.\u0026nbsp;\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eZambon AA, Muntoni F. Congenital muscular dystrophies: What is new? Neuromuscular Disorders. Elsevier; 2021;31:931\u0026ndash;42. https://doi.org/10.1016/J.NMD.2021.07.009\u003c/li\u003e\n\u003cli\u003eZambon AA, Pini V, Bosco L, Falzone YM, Munot P, Muntoni F, et al. Early onset hereditary neuronopathies: an update on non-5q motor neuron diseases. Brain [Internet]. Oxford Academic; 2023 [cited 2024 Jan 9];146:806\u0026ndash;22. https://doi.org/10.1093/BRAIN/AWAC452\u003c/li\u003e\n\u003cli\u003eKang PB, Morrison L, Iannaccone ST, Graham RJ, B\u0026ouml;nnemann CG, Rutkowski A, et al. Evidence-based guideline summary: evaluation, diagnosis, and management of congenital muscular dystrophy: Report of the Guideline Development Subcommittee of the American Academy of Neurology and the Practice Issues Review Panel of the American Association of Neuromuscular \u0026amp; Electrodiagnostic Medicine. Neurology [Internet]. Neurology; 2015 [cited 2023 Dec 13];84:1369\u0026ndash;78. https://doi.org/10.1212/WNL.0000000000001416\u003c/li\u003e\n\u003cli\u003eMaggi L, Scoto M, Cirak S, Robb SA, Klein A, Lillis S, et al. Congenital myopathies--clinical features and frequency of individual subtypes diagnosed over a 5-year period in the United Kingdom. Neuromuscul Disord [Internet]. Neuromuscul Disord; 2013 [cited 2024 Jan 8];23:195\u0026ndash;205. https://doi.org/10.1016/J.NMD.2013.01.004\u003c/li\u003e\n\u003cli\u003eSegarra-Casas A, Dom\u0026iacute;nguez-Gonz\u0026aacute;lez C, Natera-de Benito D, Kapetanovic S, Hern\u0026aacute;ndez-La\u0026iacute;n A, Est\u0026eacute;vez-Arias B, et al. Translating Muscle RNAseq Into the Clinic for the Diagnosis of Muscle Diseases. Ann Clin Transl Neurol [Internet]. John Wiley and Sons Inc; 2025 [cited 2025 Jun 25];0:1\u0026ndash;15. https://doi.org/10.1002/ACN3.70078;REQUESTEDJOURNAL:JOURNAL:23289503;JOURNAL:JOURNAL:23289503;WGROUP:STRING:PUBLICATION\u003c/li\u003e\n\u003cli\u003eLuo X, Xiao B, Liang L, Zhang K, Xu T, Liu H, et al. Blood RNA-seq in rare disease diagnostics: a comparative study of cases with and without candidate variants. J Transl Med [Internet]. BioMed Central Ltd; 2025 [cited 2025 Jun 25];23. https://doi.org/10.1186/S12967-025-06609-W,\u003c/li\u003e\n\u003cli\u003eJaramillo Oquendo C, Wai HA, Rich WI, Bunyan DJ, Thomas NS, Hunt D, et al. Identification of diagnostic candidates in Mendelian disorders using an RNA sequencing-centric approach. Genome Med [Internet]. BioMed Central Ltd; 2024 [cited 2025 Jun 25];16:1\u0026ndash;19. https://doi.org/10.1186/S13073-024-01381-W/FIGURES/6\u003c/li\u003e\n\u003cli\u003eFr\u0026eacute;sard L, Smail C, Ferraro NM, Teran NA, Li X, Smith KS, et al. Identification of rare-disease genes using blood transcriptome sequencing and large control cohorts. Nat Med. Nature Research; 2019;25:911\u0026ndash;9. https://doi.org/10.1038/s41591-019-0457-8\u003c/li\u003e\n\u003cli\u003eKremer LS, Bader DM, Mertes C, Kopajtich R, Pichler G, Iuso A, et al. Genetic diagnosis of Mendelian disorders via RNA sequencing. Nat Commun. Nature Publishing Group; 2017;8. https://doi.org/10.1038/ncomms15824\u003c/li\u003e\n\u003cli\u003eCummings BB, Marshall JL, Tukiainen T, Lek M, Donkervoort S, Foley AR, et al. Improving genetic diagnosis in Mendelian disease with transcriptome sequencing. Sci Transl Med [Internet]. American Association for the Advancement of Science; 2017 [cited 2024 Jan 16];9. https://doi.org/10.1126/scitranslmed.aal5209\u003c/li\u003e\n\u003cli\u003eMurdock DR, Dai H, Burrage LC, Rosenfeld JA, Ketkar S, M\u0026uuml;ller MF, et al. Transcriptome-directed analysis for Mendelian disease diagnosis overcomes limitations of conventional genomic testing. J Clin Invest [Internet]. American Society for Clinical Investigation; 2021 [cited 2024 Jan 16];131. https://doi.org/10.1172/JCI141500\u003c/li\u003e\n\u003cli\u003eGonorazky HD, Naumenko S, Ramani AK, Nelakuditi V, Mashouri P, Wang P, et al. Expanding the Boundaries of RNA Sequencing as a Diagnostic Tool for Rare Mendelian Disease. Am J Hum Genet [Internet]. Am J Hum Genet; 2019 [cited 2024 Jan 16];104:466\u0026ndash;83. https://doi.org/10.1016/J.AJHG.2019.01.012\u003c/li\u003e\n\u003cli\u003eLee H, Huang AY, Wang L kai, Yoon AJ, Renteria G, Eskin A, et al. Diagnostic utility of transcriptome sequencing for rare Mendelian diseases. Genetics in Medicine. Elsevier; 2020;22:490\u0026ndash;9. https://doi.org/10.1038/S41436-019-0672-1\u003c/li\u003e\n\u003cli\u003eY\u0026eacute;pez VA, Gusic M, Kopajtich R, Mertes C, Smith NH, Alston CL, et al. Clinical implementation of RNA sequencing for Mendelian disease diagnostics. Genome Medicine 2022 14:1 [Internet]. BioMed Central; 2022 [cited 2022 May 4];14:1\u0026ndash;26. https://doi.org/10.1186/S13073-022-01019-9\u003c/li\u003e\n\u003cli\u003eDe Cock L, D\u0026rsquo;haenens E, Vantomme L, Backers L, Beyens A, Claes KB, et al. Cracking rare disorders: a new minimally invasive RNA-seq protocol. NPJ Genom Med [Internet]. NPJ Genom Med; 2025 [cited 2025 Jun 25];10:45. https://doi.org/10.1038/S41525-025-00502-7,\u003c/li\u003e\n\u003cli\u003eRiquin K, Isidor B, Mercier S, Nizon M, Colin E, Bonneau D, et al. Integrating RNA-Seq into genome sequencing workflow enhances the analysis of structural variants causing neurodevelopmental disorders. J Med Genet [Internet]. BMJ Publishing Group; 2023 [cited 2025 Jun 25];61:47\u0026ndash;56. https://doi.org/10.1136/JMG-2023-109263,\u003c/li\u003e\n\u003cli\u003eDekker J, Schot R, Bongaerts M, de Valk WG, van Veghel-Plandsoen MM, Monfils K, et al. Web-accessible application for identifying pathogenic transcripts with RNA-seq: Increased sensitivity in diagnosis of neurodevelopmental disorders. Am J Hum Genet [Internet]. Cell Press; 2023 [cited 2025 Jun 25];110:251\u0026ndash;72. https://doi.org/10.1016/j.ajhg.2022.12.015\u003c/li\u003e\n\u003cli\u003eBournazos AM, Riley LG, Bommireddipalli S, Ades L, Akesson LS, Al-Shinnag M, et al. Standardized practices for RNA diagnostics using clinically accessible specimens reclassifies 75% of putative splicing variants. Genetics in Medicine. Elsevier; 2022;24:130\u0026ndash;45. https://doi.org/10.1016/J.GIM.2021.09.001\u003c/li\u003e\n\u003cli\u003eWang Y, Xie Z, Kutschera E, Adams JI, Kadash-Edmondson KE, Xing Y. rMATS-turbo: an efficient and flexible computational tool for alternative splicing analysis of large-scale RNA-seq data. Nat Protoc [Internet]. Springer Nature; 2024 [cited 2025 Jun 25];19:1083\u0026ndash;104. https://doi.org/10.1038/S41596-023-00944-2,\u003c/li\u003e\n\u003cli\u003eScheller IF, Lutz K, Mertes C, Y\u0026eacute;pez VA, Gagneur J. Improved detection of aberrant splicing with FRASER 2.0 and the intron Jaccard index. Am J Hum Genet [Internet]. Cell Press; 2023 [cited 2025 Jun 25];110:2056\u0026ndash;67. https://doi.org/10.1016/j.ajhg.2023.10.014\u003c/li\u003e\n\u003cli\u003eBrechtmann F, Mertes C, Matusevičiūtė A, Y\u0026eacute;pez VA, Avsec Ž, Herzog M, et al. OUTRIDER: A Statistical Method for Detecting Aberrantly Expressed Genes in RNA Sequencing Data. Am J Hum Genet [Internet]. Am J Hum Genet; 2018 [cited 2024 Feb 5];103:907\u0026ndash;17. https://doi.org/10.1016/J.AJHG.2018.10.025\u003c/li\u003e\n\u003cli\u003eCmero M, Schmidt B, Majewski IJ, Ekert PG, Oshlack A, Davidson NM. MINTIE: identifying novel structural and splice variants in transcriptomes using RNA-seq data. Genome Biol [Internet]. BioMed Central Ltd; 2021 [cited 2024 Feb 5];22:1\u0026ndash;26. https://doi.org/10.1186/S13059-021-02507-8/TABLES/2\u003c/li\u003e\n\u003cli\u003eMertes C, Scheller IF, Y\u0026eacute;pez VA, \u0026Ccedil;elik MH, Liang Y, Kremer LS, et al. Detection of aberrant splicing events in RNA-seq data using FRASER. Nature Communications 2021 12:1 [Internet]. Nature Publishing Group; 2021 [cited 2024 Feb 5];12:1\u0026ndash;13. https://doi.org/10.1038/s41467-020-20573-7\u003c/li\u003e\n\u003cli\u003eY\u0026eacute;pez VA, Mertes C, M\u0026uuml;ller MF, Klaproth-Andrade D, Wachutka L, Fr\u0026eacute;sard L, et al. Detection of aberrant gene expression events in RNA sequencing data. Nature Protocols 2021 16:2 [Internet]. Nature Publishing Group; 2021 [cited 2024 Jan 24];16:1276\u0026ndash;96. https://doi.org/10.1038/s41596-020-00462-5\u003c/li\u003e\n\u003cli\u003eZhao S, Macakova K, Sinson JC, Dai H, Rosenfeld J, Zapata GE, et al. Clinical validation of RNA sequencing for Mendelian disorder diagnostics. Am J Hum Genet [Internet]. Cell Press; 2025 [cited 2025 Jun 30];112:779\u0026ndash;92. https://doi.org/10.1016/j.ajhg.2025.02.006\u003c/li\u003e\n\u003cli\u003eHe W Bin, Xiao WJ, Dai CL, Wang YR, Li XR, Gong F, et al. RNA splicing analysis contributes to reclassifying variants of uncertain significance and improves the diagnosis of monogenic disorders. J Med Genet [Internet]. BMJ Publishing Group; 2022 [cited 2025 Jun 29];59:1010\u0026ndash;6. https://doi.org/10.1136/JMEDGENET-2021-108013,\u003c/li\u003e\n\u003cli\u003eStark JC, Pipko N, Liang Y, Szuto A, Ting Tsoi C, Dickson MA, et al. Clinical applications of and molecular insights from RNA sequencing in a rare disease cohort. Genome Medicine 2025 17:1 [Internet]. BioMed Central; 2025 [cited 2025 Jul 14];17:1\u0026ndash;24. https://doi.org/10.1186/S13073-025-01494-W\u003c/li\u003e\n\u003cli\u003eLewis RG, O\u0026rsquo;Shea JM, Pizzo L, Wen T, Fulmer ML, Zhao J, et al. RNA sequencing provides functional insights and diagnostic resolution in previously unsolved rare disease cases. BMC Medical Genomics 2025 18:1 [Internet]. BioMed Central; 2025 [cited 2025 Nov 16];18:182-. https://doi.org/10.1186/S12920-025-02227-Z\u003c/li\u003e\n\u003cli\u003eDobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics [Internet]. Bioinformatics; 2013 [cited 2025 Jun 25];29:15\u0026ndash;21. https://doi.org/10.1093/BIOINFORMATICS/BTS635,\u003c/li\u003e\n\u003cli\u003eDeluca DS, Levin JZ, Sivachenko A, Fennell T, Nazaire MD, Williams C, et al. RNA-SeQC: RNA-seq metrics for quality control and process optimization. Bioinformatics [Internet]. Bioinformatics; 2012 [cited 2024 Feb 19];28:1530\u0026ndash;2. https://doi.org/10.1093/BIOINFORMATICS/BTS196\u003c/li\u003e\n\u003cli\u003eRobinson JT, Thorvaldsd\u0026oacute;ttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol [Internet]. Nat Biotechnol; 2011 [cited 2025 Jun 25];29:24\u0026ndash;6. https://doi.org/10.1038/NBT.1754,\u003c/li\u003e\n\u003cli\u003eThorvaldsdottir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013;14:178\u0026ndash;92. https://doi.org/10.1093/bib/bbs017\u003c/li\u003e\n\u003cli\u003eMontgomery SB, Bernstein JA, Wheeler MT. Toward transcriptomics as a primary tool for rare disease investigation. Cold Spring Harb Mol Case Stud [Internet]. Cold Spring Harbor Laboratory Press; 2022 [cited 2025 Jun 26];8:a006198. https://doi.org/10.1101/MCS.A006198\u003c/li\u003e\n\u003cli\u003eJenkinson G, Li YI, Basu S, Cousin MA, Oliver GR, Klee EW. LeafCutterMD: an algorithm for outlier splicing detection in rare diseases. Bioinformatics [Internet]. Bioinformatics; 2020 [cited 2022 Nov 3];36:4609\u0026ndash;15. https://doi.org/10.1093/BIOINFORMATICS/BTAA259\u003c/li\u003e\n\u003cli\u003eLi YI, Knowles DA, Humphrey J, Barbeira AN, Dickinson SP, Im HK, et al. Annotation-free quantification of RNA splicing using LeafCutter. Nature Genetics 2017 50:1 [Internet]. Nature Publishing Group; 2017 [cited 2022 Aug 21];50:151\u0026ndash;8. https://doi.org/10.1038/s41588-017-0004-9\u003c/li\u003e\n\u003cli\u003eGanapathy KR, Broly M, Silverstein S, Mendoza M, Song E, Kotis B, et al. Improved Identification of Large-effect Rare Genetic Variants using Haplotype Aggregated Allele-specific Expression Data. medRxiv [Internet]. Cold Spring Harbor Laboratory Press; 2025 [cited 2026 Jan 4];2025.12.16.25341855. https://doi.org/10.64898/2025.12.16.25341855\u003c/li\u003e\n\u003cli\u003eLoh PR, Danecek P, Palamara PF, Fuchsberger C, Reshef YA, Finucane HK, et al. Reference-based phasing using the Haplotype Reference Consortium panel. Nat Genet [Internet]. Nature Publishing Group; 2016 [cited 2025 Jun 25];48:1443\u0026ndash;8. https://doi.org/10.1038/NG.3679,\u003c/li\u003e\n\u003cli\u003eByrska-Bishop M, Evani US, Zhao X, Basile AO, Abel HJ, Regier AA, et al. High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios. Cell [Internet]. Elsevier B.V.; 2022 [cited 2025 Jun 25];185:3426-3440.e19. https://doi.org/10.1016/J.CELL.2022.08.004,\u003c/li\u003e\n\u003cli\u003eCastel SE, Mohammadi P, Chung WK, Shen Y, Lappalainen T. Rare variant phasing and haplotypic expression from RNA sequencing with phASER. Nat Commun [Internet]. Nature Publishing Group; 2016 [cited 2025 Jun 25];7:12817. https://doi.org/10.1038/NCOMMS12817\u003c/li\u003e\n\u003cli\u003eMohammadi P, Castel SE, Cummings BB, Einson J, Sousa C, Hoffman P, et al. Genetic regulatory variation in populations informs transcriptome analysis in rare disease. Science (1979) [Internet]. American Association for the Advancement of Science; 2019 [cited 2025 Jun 25];366:351\u0026ndash;6. https://doi.org/10.1126/SCIENCE.AAY0256/SUPPL_FILE/AAY0256_TABLES5.XLSX\u003c/li\u003e\n\u003cli\u003eMcLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The Ensembl Variant Effect Predictor. Genome Biol [Internet]. BioMed Central Ltd.; 2016 [cited 2025 Sep 14];17:1\u0026ndash;14. https://doi.org/10.1186/S13059-016-0974-4/TABLES/8\u003c/li\u003e\n\u003cli\u003eJaganathan K, Kyriazopoulou Panagiotopoulou S, McRae JF, Darbandi SF, Knowles D, Li YI, et al. Predicting Splicing from Primary Sequence with Deep Learning. Cell [Internet]. Cell Press; 2019 [cited 2025 Nov 4];176:535-548.e24. https://doi.org/10.1016/j.cell.2018.12.015\u003c/li\u003e\n\u003cli\u003ePais LS, Snow H, Weisburd B, Zhang S, Baxter SM, DiTroia S, et al. seqr: A web-based analysis and collaboration tool for rare disease genomics. Hum Mutat [Internet]. John Wiley and Sons Inc; 2022 [cited 2025 Jun 25];43:698\u0026ndash;707. https://doi.org/10.1002/HUMU.24366,\u003c/li\u003e\n\u003cli\u003eSobreira N, Schiettecatte F, Valle D, Hamosh A. GeneMatcher: A Matching Tool for Connecting Investigators with an Interest in the Same Gene. Hum Mutat [Internet]. John Wiley and Sons Inc.; 2015 [cited 2025 Sep 14];36:928\u0026ndash;30. https://doi.org/10.1002/HUMU.22844,\u003c/li\u003e\n\u003cli\u003eUngar RA, Goddard PC, Jensen TD, Degalez F, Smith KS, Jin CA, et al. Impact of genome build on RNA-seq interpretation and diagnostics. Am J Hum Genet [Internet]. Cell Press; 2024 [cited 2025 Jul 1];111:1282\u0026ndash;300. https://doi.org/10.1016/j.ajhg.2024.05.005\u003c/li\u003e\n\u003cli\u003eKaram R, Wengrod J, Gardner LB, Wilkinson MF. Regulation of nonsense-mediated mRNA decay: Implications for physiology and disease. Biochimica et Biophysica Acta (BBA) - Gene Regulatory Mechanisms. 2013;1829:624\u0026ndash;33. https://doi.org/10.1016/j.bbagrm.2013.03.002\u003c/li\u003e\n\u003cli\u003eCassini T, Silverstein S, Behan M, Tifft CJ, Malicdan MC, Adams DR, et al. Mitochondrial trifunctional protein deficiency caused by a deep intronic deletion leading to aberrant splicing. JIMD Rep [Internet]. John Wiley and Sons Inc; 2025 [cited 2025 Jul 1];66. https://doi.org/10.1002/JMD2.12459,\u003c/li\u003e\n\u003cli\u003eSilverstein S, Cassini T, Fu J, Pusey B, Macnamara E, Frost FG, et al. RNA sequencing driven diagnosis expands the phenotypic spectrum of NBAS deficiency. Mol Genet Metab [Internet]. Academic Press Inc.; 2025 [cited 2025 Jul 1];145. https://doi.org/10.1016/j.ymgme.2025.109105\u003c/li\u003e\n\u003cli\u003eFerraro NM, Strober BJ, Einson J, Abell NS, Aguet F, Barbeira AN, et al. Transcriptomic signatures across human tissues identify functional rare genetic variation. Science. 2020;369. https://doi.org/10.1126/science.aaz5900\u003c/li\u003e\n\u003cli\u003eSoemedi R, Cygan KJ, Rhine CL, Wang J, Bulacan C, Yang J, et al. Pathogenic variants that alter protein code often disrupt splicing. Nat Genet [Internet]. Nat Genet; 2017 [cited 2025 Nov 4];49:848\u0026ndash;55. https://doi.org/10.1038/NG.3837\u003c/li\u003e\n\u003cli\u003eChong R, Insigne KD, Yao D, Burghard CP, Wang J, Hsiao YHE, et al. A Multiplexed Assay for Exon Recognition Reveals that an Unappreciated Fraction of Rare Genetic Variants Cause Large-Effect Splicing Disruptions. Mol Cell [Internet]. Cell Press; 2019 [cited 2025 Nov 4];73:183-194.e8. https://doi.org/10.1016/j.molcel.2018.10.037\u003c/li\u003e\n\u003cli\u003eGuo H, Urban AE, Wong WH. Prioritizing disease-related rare variants by integrating gene expression data. PLoS Genet [Internet]. Public Library of Science; 2024 [cited 2025 Sep 1];20:e1011412. https://doi.org/10.1371/JOURNAL.PGEN.1011412\u003c/li\u003e\n\u003cli\u003eArriaga MT, Mendez R, Ungar RA, Bonner DE, Matalon DR, Lemire G, et al. Transcriptome-wide outlier approach identifies individuals with minor spliceopathies. medRxiv [Internet]. medRxiv; 2025 [cited 2025 Sep 15]; https://doi.org/10.1101/2025.01.02.24318941\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Table","content":"\u003cp\u003e\u003cem\u003eTable 1: Results of true positive expression analysis across Muscle, Fibroblast and Cycloheximide (CHX)- Treated Fibroblasts sets. Muscle samples were run with and without an additional 100 healthy samples from GTEx.\u003c/em\u003e\u003c/p\u003e\n\u003ctable border=\"0\" cellspacing=\"0\" cellpadding=\"0\" width=\"629\"\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 139px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eCohort\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 75px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eMuscle\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 99px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eMuscle + GTEx\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 93px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eFibroblast\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 223px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eCHX\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 139px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eSample size\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 75px;\"\u003e\n \u003cp\u003e71\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 99px;\"\u003e\n \u003cp\u003e171\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 93px;\"\u003e\n \u003cp\u003e37\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 223px;\"\u003e\n \u003cp\u003e79\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 139px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eRecall\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 75px;\"\u003e\n \u003cp\u003e1/14 (7%)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 99px;\"\u003e\n \u003cp\u003e1/14 (7%)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 93px;\"\u003e\n \u003cp\u003e0/1 (0%)\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 223px;\"\u003e\n \u003cp\u003eNA\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 139px;\"\u003e\n \u003cp\u003e\u003cstrong\u003eAvg Outliers/sample\u003c/strong\u003e\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 75px;\"\u003e\n \u003cp\u003e3.08\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 99px;\"\u003e\n \u003cp\u003e4.02\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 93px;\"\u003e\n \u003cp\u003e0.8\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 223px;\"\u003e\n \u003cp\u003e1.21\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n\u003c/table\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"RNA-seq, transcriptomics, rare disease, diagnostics, pediatric neuromuscular disease, computational tools, variant interpretation","lastPublishedDoi":"10.21203/rs.3.rs-8622051/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-8622051/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003ch2\u003eBackground\u003c/h2\u003e \u003cp\u003ePediatric neuromuscular diseases are genetically and clinically heterogeneous. A substantial proportion remain without a definitive genetic diagnosis despite available clinical molecular testing. RNA-sequencing (RNA-seq) can be used to complement genome or exome sequencing to elucidate or to identify the functional impact of variants of uncertain significance, but when manually analyzed is limited to candidate DNA variants or phenotype-driven gene lists. Open-source computational tools have been developed to systematically and unbiasedly analyze RNA-seq data for aberrant splicing, expression, or allelic imbalance. However, best use practice of these tools is yet to be established.\u003c/p\u003e\u003ch2\u003eMethods\u003c/h2\u003e \u003cp\u003eTo assess the performance of selected tools, we collected RNA-seq from 97 previously diagnosed samples to establish a truth set for benchmarking. Pathogenic variants were categorized as: true positives with confirmed aberrant RNA events and true negatives with no transcriptomic effect. We assessed performance of eight commonly used tools for splicing, expression and allelic imbalance analysis. We then applied the optimal strategy to 74 undiagnosed RNA-seq samples to identify new candidate diagnoses.\u003c/p\u003e\u003ch2\u003eResults\u003c/h2\u003e \u003cp\u003eAcross 68 diagnosed probands with aberrant RNA events, tools correctly identified 28 diagnoses. Splicing analysis tools provided most of the findings, but allelic imbalance tools uniquely identified 4, underscoring their value. Conversely, the false positive rate was highest for the splice tools and lowest for expression analysis. Application of tools led to identification of candidate variants for only 9 out of 74 undiagnosed patients.\u003c/p\u003e\u003ch2\u003eConclusions\u003c/h2\u003e \u003cp\u003eInclusion of RNA-seq tools can expedite variant prioritization, characterization and interpretation in the diagnostic pipeline but remain complementary to manual analysis of loci where candidate variants were identified by DNA sequencing.\u003c/p\u003e","manuscriptTitle":"Benchmarking RNA-seq Tools for Real-World Diagnostic Applications","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-01-29 00:53:21","doi":"10.21203/rs.3.rs-8622051/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"40cabce1-9ede-41aa-ad79-625c30ce7288","owner":[],"postedDate":"January 29th, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2026-02-09T05:23:56+00:00","versionOfRecord":[],"versionCreatedAt":"2026-01-29 00:53:21","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-8622051","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-8622051","identity":"rs-8622051","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}