Extract reads with defined nucleotide at defined aligned position

[charlescoldroom]

Dear all,

I would like to extract from a BAM (or SAM) file the reads that contain a specified nucleotide at a specified location. E.g. I would like to extract all reads that have the reference allele (T) at chr1:222222.

All the tools that I could find around simply subset a BAM based on an input position, but I cannot find anything that is able to subset based on which base is at that position.

Any ideas on how this could be accomplished?

Thank a lot!

[dpryan]

You'll probably have to put something together to do this for you. If you're familiar with python, you could do this with pysam (first subset the BAM/SAM file with samtools to just get reads covering that position and then read the results into a python program, iterating over the reads and just outputing those with a T at the desired position). You might be able to do the subsetting with pysam to, I've never tried.

If you're more comfortable with C, there's the samtools C API, which can do the same things. If you're more of a java person then there's the picard API. Presumably there's something for perl too.

[charlescoldroom]

Hi dpryan,

thanks a lot for the answer! pysam pileup does exactly what I needed, we are now using it to look up one position in a bam file and then output a list of all reads names and for each one the type of nucleotide at the desired position.

Thanks a lot for your help!