Calculate Percent Spliced In (PSI)

Question

[Vikas Bansal]

Dear all,

I was wondering if there is any tool to calculate PSI from the BAM files. Basically what I want is-

I have gtf file with all the exons and I want to calculate for each exon (lets take exon filled with grey for example) how many reads mapped to junction of that exon (a and b in figure) and how many reads skip that exon (c in figure). So if we have, let's say, 10k exons in the gtf, the output should contain 10k rows with a,b and c as columns.

I used MATS and bam2ssj but did not get what I was expecting. I have gtf file with approx. 300k exons but MATS gives the output only for 3k exons. bam2ssj is calculating the values for introns rather than exons.

Any suggestions?

Thanks in advance.

[Martombo]

have a look at MISO: http://genes.mit.edu/burgelab/miso/docs/

[Vikas Bansal]

Thanks for the reply. I read about MISO but it says

Two general kinds of analyses are possible:

1. Estimate expression level of exons ("exon-centric" analysis) or,
2. Estimate expression level of whole transcripts ("isoform-centric" analysis)

I think it include the reads mapped on exons but I am interested only in junction reads.

[Martombo]

in the documentation there is also a paragraph about psi estimation.

the definition given in their paper is: 'percentage spliced in' (PSI or Ψ)11 denotes the fraction of mRNAs that represent the inclusion isoform. Reads aligning to the alternative exon or to its junctions with adjacent constitutive exons provide support for the inclusion isoform, whereas reads aligning to the junction between the adjacent constitutive exons support the exclusion isoform; the relative read density of these two sets forms the standard estimate of Ψ.

the paper is here

[derseb]

Hi Vikas,

we have just written a small protocol to do this:

http://onlinelibrary.wiley.com/doi/10.1002/0471142905.hg1116s87/abstract

Answer 1

[Pierre Lindenbaum]

I just tried to write something for this: https://github.com/lindenb/jvarkit/wiki/Biostar103303

I don't have any data to test though.

[Vikas Bansal]

Thanks a lot for the reply. Its amazing that you included it in jvarkit so quickly. I am trying to install jvarkit but it is not working. Below are my commands -

[Pierre Lindenbaum]

use jdk7 please, not 8 (java 8 is too "new" for the libraries included in htsjdk)

[Vikas Bansal]

Thanks. Now build is fine. But when I do ant biostar103303 or ant Biostar103303 or ant Biostar103303.java, it says

Buildfile: /project/specht/tools/jvarkit/htsjdk/build.xml

BUILD FAILED
Target "biostar103303" does not exist in the project "htsjdk".

Total time: 0 seconds

[Pierre Lindenbaum]

you forgot to cd ... Please run ant biostar103303 in varkit not in varkit/htsjdk

[Vikas Bansal]

Finally, I got the output file. But the columns - "exon.count_prev_and_next", "exon.count_prev_and_curr" "exon.count_curr_and_next" are all 0. I think it is not able to calculate junction reads.

[Pierre Lindenbaum]

do you have a bam (a subset) and a link to gtf so I can test my program ?

[Vikas Bansal]

Here is the link for the BAM (subset) file - https://www.dropbox.com/s/k2ibkw81t4hexq8/onlymyh7.bam

and gtf I downloaded from here - http://genome.ucsc.edu/cgi-bin/hgTables?hgsid=378435437_BVHytx2vRGH3XVmRsp39JG5acTeD&clade=mammal&org=Mouse&db=mm9&hgta_group=genes&hgta_track=refGene&hgta_table=0&hgta_regionType=genome&position=chr12%3A57795963-57815592&hgta_outputType=primaryTable&hgta_outFileName=

Please let me know if you need any other file.

[Pierre Lindenbaum]

thanks, I'll play with it next week.

[Pierre Lindenbaum]

I've updated my sources.It seems to work now. Re-download everything to use git to update the sources. Please use https://github.com/lindenb/jvarkit/issues?state=open to report others bug.

[Vikas Bansal]

Thanks. I think I found a bug and I will post it on github. Just one thing- Did you publish jvarkit? I found it quite fast. Although I never use java, I found it easy to use.

[dnaJNKY]

Will you explain what the columns represent in the results.tsv file?