Technical Question – rMATS with Genotype-Specific Custom References

[Odudu James]

Dear rMATS development team,

I am writing to seek technical guidance on whether rMATS is appropriate for my experimental setup and, if so, how to best configure it.

**Experimental context:**
I am studying alternative splicing at intron 1 of TP53 across 4 genotypes:
- Alu_free (reference level for comparisons)
- H_sense (transgene inserted in sense orientation)
- H_antisense (transgene inserted in antisense orientation)
- H_intron (wild type)

Each genotype has 3 biological replicates (H_sense has 4), giving 13 BAM files total. All BAM files were generated using STAR alignment.

**The core issue:**
Each genotype was aligned to its own custom reference genome and GTF file. The custom references were built by taking a ~90 kb region of chromosome 17 encompassing TP53 plus flanking genes, reverse complementing it into a mini-chromosome (since TP53 is on the antisense strand), and soft-masking the remainder of the reference annotation. In the non-wild type genotypes, intron 1 of TP53 was replaced with a modified sequence (transgene), which changes the length of intron 1 and shifts all downstream exon coordinates by a genotype-specific offset. This means each genotype's BAM files exist in a different coordinate space.

**My questions:**
1. Can rMATS handle pairwise comparisons where the two groups were aligned to different reference genomes with different coordinate systems? If not, is remapping all samples to a single unified reference the only solution?

2. If a unified reference is required, what would you recommend as the best strategy for constructing it given that the genotypes differ only within intron 1 of TP53?

3. Would rMATS be able to detect novel splice sites created by the transgene insertion, or is it limited to splice events present in the provided GTF?

4. Is there a recommended way to restrict the rMATS analysis to a specific genomic region (intron 1 of TP53 and its flanking exons) rather than running it genome-wide?

Thank you very much for your time. I am happy to provide additional details about the reference construction or file structure if helpful.

[kutsc...@gmail.com]

rMATS assumes a single coordinate system. I don't think there is an easy way for rMATS itself to handle two groups where the genome region of interest has a large change. If you were to align everything to the same reference sequence then I expect that many of the reads would be aligned with insertions or deletions or fail to align. rMATS actually filters out alignments with insertions or deletions

rMATS can detect novel splice sites if it's run with --novelSS

You can restrict rMATS to a specific region by using a --gtf that only includes the exons in the region of interest. Any alignments which don't overlap an exon from the --gtf will be filtered out

One thing you could try is processing each group separately using the specific reference genome and GTF for that group for both aligning the reads and running rMATS. Then you could attempt to match up events from the output files like SE.MATS.JC.txt based on some translation of the coordinates. For each reference genome you could translate the coordinate columns in the rMATS output files to be offsets from the intron that differs. The coordinate number could be converted to something like {offset}_before, {offset}_within, or {offset}_after depending on whether the coordinate is before, within, or after the intron

If you do manage to find the same event in two groups then you can run the rMATS statistical model by using the count columns like IJC_SAMPLE_1 and the isoform length columns like IncFormLen. You can run rMATSexe directly on a file with the counts like in this post: https://groups.google.com/g/rmats-user-group/c/2PJ6DWFu1m8/m/0J0eY3XlAAAJ

Eric

[kutsc...@gmail.com]

The statistical model needs counts from two groups so it only makes sense to run it on matched events. You won't be able to use cp_with_prefix.py or prepare_stat_inputs.py if you ran the groups separately with each group's specific genome

You can run --task stat if you create fromGTF.SE.txt, JC.raw.input.SE.txt, and JCEC.raw.input.SE.txt files (also for MXE,A3SS,A5SS,RI). When you match an event between two groups you can add it to your new fromGTF file. The event ID likely won't match up, but you can number the new file starting from 0. For each matched event, the corresponding new JC.raw.input.SE.txt and JCEC.raw.input.SE.txt files can use IJC_SAMPLE_1 and SJC_SAMPLE_1 from one group directly, and then use the values from the other group as the _SAMPLE_2 columns

