Event observed in sashimi, not declared in rmats output

[neuroRNA]

Hi everyone!

I wanted to share an issue I’ve encountered while using rMATS.

 

I am interested in the alternative splicing of exon 7 in SMN2 gene (coordinates hg38 chr5:70076520-70076574). In the bam generated when running rmats from FASTQs with default parameters I can clearly see that the exon under study is more included in our contrast group (SMN) than in the control (NT). * See attached image. 

 

However, in the rMATS output, this event does not appear as differentially spliced. When I check the output file (SE.MATS.JC.txt), my exon of interest has the following results:

 

ID        .       GeneID                       geneSymbol        chr     strand        exonStart_0base exonEnd        upstreamES        upstreamEE        downstreamES        downstreamEE        

95692        "ENSG00000205571.14"        "SMN2"         chr5        +        70076520        70076574        70070640        70070751        70077018        70077160

 

IJC_SAMPLE_1        SJC_SAMPLE_1        IJC_SAMPLE_2        SJC_SAMPLE_2        IncFormLen        SkipFormLen        PValue        FDR        IncLevel1        IncLevel2        IncLevelDifference

178,197,320        0,5,1                              164,131,189        5,3,9        203                              149        1        1.0        1.0,0.967,0.996          0.96,0.97,0.939        0.031

 

Here, sample 1 corresponds to the SMN group. 

 

Do you have any explanation why this exon is so cleary differentially spliced according to the bam (and indeed validated by qPCR in the same samples), but it doesn't appear as differentially spliced in rmats output?

 

 I would greatly appreciate any advice or suggestions.

 

Thanks!

 

María

[kutsc...@gmail.com]

From the sashimiplot it looks like there are around 80 reads for the skipping junction in the 1st sample, but the rMATS output has 0 for the 1st sample in SJC_SAMPLE_1. The issue is that rMATS is not counting most of the reads for the skipping junction. This post describes the differences in counting between rMATS and rmats2sashimiplot: https://github.com/Xinglab/rmats2sashimiplot/issues/33#issuecomment-656720560

In this case I think the issue is multimapped reads. rmats2sashimiplot will count multimapped reads, but rMATS requires reads to be uniquely mapped (NH=1). The sequence around the skipping junction may be similar for SMN1 and SMN2 and result in multimapping

You could edit your bam files to always set the NH tag to 1 and then rMATS would use all the reads. That could lead to double counting some reads. Another option is to filter for the primary alignments when setting NH to 1: https://github.com/Xinglab/rmats-turbo/issues/298#issuecomment-1568412497

Eric

[neuroRNA]

Thanks so much, Eric, for the detailed response and the links! 

The explanation of how rmats2sashimi counts reads is really clarifying.  We'll try the recommended advice of editing the BAM files to set the NH tag to 1 and running rMATS again. 

Really appreciate the help!

María