Is it possible to remove reads from BAM which match the reference at a position
1
0
Entering edit mode
2.5 years ago
yryan ▴ 10

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 • 1.4k views
ADD COMMENT
2
Entering edit mode
2.5 years ago

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
ADD COMMENT
0
Entering edit mode

That worked perfectly, thanks so much for that!

ADD REPLY

Login before adding your answer.

Traffic: 836 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6