How to extract reads with no INDEL?

Question: How to extract reads with no INDEL?

3.8 years ago by

Jeason Rad • 30 wrote:

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?

sam samtools indel • 2.4k views

ADD COMMENT • link •

modified 3.8 years ago by dariober ♦ 11k • written 3.8 years ago by Jeason Rad • 30

3.8 years ago by

dariober ♦ 11k

WCIP | Glasgow | UK

dariober ♦ 11k wrote:

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

ADD COMMENT • link written 3.8 years ago by dariober ♦ 11k

Thanks, dariober, it working properly!

ADD REPLY • link written 3.8 years ago by Jeason Rad • 30

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

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

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?

ADD REPLY • link written 7 weeks ago by shimbalama • 10

@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.

ADD REPLY • link written 7 weeks ago by dariober ♦ 11k

ask as a new question please

ADD REPLY • link written 7 weeks ago by Pierre Lindenbaum ♦ 134k

3.8 years ago by

Pierre Lindenbaum ♦ 134k

France/Nantes/Institut du Thorax - INSERM UMR1087

Pierre Lindenbaum ♦ 134k wrote:

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);'

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

Thx for your help, pierre!

ADD REPLY • link written 3.8 years ago by Jeason Rad • 30

Please log in to add an answer.

Similar posts • Search »