Calculate reads and paired inserts covering a junction

Question: Calculate reads and paired inserts covering a junction

3.8 years ago by

craigmonger • 10 wrote:

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.

example

Is this something any existing tool can calculate?

Regards,

Craig

rna-seq • 1.3k views

ADD COMMENT • link •

modified 3.8 years ago by lakhujanivijay ♦ 5.4k • written 3.8 years ago by craigmonger • 10

3.8 years ago by

Pierre Lindenbaum ♦ 134k

France/Nantes/Institut du Thorax - INSERM UMR1087

Pierre Lindenbaum ♦ 134k wrote:

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

ADD COMMENT • link written 3.8 years ago by Pierre Lindenbaum ♦ 134k

Thanks Pierre for your answer and tools!

ADD REPLY • link written 3.8 years ago by craigmonger • 10

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

ADD REPLY • link written 3.8 years ago by Pierre Lindenbaum ♦ 134k

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.

ADD REPLY • link written 3.8 years ago by craigmonger • 10

3.8 years ago by

Rohit • 1.4k

California

Rohit • 1.4k wrote:

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.

ADD COMMENT • link written 3.8 years ago by Rohit • 1.4k

3.8 years ago by

lakhujanivijay ♦ 5.4k

India/Ahmedabad

lakhujanivijay ♦ 5.4k wrote:

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

enter image description here

ADD COMMENT • link written 3.8 years ago by lakhujanivijay ♦ 5.4k

Please log in to add an answer.

Similar posts • Search »