How to extract reads with INDELs > a given size?

Question

I'm trying to modify this https://www.biostars.org/p/253774/

To get reads with deletions > 20bp

I think this gives reads with exactly 20bp dels:

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

is it possibe to do >?

Thanks!

Answer 1

There's no simple way to do that with awk, since you would need to parse the CIGAR string and iterate over the operations. In python using pysam you would iterate over the reads and filter as follows:

keep = False
for OP, LEN in b.cigartuples:
    if OP == 2 and LEN > 20:
        keep = True
        break

if keep:
    ...do stuff...