Question: Tool to unmark duplicates
3
gravatar for Mark Ebbert
4.1 years ago by
Mark Ebbert70
United States
Mark Ebbert70 wrote:

Hi,

I received some bams where duplicates were marked using Picard. I'd like to unmark all duplicate reads for a study we're doing. Does anyone know of an existing tool? I haven't found anything in google land.

If there isn't an existing tool, could you provide some guidance how to unmark them myself? I've read the SAM spec (http://samtools.github.io/hts-specs/SAMv1.pdf) and understand I need to change the bit flag, but I'm not sure the best way to approach it. 

Thanks!

Mark

ADD COMMENTlink modified 4.1 years ago by Pierre Lindenbaum111k • written 4.1 years ago by Mark Ebbert70

Hello me.mark!

It appears that your post has been cross-posted to another site: SEQanswers

This is typically not recommended as it runs the risk of annoying people in both communities.

ADD REPLYlink modified 4.1 years ago • written 4.1 years ago by Devon Ryan82k

Thank you for the reminder. My colleague and I were both looking for a solution and it looks like he also posted a question about it. Please excuse the duplicate post.

Edit: Thank you everyone for your helpful responses!

ADD REPLYlink modified 4.1 years ago • written 4.1 years ago by Mark Ebbert70

That explains the very different usernames, no worries!

ADD REPLYlink modified 4.1 years ago • written 4.1 years ago by Devon Ryan82k

This post explains how to do this using Picard: Remove flags of MarkDuplicates (picard)

ADD REPLYlink written 3.8 years ago by Malachi Griffith16k
4
gravatar for Pierre Lindenbaum
4.1 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum111k wrote:

as said Devon, writing this tool is very easy ;-)

I wrote it and put it here: https://github.com/lindenb/jvarkit/wiki/Biostar106668

Example:

$ curl -s "ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/phase1/technical/ncbi_varpipe_data/alignment/NA12751/NA12751.chrom20.ILLUMINA.mosaik.CEU.low_coverage.20101123.bam" |\
samtools view -h - |\
head -n 10000 |\
java -jar dist/biostar106668.jar -b > out.bam

 

 

ADD COMMENTlink written 4.1 years ago by Pierre Lindenbaum111k
1

I guess my "has written" should have been "will shortly write" :)

ADD REPLYlink written 4.1 years ago by Devon Ryan82k
4
gravatar for Ashutosh Pandey
4.1 years ago by
Philadelphia
Ashutosh Pandey11k wrote:

Here is something you can try.  I haven't tested it but I guess it should work. Let me know if you get any error while reconverting to bam.     samtools view -h input.bam | awk '{printf "%s\t", $1; if(and($2,0x400)){t=$2-1024}else{t=$2}; printf "%s\t" , t; for (i=3; i<NF; i++){printf "%s\t", $i} ; printf "%s\n",$NF}'| samtools view -Sb - > input.unmarked.bam

ADD COMMENTlink modified 4.1 years ago • written 4.1 years ago by Ashutosh Pandey11k

It should be noted that the awk and() function doesn't always exist (e.g., on OS X).

ADD REPLYlink written 4.1 years ago by Devon Ryan82k

It should also be noted that "you can't do binf on a mac"

ADD REPLYlink written 3.5 years ago by Aaron Statham1.1k

I wouldn't go that far. OS X is essentially a BSD variant with a really nice user interface. The and() function is just a difference between gawk and mawk and one can always install the other one.

ADD REPLYlink written 3.4 years ago by Devon Ryan82k

you could shorten this by just setting $2-=1024 in your if and then printing $0.

ADD REPLYlink written 4.1 years ago by brentp22k
1
gravatar for Devon Ryan
4.1 years ago by
Devon Ryan82k
Freiburg, Germany
Devon Ryan82k wrote:

I'm not familiar with any existing tool (I wouldn't be surprised if Pierre Lindenbaum has written such a tool), but this should be simple enough to code. Depending on the version of awk you have on you system, you may simply be able to use it. In some implementations of awk (not all of them), the command if(and($2,1024)){$2-=1024} will unmark a duplicate. You would simply samtools view -h foo.bam | awk ...stuff... | samtools view -Sb - > foo.unmarked.bam.

The other simple way to do this is with python using pysam. Pysam allows reading from and writing to SAM/BAM files from within python, so you could perform a similar operation there.

ADD COMMENTlink written 4.1 years ago by Devon Ryan82k
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: 1019 users visited in the last hour