Bug Report: rMATS Reports MXE Events Not Supported by Splice Junctions in BAM Files or by Junctions Annotated in the GTF File

[santataRU]

Dear rMATS Developers,

I am using rMATS v4.3.0 on macOS (ARM64) and appreciate the effort you've put into this powerful tool. I understand that identifying MXE (mutually exclusive exon) events can be challenging. However, I have encountered a more serious issue: some reported MXE events do not appear to be supported by splice junctions in the BAM files or by Junctions Annotated in the GTF File.

Below is an example (hg38 reference genome):

ID	GeneID	geneSymbol	chr	strand	1stExonStart_0base	1stExonEnd	2ndExonStart_0base	2ndExonEnd	upstreamES	upstreamEE	downstreamES	downstreamEE	ID	IJC_SAMPLE_1	SJC_SAMPLE_1	IJC_SAMPLE_2	SJC_SAMPLE_2	IncFormLen	SkipFormLen	PValue	FDR	IncLevel1	IncLevel2	IncLevelDifference
2786	ENSG00000117448.14	AKR1A1	chr1	+	45552435	45552563	45561788	45561878	45550965	45551155	45567981	45568177	2786	107,93,129	146,152,188	161,136,158	101,92,86	277	239	0	0	0.387,0.346,0.372	0.579,0.561,0.613	-0.216

In this case, there is no junction supporting splicing from the 1stExon to the downstream exon (I run with --novelSS turned off), as illustrated in the attached screenshot. This raises concerns about many more potential false positives.

