hi, all I wrote a C program to stat the depth of each position. The input is a BAM file. Before I stat it, I used Picard Markduplicates to flag the duplicate reads instead of remove them directly(Don't ask me why, my boss decided it). But I wanna stat the rmdup depth, like what the samtools depth did. I planned using (flag &1024) to screen the dup reads. But all the reads screened. Yes, I know it means every read at least have one optical dup or PCR dup, and they are duplicate reads for each other. but I just wanna stat each dup group reads once. I don't wanna check the quality value of each read, and keep the best one. It's waste a lot of CPU times and mem. Anyone can tell me a easier way? THX~
0 answers
No answers yet.
Log in to answer this question.
-
Divergent number of replicates between FastQC and samtools
written by user31888
-
Where to find PCR duplicate reads in bam file?
written by bioinforesearchquestions
-
Remove Duplicates Workflow
written by Anjan
-
Picard Mark Duplicate Reads In Galaxy
written by Tonyzeng •
-
Marking Duplicates With Molecular Tag (Umi)
written by brentp
-
Does Samtools Tview Hide Pcr And Optical Duplicates?
written by Matthieu Miossec
-
Denovogear Support Marked Dup And Sorted Bam Files?
written by kulkarni.aishwarya87 •
-
Bam File: Which Column Change When We Flag The Duplicate Reads
written by star.bioinformatics •
-
Samtools Dedup Documentation
written by Zev.Kronenberg
-
Samtools Depth And Gaps
written by Juliofdiaz
I'm not sure how to interpret "..all the reads screened. Yes, I know it means every read at least have one optical dup or PCR dup, and they are duplicate reads for each other". Perhaps you could clarify what you mean by that. I recall that MarkDuplicates leaves one read (or pair) from a group unmarked (i.e., bit 1024 isn't set, I think this is done for the read or pair with the highest MAPQ), which would seem to facilitate what you're trying to do. Yes, you still have to iterate through all of the reads to check their flag, but that's relatively quick and there's really no by-passing it.
Thanks for you comment, and sorry about my English.. I use Picard mark the duplicate reads. but all the reads flag be marked (bit 1024 is set), not just the highest MAPQ one. have you meet this problem? I think should I check the MAPQ of each read and find the best one. or have an easier way. I just add a picture of my result for you information.
As dpryan79 mentioned above that Picard should mark all but one read (with highest MAPQ) as duplicates. Check it again carefully. Also, paste a few line from the bam file representing all the reads with identical start position where all of them have been marked duplciates.