Thanks a lot for your answer !!
I would like to extract intron-exon reads (unspliced reads) from bam file. I tried use the with following command command, bedtools intersect -s -F 1 -split -a input.bam -b junction.bed >Intron_out.bam
Where junction.bed consists of last bp intron and first bp exonic. It doesn't to seem to be working, because, I am also getting the split reads that I am not expecting.
I would really appreciate, if someone can help me with the same.
2 answers
I would think that a multi-step process would work:
- Filter for reads overlapping an exon
- Filter those passing step 1 for those overlapping an intron
- Filter those above for more than 1 CIGAR operation (you could filter out spliced alignments instead, but there's not then a simple
samtools viewoption...this is quick and easy and close enough).
Hi Devon, Could you please tell me the best way to filtered the reads overlapping an exon. I was trying to bedtools intersect with (v) options. Is it correct in your opinion ? My Best Regards, Prasoon
You're not going to want the -v option to select for reads overlapping exons.
Sorry for the confusing language, I want the reads overlapping intron-exon. Therefore, My question about the step 1 that you have suggested. How can I filter reads overlapping exon?
bedtools intersect with a BED file of exons for step 1. bedtools intersect with a BED file of introns for step 2.
I'm not sure if I understand well, your problem: you want the cigar section ('M' or 'X' or '=') overlapping the BED file isn't it.
Here is a solution using http://lindenb.github.io/jvarkit/SamJdk.html
private List<Interval> intervals=null;
public Object apply(final SAMRecord record) {
if(this.intervals==null)
{
/** first call: we load the BED file "junction.bed" */
try
{
java.io.BufferedReader r=new java.io.BufferedReader(new java.io.FileReader("junction.bed"));
this.intervals = r.lines().
filter(S->!(S.isEmpty() || S.startsWith("track")|| S.startsWith("browser"))).
map(S->{
final String tokens[]=S.split("[\t]");
return new Interval(tokens[0],1+Integer.parseInt(tokens[1]),Integer.parseInt(tokens[2]));
}).
collect(Collectors.toList());
r.close();
}
catch(Exception err)
{
throw new RuntimeException(err);
}
}
// ignore unmapped
if(record.getReadUnmappedFlag()) return false;
final Cigar cigar = record.getCigar();
if(cigar==null || cigar.isEmpty()) return false;
// loop over the cigar string
int pos = record.getStart();
for(final CigarElement ce: cigar)
{
final CigarOperator op= ce.getOperator();
if(op.isAlignment()) // it's 'M' or 'X' or '='
{
// search overlapping our bed list
for(final Interval r:this.intervals)
{
//found !
if(r.getContig().equals(record.getContig()) &&
pos <=r.getStart() &&
(pos+ce.getLength()) >= r.getEnd())
{
// accept this sam record
return true;
}
}
}
if(op.consumesReferenceBases())
{
pos+=ce.getLength();
}
}
// we don't accept
return false;
}
usage:
java -jar dist/samjdk.jar -f biostars292561.code --body input.sam
Thanks a lot for your answer.I would like to extract the reads transverse of intron-exon.
so , my solution should work.
P.
Hi Pierre, Thanks a lot. this program seems to be quite slow. Can you please suggest me if I can do something to speed up ? My Best Regads, Prasoon
sure, you can reduce the input bam using samtools view and piping it into my tool. Something like:
samtools view -bu your.bam "chrom:1234-456" |java -jar dist/samjdk.jar -f biostars292561.code --body
or
samtools view -bu -L intervals.bed your.bam |java -jar dist/samjdk.jar -f biostars292561.code --body
Log in to answer this question.