Calculate reads and paired inserts covering a junction

Question

[craigmonger]

Hi all,

Is there a way to calculate the number of read alignments which cover an exon junction, as well as the number of read pairs which align to flanking exons whose insert covers the junction?

E.G. in the following example the coverage would be 3.

Is this something any existing tool can calculate?

Regards,

Craig

Answer 1

[Pierre Lindenbaum]

Using samjs:

the first exon ends at "rotavirus:100", the second exon starts at "rotavirus:125"

• the first samtools get the reads in that region
• samjs remove the unmapped reads or the reads on a bad contig, we acceot the reads starting before exon1_end and ending after exon2_start. If the reads are paired we accept the reads where the left segment is before exon1_end and the right segment is after exon2_start.
• the last samtools count the reads

    $ samtools view -b S1.bam "rotavirus:100-1000" |\
     java -jar dist/samjs.jar -e 'function accept(rec){var exon1_end=100,exon2_start=125;if(rec.getReadUnmappedFlag() ) return false; if(rec.getReferenceName()!="rotavirus") return false; if(rec.getAlignmentStart() <= exon1_end && rec.getAlignmentEnd()>=exon2_start) return true; if(rec.getReadPairedFlag() && !rec.getMateUnmappedFlag() && Math.min(rec.getAlignmentStart(),rec.getMateAlignmentStart()) <=exon1_end && Math.max(rec.getAlignmentEnd(),rec.getMateAlignmentStart())>=exon2_start ) return true; return false;}accept(record)'  | \
    samtools view -c

[craigmonger]

Thanks Pierre for your answer and tools!

[Pierre Lindenbaum]

less than one minute between my answer and the accepted status. Thanks, but I can't imagine you have just tested this.

[craigmonger]

Apologies, I should have tested before taking your word for it that it would work.

I have now tested the proposed answer and it seems to have worked.

Thank you once again.

Answer 2

[Rohit]

If you have all the junctions in a bed file, and by flanking exons covering the junction you expect +/-50 bp for each junction, you cancalculate coverage using bedtools and try to overlap your coordinates with the exon junctions modified to include flanking regions.

Answer 3

[lakhujanivijay]

You can try HTSeq. Check-out the different modes available.