To help troubleshoot, I’ve shared a folder containing the rMATS MXE output, GTF files, and the BAM file corresponding to this region: (https://drive.google.com/drive/folders/1BtNkHl2UtYMAhnDxbpb6hvJaDWTWgbZ4?usp=sharing).

Please let me know if additional information would help. I’d be grateful for any insights into this behavior.

Best regards,
Xiao Lei

[EricKutschera]

Thanks for providing the files to troubleshoot

I was able to get rMATS v4.3.0 to detect that event using two reads and two transcript annotations from the files you shared

@HD	VN:1.0	SO:coordinate
@SQ	SN:chr1	LN:248956422
AV241602:4_7_2025_B_Xiao:2437668272:2:10605:1980:1444	83	chr1	45552429	60	135M15418N15M	=	45551087	-212	TGCCAAGGTTTCAAGTGATCCTCCCGCCTCAGCCTGCCCAGGTGCTGAGATTACATGTATGAGCCACTGCACCTGGAAAGGAGCCAGAAATGTGAAGTGCTAGCTGAAGGATGAGCAGCAGCTAGCCAGGCAAAGGCGGGGAGACAACCC	MLMMMMMMMMMMMNNLNNNNNNMN=NNNNNNNOOOOOOOOOONOOONONOONOOONOOONOOKOPOP@POOPPPPPPPPPOPPPPPNPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPOPPPPPPPPPPPPPPOPPPPMMLMMMMG	AS:i:-7	ZS:i:-11	XN:i:0	XM:i:1	XO:i:0	XG:i:0	NM:i:1	MD:Z:4C145	YS:i:0	YT:Z:CP	XS:A:+	NH:i:1
AV241602:4_7_2025_B_Xiao:2437668272:1:21804:0645:3001	163	chr1	45561816	60	63M6103N87M	=	45567979	210	CTGCACACTGGGCAGAAGATGCCTCTGATTGGTCTGGGTACCTGGAAGAGTGAGCCTGGTCAGGCGGGGAGACAACCCCTTCCCCAAGAATGCTGATGGGACTATATGCTACGACTCCACCCACTACAAGGAGACTTGGAAGGCTCTGGA	EEEFFFFFJGJJIJJMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMGHMMMMMMMMMMMMMLMMMMMMMMMMMMMMMMMMMMMLMMMMMJMMMMLMMKLMJMMMMMMMLMMMMKMMMMMMMMDMMM>M@MMMMMJMMMMIIMNML	AS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:150	YS:i:0	YT:Z:CP	XS:A:+	NH:i:1

chr1	StringTie	transcript	45550804	45570049	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; gene_name "AKR1A1";
chr1	StringTie	exon	45550804	45550876	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "1"; gene_name "AKR1A1";
chr1	StringTie	exon	45550966	45551155	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "2"; gene_name "AKR1A1";
chr1	StringTie	exon	45552436	45552563	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "3"; gene_name "AKR1A1";
chr1	StringTie	exon	45561789	45561878	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "4"; gene_name "AKR1A1";
chr1	StringTie	exon	45566569	45566688	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "5"; gene_name "AKR1A1";
chr1	StringTie	exon	45566869	45567020	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "6"; gene_name "AKR1A1";
chr1	StringTie	exon	45567982	45568177	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "7"; gene_name "AKR1A1";
chr1	StringTie	exon	45568485	45568684	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "8"; gene_name "AKR1A1";
chr1	StringTie	exon	45568927	45568999	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "9"; gene_name "AKR1A1";
chr1	StringTie	exon	45569143	45569229	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "10"; gene_name "AKR1A1";
chr1	StringTie	exon	45569891	45570049	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.3"; exon_number "11"; gene_name "AKR1A1";
chr1	StringTie	transcript	45550826	45570049	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; gene_name "AKR1A1";
chr1	StringTie	exon	45550826	45550876	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "1"; gene_name "AKR1A1";
chr1	StringTie	exon	45550966	45551155	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "2"; gene_name "AKR1A1";
chr1	StringTie	exon	45561789	45561878	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "3"; gene_name "AKR1A1";
chr1	StringTie	exon	45566569	45566688	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "4"; gene_name "AKR1A1";
chr1	StringTie	exon	45566869	45567020	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "5"; gene_name "AKR1A1";
chr1	StringTie	exon	45567982	45568177	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "6"; gene_name "AKR1A1";
chr1	StringTie	exon	45568927	45568999	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "7"; gene_name "AKR1A1";
chr1	StringTie	exon	45569143	45569229	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "8"; gene_name "AKR1A1";
chr1	StringTie	exon	45569891	45570049	1000	+	.	gene_id "ENSG00000117448.14"; transcript_id "MSTRG.991.6"; exon_number "9"; gene_name "AKR1A1";

I used this command: rmats.py --readLength 150 --allow-clipping --libType fr-firststrand --b1 ./b1.txt --gtf ./transcripts.gtf -t paired --nthread 1 --od ./out_dir --tmp ./tmp_dir --statoff

The 1-based coordinates for the exons in the event are:
upstream (45550966, 45551155)
1st (45552436, 45552563)
2nd (45561789, 45561878)
downstream (45567982, 45568177)

The transcript MSTRG.991.3 has upstream as exon 2, 1st as exon 3, 2nd as exon 4, and downstream as exon 7.
The transcript MSTRG.991.6 has upstream as exon 2, 2nd as exon 3, and downstream as exon 6.

As mentioned in this related issue (https://groups.google.com/g/rmats-user-group/c/aVdcBpcoJoQ) rMATS requires transcript evidence for each isoform in the MXE event (upstream->first->downstream, and upstream->second->downstream): https://github.com/Xinglab/rmats-turbo/blob/v4.3.0/rMATS_pipeline/rmatspipeline/rmatspipeline.pyx#L1546

There are two ways that rMATS can have transcript evidence for the 3 exon sequence. One way is to have a transcript in the --gtf that uses those 3 exons directly one after the other. The other way is from a supporting read from an alignment file. The read support for a 3 exon sequence is not from a read with 2 junctions, but from a single junction of a read combined with an annotated transcript: https://github.com/Xinglab/rmats-turbo/blob/v4.3.0/rMATS_pipeline/rmatspipeline/rmatspipeline.pyx#L381

For each individual junction in the read, rMATS will check the annotated transcripts in the gene for a transcript that uses both splice sites of the junction. If any transcript uses the exact junction then rMATS will not add a novel junction. If the exact junction isn't used, but some transcript(s) have the two splice sites on exons that are separated by 1 or more other exons then rMATS will add the novel junction and consider this read as supporting that transcript but with the exons that separated these two splice sites removed

In this example the GTF has upstream->1st->2nd->x->x->downstream and also a read for the 1st->downstream junction which is not in the GTF so rMATS uses that read as evidence for upstream->1st->downstream. Similarly, the GTF has upstream->2nd->x->x->downstream and also a read for the 2nd->downstream junction which is not in the GTF so rMATS uses that read as evidence for upstream->2nd->downstream

[santataRU]

Hi, Eric,

Thank you very much! You mentioned two reads as examples supporting the MXE events, based on my understanding of the counts:

IJC_SAMPLE_1	SJC_SAMPLE_1	IJC_SAMPLE_2	SJC_SAMPLE_2
107,93,129	146,152,188	161,136,158	101,92,86

there should be numerous reads spanning the same junctions as those example reads in each of the six samples (three per group). However, I did not observe such reads in certain samples. I will double-check this.

Best regards,
Xiao

[santataRU]

Hi Eric,

I did not observe any reads in the sample HT1080_shDDX42-8_3_VSVNL43_sorted_region_45550045_45571581.bam that support the junctions from the 1st->downstream exon or from the 2nd->downstream exon for the reported MXE event. However, the event still has junction counts in the output report, which seems unexpected.

Best regards,
Xiao

[EricKutschera]

rMATS first detects events (https://github.com/Xinglab/rmats-turbo/blob/v4.3.0/rMATS_pipeline/rmatspipeline/rmatspipeline.pyx#L4020) then it counts the reads that support the isoforms in those events (https://github.com/Xinglab/rmats-turbo/blob/v4.3.0/rMATS_pipeline/rmatspipeline/rmatspipeline.pyx#L4032)

In order to detect the MXE event it needs evidence for the two isoforms. Then the read counting is done as shown in the diagram from the README: https://github.com/Xinglab/rmats-turbo/tree/v4.3.0?tab=readme-ov-file#output

Running with --individual-counts gives these columns:

upstream_to_first_count    107,93, 129,161,136,158
first_to_downstream_count  0,  0,  0,  0,  0,  0
first_count                0,  0,  0,  0,  0,  0
upstream_to_second_count   145,152,188,101,92, 86
second_to_downstream_count 1,  0,  0,  0,  0,  0
second_count               0,  0,  0,  0,  0,  0

Almost all of the counts are for upstream_to_first or upstream_to_second

[santataRU]

Thank you very much, Eric.
Your analysis above indicates that there are no counts for the first_to_downstream junction. However, yesterday you identified a read that supports the 1st → downstream junction (#512 (comment)). I was wondering how to explain this discrepancy?

[EricKutschera]

The alignment starting at 45552429 with the CIGAR 135M15418N15M is for these coordinates:
(45552429, 45552563), (45567982, 45567996)

The junction is from the end of the 1st exon to the start of the downstream exon. The read ends within the downstream exon, but starts 7bp before the 1st exon start (45552429 compared to 45552436). The read has the required splice sites for event detection, but the read start coordinate prevents it from counting after event detection

The code I linked to before, that uses a read to support a 3 exon sequence for MXE detection is only checking the splice sites against the annotated transcripts. The code for counting a read toward the 1st to downstream junction of a detected MXE event also checks the 1st exon start coordinate: https://github.com/Xinglab/rmats-turbo/blob/v4.3.0/rMATS_pipeline/rmatspipeline/rmatspipeline.pyx#L2049

[santataRU]

Thank you very much, Eric.

Based on discussions in the rMATS Google forum, I understand that identifying MXE events can be challenging. I’ve been thinking that a simple yet effective rule might help improve accuracy in certain cases. For example, in this instance, using the IGV browser, I can clearly observe that the junction counts supporting splicing from the 1st exon to the 2nd exon are much higher than those supporting 1st → downstream or 2nd → downstream splicing.

By definition, mutually exclusive exons cannot co-exist. Therefore, a possible criterion to validate an MXE event could be: the junction counts supporting 1st exon → 2nd exon must be lower than at least one of the following:

1. upstream_to_first_count

2. first_to_downstream_count

3. upstream_to_second_count

4. second_to_downstream_count

Based on this rule, the event in question would not qualify as a true MXE event.

Interestingly, rMATS did successfully identify the 1st exon as a skipped exon (SE) event with significance, which seems better supported by the data. Here's the corresponding SE record:

ID	GeneID	geneSymbol	chr	strand	exonStart_0base	exonEnd	upstreamES	upstreamEE	downstreamES	downstreamEE	ID	IJC_SAMPLE_1	SJC_SAMPLE_1	IJC_SAMPLE_2	SJC_SAMPLE_2	IncFormLen	SkipFormLen	PValue	FDR	IncLevel1	IncLevel2	IncLevelDifference
28162	ENSG00000117448.14	AKR1A1	chr1	+	45552435	45552563	45550542	45551155	45561788	45561878	28162	179,134,186	342,324,353	255,234,256	219,212,197	277	149	0	0	0.22,0.182,0.221	0.385,0.373,0.411	-0.182

This SE event is well supported by read data and appears to be the more biologically plausible interpretation.

What do you think?

Best regards,
Xiao

[EricKutschera]

Your suggestion to validate a MXE event based on the different counts seems reasonable. rMATS doesn't report the 1st exon to 2nd exon junction count for a MXE event. In this specific example you could get that count from the SE event that uses that junction. Maybe in the future the code could be updated to output that value for all MXE events. For now you could run rMATS with a custom --fixed-event-set such that the output would have all of the junction counts you want to use to validate any MXE events of interest

[santataRU]

Thank you very much, Eric.

I've observed that significant MXE events include the same 1st and/or 2nd exons involved in significant SE events. Recently, I found that many MXE calls are likely not well supported by reads.

For example, I used rMATS to analyze RNA-seq data comparing control (shNT) to knockdown (shDDX42) using two different short hairpin RNAs (shDDX42_2 and shDDX42_8). When comparing the overlap of splicing events between the two knockdown conditions, I noticed that SE, RI, A5SS, and A3SS events showed good agreement between the two shRNA treatments. In contrast, the MXE events showed poor overlap (see attached image). I performed a filter first so that all events in the overlap analysis contain at least 3 counts for each event in each sample.

One likely explanation is that many of the identified “MXE” events do not meet the stricter criteria we discussed before and may be incorrectly annotated as MXEs. If rMATS could incorporate a filter to remove these questionable cases, I suspect the overlap of MXE events between the two shRNA treatments would improve substantially.