1
$\begingroup$

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!

$\endgroup$
2
$\begingroup$

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...
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.