This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Is it possible to remove reads from BAM which match the reference at a position

I'm interested in looking for modified bases and for the downstream analysis I'm doing I'd like to remove the reads which match the reference base at a certain position. Is there an easy way to do this?

I'm using a sorted bam file, and in my reference at position 600 the base is an A, is there an easy way to remove all reads from the bam file with an A at this position, and only leave behind the modified bases?

alignment sequencing samtools
• 2,797 views • link
written by yryan • yryan

1 answer

Using samjdk: http://lindenb.github.io/jvarkit/SamJdk.html

and the following snippet.

final boolean accept=true;
final String contig="RF02";
final byte base='A';
final int position=877;


if(record.getReadUnmappedFlag()) return accept;
if(!record.getContig().equals(contig)) return accept;
if(record.getAlignmentEnd() < position) return accept;
if(record.getAlignmentStart() > position) return accept;
final byte bases[] = record.getReadBases();
if(SAMRecord.NULL_SEQUENCE==bases || bases==null) return accept;
final Cigar cigar = record.getCigar();
if(cigar==null || cigar.isEmpty()) return accept;
final int readpos1 = record.getReadPositionAtReferencePosition(position,false);
if(readpos1 < 1) return accept;
final int readpos0 = readpos1-1;
if(readpos0 < 0 || readpos0 >= bases.length) return accept;
if(bases[readpos0]!=base) return accept;
return !accept;

usage:

java -jar dist/samjdk.jar -f snippet.code in.bam
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

That worked perfectly, thanks so much for that!

• 0 views • link
written by yryan • yryan

Log in to answer this question.

More posts like this