Dear group, 
Given RNA-Seq data for either 1 or more than 1 different samples/conditions, I am interested in finding 

A. an exon-intron or intron-exon junction with reads spanning these junctions.  
B. Another similar question is identify intron-retention, where reads map exon-intron-exon. 

A differs from B, where A coverage ends halfway in intron. There are algorithms such as iREAD for identifying intron-retention similar to B above. 

Is there a way, say samtools or other tool where I can find intron-exon or exon-intron read overhang extension.  

As an example, image is given for easy understanding. Thanks for your help. 

-A



Image showing intron-exon read overhang for top sample, whereas for the rest of samples below no such read overhang or coverage is observed.

![enter image description here][1]


  [1]: http://funkyimg.com/i/2JmgG.jpg
Steps via `awk` and [BEDOPS][1].

Build a BED file of exons:

    $ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gtf.gz \
        | gunzip -c - \
        | awk '($3=="exon")' - \
        | gtf2bed - \
        | cut -f1-6 - \
        > gencode.v28.annotation.exons.c1t6.bed

Convert transcripts to merged exons:

    $ awk -f transcripts2mergedExons.awk gencode.v28.annotation.exons.c1t6.bed > gencode.v28.annotation.mergedExons.c1t6.bed

Make an exon-intron list:

    $ awk -f mergedExons2exonIntronList.awk gencode.v28.annotation.mergedExons.c1t6.bed > gencode.v28.annotation.exonsAndIntrons.c1t6.bed

Convert these to junctions:

    $ awk -f exonIntronList2JunctionList.awk gencode.v28.annotation.exonsAndIntrons.c1t6.bed > gencode.v28.annotation.exonsIntronJunctions.c1t6.bed

Optionally, pad them, e.g. by 25 nt around the junction (adjust as needed):

    $ bedops --everything --range 25 gencode.v28.annotation.exonsIntronJunctions.c1t6.bed > gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed

Map reads to the padded junctions:

    $ bedmap --echo --count --delim '\t' gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed <(bam2bed < reads.bam) > answer.bed

The file `answer.bed` will contain the junction and the number of reads that map to — overlap with — the optionally-padded junction, by one or more bases.

If you want the actual reads that map to the exon-intron junctions:

    $ bedmap --echo-map --delim '\t' gencode.v28.annotation.exonsIntronJunctions.pad25.c1t6.bed <(bam2bed < reads.bam) | sort-bed - > answer.bed

In this second example, the file `answer.bed` will contain the reads that map to the junction, by one or more bases.

Some links to Github Gists with listed `awk` scripts:

* (`transcripts2mergedExons.awk`) https://gist.github.com/alexpreynolds/7f68953db11c7c00569141b577c01c81 
* (`mergedExons2exonIntronList.awk`) https://gist.github.com/alexpreynolds/431216b2c3381b476bf91b149d35b2cc 
* (`exonIntronList2JunctionList.awk`) https://gist.github.com/alexpreynolds/72c8deb54e85e4028698563d8c49f470 


  [1]: https://github.com/bedops/bedops