Question: Calculate reads and paired inserts covering a junction
0
gravatar for craigmonger
3.8 years ago by
craigmonger10
craigmonger10 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 COMMENTlink modified 3.8 years ago by lakhujanivijay5.4k • written 3.8 years ago by craigmonger10
3
gravatar for Pierre Lindenbaum
3.8 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum134k 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 COMMENTlink written 3.8 years ago by Pierre Lindenbaum134k

Thanks Pierre for your answer and tools!

ADD REPLYlink written 3.8 years ago by craigmonger10

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

ADD REPLYlink written 3.8 years ago by Pierre Lindenbaum134k
1

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 REPLYlink written 3.8 years ago by craigmonger10
0
gravatar for Rohit
3.8 years ago by
Rohit1.4k
California
Rohit1.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 COMMENTlink written 3.8 years ago by Rohit1.4k
0
gravatar for lakhujanivijay
3.8 years ago by
lakhujanivijay5.4k
India/Ahmedabad
lakhujanivijay5.4k wrote:

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

enter image description here

ADD COMMENTlink written 3.8 years ago by lakhujanivijay5.4k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1567 users visited in the last hour
_