Routine RNA-based analysis of potential splicing variants facilitates genomic diagnostics and reveals limitations of in silico prediction tools

Abstract

DNA variants affecting pre-mRNA splicing are an important cause of genetic disorders and remain challenging to interpret without experimental data. Although variant classification guidelines recommend experimental characterization of variant splicing effects, the added value of routine diagnostic investigation of patient mRNA splicing has not been systematically described. Here, we assessed the utility of pre-mRNA splicing analysis in a diagnostic setting for 202 suspected splice-altering variants from individuals referred for genetic testing. Pre-mRNA splicing was assessed in patient cells by RT-PCR, followed by agarose gel electrophoresis and Sanger sequencing and/or exon trapping assays. An effect on pre-mRNA splicing was demonstrated in 63% (n = 128/202) of the tested variants. Among the 177 variants initially classified as variants of uncertain significance (VUS), 54% (n = 96/177) were reclassified based on pre-mRNA splicing analysis, including 48% (n = 85/177) that were upgraded to likely pathogenic or pathogenic. We benchmarked the splice prediction algorithms SpliceAI, SQUIRLS, SPiP, and Pangolin, the tools integrated in Alamut on this clinically relevant and experimentally validated dataset, and the CAGI6 splicing VUS dataset and found variable performance dependent on variant type and location. No single tool classified all variants equally well. We describe several examples of hard-to-predict effects and unexpected results highlighting the limitations of prediction tools, including a not previously described variant type affecting U12-splice site subtype. In summary, we provide a framework for RNA-based analysis in a molecular diagnostic setting, demonstrate the added value of routine testing of RNA from individuals with suspected splice-altering variants, and highlight the limitations of in silico prediction tools.

Keywords: RNA splicing; diagnostics; genetic disorders; splice prediction algorithms.

Figure 1

Diagnostic analysis of mRNA splicing identifies various mRNA splicing effects and facilitates DNA variant reclassification (A) Schematic of the experimental setup. Patient tissue is collected (left), mRNA is isolated (second image from the left) and cDNA is generated using reverse transcriptase (second image from the right). A PCR fragment is subsequently amplified using primers that flank the variant (right). mRNA splicing effects are visualized after agarose gel electrophoresis and Sanger sequencing. (B–E) Doughnut plots showing characteristics of the tested DNA variants. (B) Disease categories to which the tested DNA variants were allocated. N.A., not available. (C) Variant classification of the tested DNA variants prior to mRNA splicing analysis. VUS, variant of uncertain significance. (D) Types of aberrant splicing events identified for samples in which aberrant splicing was found (n = 128/202). SS, splice site. (E) Effects on reading frame for samples in which aberrant splicing was found (n = 128/202). Unclear: in 1 case, the aberrant splicing product could barely be detected; it is likely degraded through NMD. Other: in both cases, the first coding exon of the gene (including the start codon) is skipped. (F) DNA variant reclassification of all VUS tested in this study (n = 177/202) after mRNA splicing analysis, combining the mRNA splicing results with other sources of relevant information, including variant frequency, segregation, and clinical phenotype.

Figure 2

SpliceAI and Pangolin most accurately predict effects on mRNA splicing (A) Bar chart depicting the positions of the DNA variants tested in this study relative to canonical splice donor (left) and splice acceptor (right) sites. Bars are color-coded, depicting the number of variants for which aberrant mRNA splicing was identified. (B) Violin plots showing algorithm scores for Splice AI, Pangolin, SPiP, and SQUIRLS stratified by experimental result (no effect on RNA splicing, in green, and effect on RNA splicing, in red), with each dot representing a single variant from the current dataset (lighter tone) or the CAGI6 dataset (darker tone). ∗∗∗∗p ≤ 0.0001, Kruskal-Wallis test with Dunn's multiple comparisons test. (C and D) (C) Area under the receiver operator characteristic curve (AUROC) and (D) area under the precision-recall curve (AUPRC) for Splice AI, Pangolin, SPiP, and SQUIRLS.

Figure 3

