This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to extract reads with no INDEL?

Hi, all. I am using samtools view .bam [region] to extract all reads overlapping the region. I find some reads with INDEL and I want to filter them, only getting the reads with no INDEL. Which samtools options I should use?

samtools indel sam
• 6,518 views • link
updated by dariober dariober • written by Jeason Rad • Jeason Rad

2 answers

Isn't it enough to keep reads without I or D operator in the cigar? For this you can use awk like:

samtools view -h in.bam chr1:12340-200000 \
| awk '$1 ~ "^@" || $6 !~ "I|D"' \
| samtools view -b - > out.bam
• 0 views • link
written by dariober dariober

Thanks, dariober, it working properly!

• 0 views • link
written by Jeason Rad • Jeason Rad

please check the green marks on the left to validate+close the question

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Hi all. Any way to modify this so it finds read with deletions > 20bp?

samtools view -h in.bam chr1:12340-200000 \
| awk '$1 ~ "^@" || $6 ! "20D"' \
| samtools view -b - > out.bam

Would give reads with exactly 20 bp, right? Is it possible to do > 20bp?

• 0 views • link
written by shimbalama • shimbalama

@shimbalama - If you really wanted to you could write an awk script that parses the cigar string and checks for a D operator > 20. But I think it would be easier to write a python script using the pysam package to parse the cigar string or see if Pierre's tools can do it (probably they can). If you get stuck, maybe open a new question for it.

• 0 views • link
written by dariober dariober

using samjs http://lindenb.github.io/jvarkit/SamJavascript.html

samtools -bu view your.bam chr1:12340-200000 |\
java -jar dist/samjs.jar -e 'function accept(r) { if(r.getReadUnmappedFlag()) return false; var cigar=r.getCigar();if(cigar==null) return false; for(var i=0;i< cigar.numCigarElements();++i) {if(cigar.getCigarElement(i).getOperator().isIndelOrSkippedRegion()) return false; } return true;} accept(record);'
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this