I am handling a few NGS datasets. While using samtools rmdup, sometimes it removes 10-15% of reads and sometimes 70-80% of reads. It is quite obvious that the second situation might have arisen due to wrong sequencing chemistry. But really I want to know, what is the basis of removing PCR duplicates? I have found some explanations on SEQanswers but they were not what I needed. Looking for your suggestions, and if possible also links about SAMTools rmdup.
I will just try to put the conceptual idea behind it as far as I have understood.
In a paired-end data, if there are few pair of reads aligned at same position, samtools mark them as PCR duplicates and retain only one pair (mostly the first one).
Its simple as that in paired end sequencing, the only possibility to get many reads (as a pair which represents the whole fragment) aligned at same position is highly likely due to overamplicification of the same fragment in PCR (possibly due to more number of PCR cycles or by starting with fewer/less diverse starting material of DNA). Beacuse, DNA getting sonicated at same start and end point in different DNA is very least probable. So, rmdup marks these ones as duplicates to avoid reads projecting less confident variations as a higher one based on number of reads coming from the same fragment which got amplified in PCR.
Samtools is not considering the read sequence aligned at those position, it just looks at the position, thats it. If there are reads aligned at same position but their sequences are different, still it will mark that as duplicates. The facts to keep in mind about this issue is that the reasons for having different sequence reads aligning to same position.
1. Some reads could have one or two sequencing errors.
2. False alignment which is less probable.
In both the cases, its good to be removed. So no need to worry about that.
I am not quite sure, but I think the same applies to single-end. But it wont be accurate enough because the tool only knows start point of the fragment but not the end where you could loose many reads.
As a prelimiary step, you can check the raw data in FastQC which roughly tells about the duplication level in reads. Though the tool counts for first 200,000 sequences only, we will get to know its level roughly.
I will just try to put the conceptual idea behind it as far as I have understood.
In a paired-end data, if there are few pair of reads aligned at same position, samtools mark them as PCR duplicates and retain only one pair (mostly the first one).
Its simple as that in paired end sequencing, the only possibility to get many reads (as a pair which represents the whole fragment) aligned at same position is highly likely due to overamplicification of the same fragment in PCR (possibly due to more number of PCR cycles or by starting with fewer/less diverse starting material of DNA). Beacuse, DNA getting sonicated at same start and end point in different DNA is very least probable. So, rmdup marks these ones as duplicates to avoid reads projecting less confident variations as a higher one based on number of reads coming from the same fragment which got amplified in PCR.
Samtools is not considering the read sequence aligned at those position, it just looks at the position, thats it. If there are reads aligned at same position but their sequences are different, still it will mark that as duplicates. The facts to keep in mind about this issue is that the reasons for having different sequence reads aligning to same position.
1. Some reads could have one or two sequencing errors.
2. False alignment which is less probable.
In both the cases, its good to be removed. So no need to worry about that.
I am not quite sure, but I think the same applies to single-end. But it wont be accurate enough because the tool only knows start point of the fragment but not the end where you could loose many reads.
As a prelimiary step, you can check the raw data in FastQC which roughly tells about the duplication level in reads. Though the tool counts for first 200,000 sequences only, we will get to know its level roughly.
The rationale is that given the low coverage of the reads w.r.t. the genome it is extremely improbable for two or more reads to fall in the same exact position by chance (or, in other words, to be generated by the same position of the genome). Thus, when this happens it is considered an amplification artifact. Obviously this reasoning is meaningful if you are working with big genomes and low coverage, e.g. a ChIP-Seq of 20x10^6 reads on the human genome. If you are working with small genomes, e.g. the same ChIP-seq on yeast than the situation is much different.
I am not a specialist but you can have a look at the Broad Institute GATK workshop: http://www.broadinstitute.org/videos/broade-mapping-and-duplicate-marking
It's a presentation about aligning the reads and removing the duplicates. They are using Picard, but I think the basic idea is the same.
One point to add is that apparently samtools does not handle inter-chrosomosome paired reads, while picard does. I have been using picard quite a bit and it seems to work well for me in the case of WGS. Normally, you may see up to 10% of PCR duplicates for Illumina paired-end data. But I have also seen some really serious problem with PCR duplicates in the case of paired-mate data. Having the PCR duplicates removed before detecting the variants, especially those that utilize some sort of frequency (or number of supporting reads) as cutoff, is very important.
One thing in addition to the answers here: samtools rmdup removes all but one read-pairs. The remaining read-pair is the one having the highest mapping quality. This implicates to have the lowest number of erros/variations possible.
Not sure how it exactly works on single-end data, presumably on start and stop coordinates of the read...
Thanks a lot to everyone. Yes, I am dealing with paired end data. So all your inputs really helped me to understand the issue. Mohamed Ashick mentioned about sequencing errors. How to rule out that those are sequencing errors but not variants even if those paired reads aligned to the same start and end position w.r.t. genome. Does phred score of bases has any influence to support the comment. Also I want to mention that, he dataset from which 70-80% reads are removed is the first run of Illumina HiSeq1000. Can it be the reason for the same.
Thanks in advance.
To answer your first question, its highly unlikely that paired reads aligning to same position and having a variation other than the fact that it could be sequencing error. That is why samtools doesnt consider sequences other than positions. If there are chunk of reads aligned exactly at same position bt having difference in bases, its supposed to be Sequencing error or mis-alignment. Even if there is a slight possibility, there are Mapping quality, SNP quality and read depth thresholds which will take care to pick actual variations.
Coming to your second question, HiSeq has high thoughput compared to GAIIx. So, there could be extra PCR cycles added to meet the throughput. I am not sure about this, but possible.
- Badges
- Science topic
- Similar topics
- Bioinformatics
- Bioinformatics and Computational Biology