Homozygous in-frame deletion in CC2D2A as a cause for Joubert syndrome (A) Alamut splice prediction scores for wild type (top) and NM_001080522.2(CC2D2A):c.3976-3C>A, p.(?) (bottom). Note that splice acceptor prediction scores for exon 32 are only moderately decreased for the variant. (B) Agarose gel electrophoresis of RT-PCR products using primers for CCD2D2A exons 30 and 34. RT-PCR product size is decreased in the RNA isolated from fibroblasts from individual 1, indicated by the red arrow, both in the presence and absence of CHX treatment. Numerals on the left represent size marker (in nucleotides, first lane). NC, negative control. (C) Sanger sequencing of the amplified RT-PCR products from control (top) and individual 1 (bottom) tissues. Sequence analysis reveals skipping of exon 32 in ∼95% of amplified transcripts. (D) Schematic representation of the effect on mRNA splicing of the NM_001080522.2(CC2D2A):c.3976-3C>A variant in the tissue of individual 1. Genomic coordinates for chromosome 4 (GRCh37 (hg19)) are shown. (E and F) Brain MRI of individual 1 depicting a molar tooth sign (white arrows), typical of Joubert syndrome.

Figure 4

Pseudo-exon inclusion in ATM directs therapy choice (A) Alamut splice prediction scores for wild type (top) and NM_000051.3(ATM):c.3994-160T>C (bottom). Note the prediction of a non-canonical splice donor at position c.3994-161, which is lost in the variant, and the presence of both a strong non-canonical donor and an acceptor. (B) Agarose gel electrophoresis images after RT-PCR using primers for ATM exons 25 and 28. RT-PCR reveals 2 products of increased size in RNA isolated from blood from individual 2, indicated by 2 red arrows. Numerals on the left represent size marker (in nucleotides, first lane). (C) Alamut splice enhancer (exonic splice enhancer [ESE]) prediction for NM_000051.3(ATM):c.3994-160T>C (bottom) shows a gain of an SF2/ASF and an SF2/ASF(IgM-BRCA1) ESE. (D) Schematic representation of the effect on mRNA splicing in the blood of individual 2. Genomic coordinates for chromosome 11 (GRCh37 (hg19)) are shown.

Figure 5

Skipping of BLTP1 exon 75 through loss of an ESE (A) Volcano plot at gene level (z threshold = 4; p < 0.0025). BLTP1 is among the most strongly downregulated genes. (B) BLTP1 exon 75 ranking plot: individual 5 (red dot, indicated by the arrow, ∼55 counts) has a decreased expression of exon 75 compared to control samples (gray dots, average count ∼100). (C) Alamut ESE prediction for variants NM_015312.3(BLTP1):c.13212C>T, p.(Val4404 = ) and NM_015312.3(BLTP1):c.13222A>G, p.(Ile4408Val) combined (bottom) shows gain of an SF2/ASF ESE. (D) Integrative Genomics Viewer (IGV) Sashimi plots of RNA-seq data of untreated and CHX-treated samples of individual 4 and a control at the BLTP1 exon 74–76 region. Note the skipping of exon 75 in the individual’s sample, which is stabilized upon CHX treatment. (E) IGV plot of RNA-seq data (mapped reads) of untreated and CHX-treated samples of individual 4 and a control, showing allelic imbalance at positions NM_015312.3(BLTP1):c.13212C>T, p.(Val4404 = ) and NM_015312.3(BLTP1):c.13222A>G, p.(Ile4408Val). Calls are skewed toward the maternal transcript (C and A nucleotides, respectively) in the untreated cells, suggesting that exon 75 is skipped from the paternal transcript. (F) Exon trapping assay confirms exon skipping by variant NM_015312.3(BLTP1):c.13212C>T but not by NM_015312.3(BLTP1):c.13222A>G. Numerals on the left represent size marker (in nucleotides, first lane). (G) Schematic representation of the effect on mRNA splicing in the tissue of individual 4. Genomic coordinates for chromosome 4 (GRCh37 (hg19)) are shown.

Figure 6

U12 splice site donor/acceptor combination subtype switch induced by a splice donor variant in THOC2 (A) Alamut splice prediction scores for wild type (top left and right) and NM001081550.1(THOC2):c.4754+1A>G (bottom left). This variant is located at the intron 37 splice donor. The intron 37 splice acceptor is also depicted to illustrate its sequence context in relation to the U12 subtypes (bottom). (B) Agarose gel electrophoresis images after RT-PCR using primers for THOC2 exon 35 and the 3′ UTR. Numerals on the left represent size marker (in nucleotides, first lane). (C) Sanger sequencing of the amplified RT-PCR products from control (top) and the blood of individual 3 (bottom) reveals a 7-nt deletion at the exon 38 splice acceptor. (D) Schematic representation of the effect on mRNA splicing in the blood of individual 3. Genomic coordinates for chromosome X (GRCh37 (hg19)) are shown.