Extracting reads with a particular mismatch from bam file

[DerSeb]

Hi @ll,

I want to do the following and have not found any suggestions on the web yet:

I would like to extract reads from a bam file that contain a certain mismatch to the reference.

This means irrespective of the remaining sequence, if a read contains for example a "G" instead of the reference allele "A" at any position I would like to write it into a separate bam file.

Do you think that is possible?

I do not want to call SNPs or get the number of non-ref bases in a region. I really want to extract the entire read with at least one specific mismatch (e.g. only A->G), no matter if there is a SNP at that position or not.

Thanks for you attention,
Seb

[dpryan]

I doubt there's a premade tool for that, but you could do that with pysam by first reading the genome into memory, parsing the MD string to find mismatches, seeing if it's of the type you want, and then writing the read to a new file.