1
$\begingroup$

I have 3 biological replicates of RNA seq data for a particular condition. I want to find out intron retention events from those biological replicates for a given condition. There is no comparison I am doing here. Just from the three replicates, I am trying to find out what introns are retained in my sample. There is only one condition here. Tools like rMATS require 2 conditions. Is there any tool I can use to find out intron retention events given that I only have one condition (with replicates).

Ideally I would like to use a method similar to that used in this paper, though I can't fully grasp their methodology. This would allow me to get p-values for enrichment of specific introns.

$\endgroup$
0

3 Answers 3

1
$\begingroup$

All introns will have some level of retention and I guess no intron will be 100% retained. In the absence of a comparator, you will need to work out how much intron retention an intron needs to be called retained. This is a human judgement, there is no automated way of deciding this.

What you need to do this is some sort of PSI calculation (percent-spliced-in). MISO will provide this, but it doesn't do replicates, so I guess you'll need to run it on each sample and then average the PSIs and determine if it surpasses some sort of threshold you have set.

$\endgroup$
1
  • 1
    $\begingroup$ I second this, but would suggest using SUPPA which handles replicates out of the box and is actively developed. As for deciding which introns are retained, I expect that most biological conditions most of them will have a low PSI (which a nice, easily interpretable metric), so I have in the past plotted the distribution and used some upper bound of the distribution to be on the safe side. That said, some condition to use as a reference would be better. $\endgroup$ Jan 22, 2018 at 8:46
0
$\begingroup$

The method you linked to is interesting, though I suspect that Ian's answer is both easier and generally sufficient to get what you want.

Having said that, the general idea behind the method you linked to is to use DESeq2 to test for differentially expressed introns. To generate the control samples, they're effectively taking the number of reads across all introns in a gene and distributing them over the introns. They weight the distribution such that longer introns receive more reads (they use an effective length for this, but I'm skeptical that this makes that much of a difference). DESeq2 is then giving the p-values you're after. Unless you can get the authors to send you there code, you'll have to write a bit of code to do all of this...but it shouldn't be too bad.

BTW, you can use the GenomicFeatures package from bioconductor, specifically the exonsBy(..., by="gene) follows by gaps().

$\endgroup$
7
  • $\begingroup$ Hi Devon, Thanks for your help.The paper compares insilco replicate counts with actual counts and that is what is input for DESeq2. From your explanation, for a gene, they took all the reads which map introns and then lets say we get 1000 reads for all the possible introns for this gene and then for each intron, they give a weightage like 1000/(len of intron). I am looking at GenomicFeatures and I have already created the intron file coordinates which you mentioned at seqanswers. If possible and if you have time can you guide me in writing the above code. Thanks for your help always Regards $\endgroup$ Jan 9, 2018 at 20:59
  • $\begingroup$ It's not 1000/(len of intron), but rather 1000/(len of intron / sum of all intron lengths in gene). $\endgroup$
    – Devon Ryan
    Jan 9, 2018 at 21:10
  • $\begingroup$ I only have time to provide vague guidance on the code. You can use findOverlaps() (from GenomicRanges in bioconductor) for the counts. Aside from that it's mostly split() and mendoapply() that you'll need. $\endgroup$
    – Devon Ryan
    Jan 9, 2018 at 21:13
  • $\begingroup$ so that would give us insilco replicate1 count for a particular intron and similarly insilco replicate2 count for 2nd replicate and then we compare these 2 with actual counts we are getting in rep1 and rep2 via DESeq2 $\endgroup$ Jan 9, 2018 at 21:15
  • 1
    $\begingroup$ D'oh! I've actaully done this before (that'll teach me not to look at the links in the question). I used feature counts to count the reads in introns, then a python script to calculated the effective length from mappability tracks downloaded from UCSC, and then used an R script to find the "detained introns". See these gists: gist.github.com/IanSudbery/df92ebf4f84d76027394750a6516db1b , gist.github.com/IanSudbery/eb36ba5d59b4455f9cd7b2085447eb3d. I offer no guarantees with these, they were pulled from my archives. $\endgroup$ Jan 11, 2018 at 11:14
0
$\begingroup$

IRFinder can be used to detect the IR events and can easily reach your goal.

Here is an example. IRFinder can calculate the IR ratio to detect IR event. mimirna.centenary.org.au/irfinder/example1.html. (You just need run IRFinder for one condition and each replicate will output a folder). Two conditions can be used to detect the differential IR events. If you have only one condition, you can select the IR events from the common events of the three replicates or select IR events exist at least two of three replicates based on the main output files( named IRFinder-IR-dir.txt). The IR events can be selected by filtering the IR ratio, Intron coverage etc

$\endgroup$
1
  • 2
    $\begingroup$ Hi Renhua, this answer doesn't seem to really address the issue. Could you please elaborate how IRFinder helps when there is one condition with replicates? How would that be done with IRFinder? Many thanks for your time $\endgroup$
    – llrs
    Jan 22, 2018 at 9:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.