This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Get Junctions File Using Gsnap

Hi Everyone

So I am working on rna seq data. I am using gsnap for my analysis(though I have also used tophat). With tophat , a junctions.bed file is produced which actually tells us about the splicing patterns going on.

Using gsnap I am able to produce sam/bam files. Is there a way I can convert these sam/bam files into something similar to junctions.bed files.

My main goal is to see what is the difference between junctions produced from tophat and gsnap. Tophat produces junctions.bed file but gsnap doesnt. So is there a way I can get a junctions file from sam/bam output of gsnap.

Hope to hear from you guys soon

Regards
Varun

gsnap rna-seq
• 4,287 views • link
updated by Ram Ram • written by Varun Gupta Varun Gupta

Are you comfortable writing perl scripts?

• 0 views • link
written by Arun Arun

great, then it is relatively straightforward to obtain the junctions given a SAM/BAM file. All you need is the start position of a read and its CIGAR. Suppose a read starts at 10000 and cigar is 15M20N75M, then the read is present from 10000-10014 and 10035-10109. So, the read's junction coordinates are 10015-10034. It should be a bit more involved to take into account the effect of Insertions and deletions. But this is the basic principle. Let me know if you run into problems.

• 1 views • link
written by Arun Arun

Hey Arun

What if the CIGAR string is something like this 1S14M20N75M

• 0 views • link
written by Varun Gupta Varun Gupta

1 answer

If you're comfortable with R/BioConductor

and have a current installation

then

you can the the approach given in my recent update to Obtaining Junctions.Bed File From Aligners Such As Gsnap

• 0 views • link
written by Malcolm.Cook Malcolm.Cook

Log in to answer this question.

More posts like this