Question: PCR duplicate reads - samtools
gravatar for dina91n
6 weeks ago by
dina91n0 wrote:

Dear All,

I am trying to find the PCR duplicate reads from the bam/sam. I have used picard to mark duplicates in the bam/sam file and then samtools flagstats. I obtained the following output:

196499952 + 0 in total (QC-passed reads + QC-failed reads)
125920476 + 0 secondary
0 + 0 supplementary
11613716 + 0 duplicates
195906701 + 0 mapped (99.70% : N/A)
70579476 + 0 paired in sequencing
35289738 + 0 read1
35289738 + 0 read2
68711404 + 0 properly paired (97.35% : N/A)
69739378 + 0 with itself and mate mapped
246847 + 0 singletons (0.35% : N/A)
127264 + 0 with mate mapped to a different chr
31498 + 0 with mate mapped to a different chr (mapQ>=5)

Is the number of duplicate reads calculated with respect to the number of the QC-passed reads?

Thank you in advance

ADD COMMENTlink modified 6 weeks ago by swbarnes27.5k • written 6 weeks ago by dina91n0

Yes, this is obvious as you do not have any QC-failed reads.

ADD REPLYlink written 6 weeks ago by ATpoint31k

To make it clear, the "+ 0" on every line refer to the number of Qc-failed reads in each category. And the number to the left of the "+" always refer to Qc-passed reads.

ADD REPLYlink written 6 weeks ago by Carlo Yague4.9k
gravatar for swbarnes2
6 weeks ago by
United States
swbarnes27.5k wrote:

Your "196499952 + 0 in total" is purely a count of the number of lines in the bam. You also have 125920476 lines of your bam marked as secondary alignments; that is, your aligner assigned the same read to multiple possible positions. So you have the same reads in there over and over again, and they might be marked as duplicates over and over again.

You might want to filter secondary alignments away to get a better count.

ADD COMMENTlink written 6 weeks ago by swbarnes27.5k

Thank you for your answer. Therefore, should I filter the secondary alignments to get the real number of duplicate reads? How can I do this? Is there a faster way to get the number of duplicate reads? Thank you in Advance!

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by dina91n0

Why doesn't samtools view work for you?

ADD REPLYlink written 6 weeks ago by swbarnes27.5k

Thank you for your hint! I have two questions: 1) I used the following command to filter out secondary alignments and to get the number of reads aligned (number of lines in the file) :

samtools view -h -F 256 file.bam > primaryAlign.bam
grep -v '^@'  primary_align.bam | wc -l

the number of lines is 70579476; then I extracted the duplicate reads and get the number of the rows in the file:

samtools view -h -f 1024  primaryAlign.bam > dupl.bam
 grep -v '^@'  dupl.bam  | wc -l

the number of rows is 11613716. So does that mean that all of my duplicate read counts are in the first alignment?

2) when I aligned my reads with bowtie2, I got the following output:

35289738 reads; of these:
  35289738 (100.00%) were paired; of these:
    934036 (2.65%) aligned concordantly 0 times
    25038058 (70.95%) aligned concordantly exactly 1 time
    9317644 (26.40%) aligned concordantly >1 times
    934036 pairs aligned concordantly 0 times; of these:
      330737 (35.41%) aligned discordantly 1 time
    603299 pairs aligned 0 times concordantly or discordantly; of these:
      1206598 mates make up the pairs; of these:
        593251 (49.17%) aligned 0 times
        147482 (12.22%) aligned exactly 1 time
        465865 (38.61%) aligned >1 times
99.16% overall alignment rate

How is it possible that the number of reads aligned is 70579476, if I got 99.16% of the alignment rate?

Thank you in advance!

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


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