I have an indexed BAM file containing long-read sequencing data and I'd like to split the reads contained within into those with a known deletion and those without the deletion (I have the deletion coordinates available to me) when mapped against the hg38 genome. Ideally I'd like to end up with two BAM files (one with reads containing the deletion and one without) for downstream analysis. I've searched a bit online and haven't really found any tools that do exactly that, but this seems like a common enough task that there has to be something.

Thanks in advance for any help/guidance!


In addition to @gringer's great answer, I decided to write a small python script myself that solves my issue using pysam, the fruits of which can be found here. If anyone has any feedback or suggestions on it I'll happily take them, as I'm still very new to this. It's undergone very little testing, so user be warned.


I was asked to write some code sometime near the end of last year to do this, so that viral sequences could be split based on variants at a particular location. The code I wrote doesn't actually do the splitting - it creates SAM read groups based on those variants - but that splitting can be done after processing with this script using samtools split.

As with most of my scripts, it hasn't had much testing / use:



samtools view -h mapped_reads.bam | ./samVarSplitter.pl [-ref <ref>] [-pos <int>] [options]


samtools view -h mapped_reads.bam | ./samVarSplitter.pl -ref Wuhan-Hu-1 -pos 11083

If the input file contains a lot of reads that don't include the variant location, I'd recommend pre-filtering reads to only include those that cover that location, because expanding the CIGAR string and storing a record of the positions can take some time.

Note: currently the position argument only supports a single location. I could extend this to multiple positions (i.e. a position range) if that were a desired feature.

  • $\begingroup$ Thanks so much, this is pretty much what I was looking for! If I understand it correctly, the script only works with a single position though correct? The deletion I'm looking at is 9bp in size, but I suppose I can just use one of the positions as a stand-in for the entire deletion. $\endgroup$
    – jazzbo
    Aug 31 '20 at 18:54

The CIGAR string should tell you if a read has a deletion in it. So filter the bam files based on that.

  • 1
    $\begingroup$ The OP is looking for a known deletion, so if you look to identify reads in the region with 'D' in the cigar, you might include other deletions $\endgroup$
    – StupidWolf
    Aug 27 '20 at 8:26

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.