Thanks Pierre for your answer and tools!
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
3 answers
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
less than one minute between my answer and the accepted status. Thanks, but I can't imagine you have just tested this.
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.
Log in to answer this question.