Eric

[Odudu James]

Subject: Follow-up: rMATS Developer Approach Implementation for TP53 Intron 1 Analysis

Good afternoon Eric,

I hope you are doing well. I would like to follow up on the recommendation you provided concerning the analysis of my dataset. I am attaching the scripts I implemented based on your guidance (please note the .py files will need to be viewed with a text editor such as Notepad).

Here is my understanding of what was implemented at each stage:

Stage 1:
Each genotype was processed independently by rMATS using only its own BAM files and matching GTF. The --statoff flag was used because no statistical comparison is made at this stage; only read counting and event detection are performed.

Stage 2:
Because intron 1 lengths differ between genotypes, the same biological event has different raw coordinates in each genotype's output. All coordinates were translated into a standardized framework using the known intron 1 boundaries and offsets relative to H_Intron (wild type). Coordinates were converted into three standardized labels: upstream_{position} for coordinates before intron 1 start, intron1_rel_{offset} for coordinates within intron 1, and downstream_wt_{position} for coordinates after intron 1 end expressed in H_Intron coordinate space.

Stage 3:
For each pairwise comparison and event type, an inner join was performed on the standardized coordinate columns to find matching events between the two genotypes. This produced the three input files required by --task stat: fromGTF.{event}.txt, JC.raw.input.{event}.txt, and JCEC.raw.input.{event}.txt. The experimental genotype's counts were used as SAMPLE_1 and the reference genotype's counts as SAMPLE_2, as per your recommendation.

Stage 4:
run_rmats --task stat was run on each comparison directory containing the matched input files built in Stage 3.

In terms of matched events, Stage 3 produced the following:
- Comparison 1 (H_Sense vs Alu_free): ~718 matched events
- Comparison 2 (H_Antisense vs Alu_free): ~709 matched events
- Comparison 3 (H_Antisense vs H_Sense): ~261,220 matched events (both genotypes share identical intron 1 length)
- Comparison 4 (Alu_free vs H_Intron): ~715 matched events

I have a few specific questions I would appreciate your feedback on:

1. Does the coordinate translation in Stage 2 match what you had in mind, particularly the downstream_wt conversion using the intron 1 offset?
2. Does the event matching approach in Stage 3 via inner join on standardized coordinate columns look correct?
3. We noticed a large FDR discrepancy in Comparison 4 between this approach and a direct pairwise run. Several events that are highly significant in the direct approach (FDR as low as 5.37e-11) show FDR = 1.0 in the Developer's Approach. We believe this is due to coordinate artifacts in the direct approach since the H_Intron intron 1 end (position 40,868) has no equivalent annotation in Alu_free. Does this interpretation make sense?

Thank you very much for your time and I look forward to your feedback.

Sincerely,
Odudu

--
You received this message because you are subscribed to the Google Groups "rMATS User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rmats-user-gro...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/rmats-user-group/8cbe3d53-1a5d-4dda-886c-ee396b0a4dfan%40googlegroups.com.

[kutsc...@gmail.com]

1. The code for the coordinate translation in Stage 2 looks good to me
2. The event matching code seems fine until the de-duplication. I don't think there should be multiple events in a single rMATS run with the same key. The join should only be able to find at most 1 match per key when checking a pair of files. Can you find an example of events with the same key in one of the original rMATS output files?
3. How did you do a direct pairwise comparison with samples that have different genome coordinates?

The matched event counts seem strange. Comparisons 1,2,4 have around 700 but comparison 3 which has the same intron length has 261k. It might be that the comparisons with different coordinate systems are only matching events before the intron region and the other comparison is able to match events at any position. You could try to manually find an event downstream of the intron in one of the other comparisons where the event exists for both samples in the untranslated output files. Also, 261k seems like more events than just for a single chromosome. Are there events for other chrs and is the coordinate translation being applied to other chrs?

Eric