Question: Is it possible to remove reads from BAM which match the reference at a position
0
gravatar for yryan
5 weeks ago by
yryan0
yryan0 wrote:

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?

sequencing samtools alignment • 192 views
ADD COMMENTlink modified 5 weeks ago • written 5 weeks ago by yryan0
2
gravatar for Pierre Lindenbaum
5 weeks ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum126k wrote:

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 COMMENTlink written 5 weeks ago by Pierre Lindenbaum126k

That worked perfectly, thanks so much for that!

ADD REPLYlink written 5 weeks ago by yryan0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 796 users visited in the last hour