Calculating p-value for introns which are retained (or not) between 2 conditions
Here is how we do this:

First break your gene model down into non-overlapping chunks, such that any possible gene model can be built from a combination of those chunks. E.g.

transcript 1 |&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;|----------|&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;|---|&gt;&gt;&gt;&gt;|
transcript 2 |&gt;&gt;&gt;&gt;&gt;&gt;|-------------|&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;|---|&gt;&gt;&gt;&gt;|
transcript 3 |&gt;&gt;&gt;&gt;&gt;&gt;|-------------|&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;&gt;|
chunks       |+++1++|+2|+++++3++++|+++4++++|+5+|+6++|


this can be achieved using cgat gtf2gtf --method=genes-to-unique-chunks from the cgat-apps toolkit if you don't want to code it yourself.

I call
chunks 1, 4 and 6 constitutive exon
chunk  3 constitutive intron
chunk  2 alternate intron
chunk  5 annotated retained intron.

Now use featureCounts to count the number of reads in each chunk.

Test the chunk counts for usage using DEXSeq. 

Filter the results to just include chunks of type 2, 3 and 5 and recompute the FDRs.