Sort reads in BAM file based on presence of specific deletion?

Question

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!

UPDATE

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.

Answer 1

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:

https://gitlab.com/gringer/bioinfscripts/-/blob/master/samVarSplitter.pl

usage:

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

example:

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.

Answer 2

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