SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Can Cuffdiff treat paired-end and single-end reads at the same time? zun RNA Sequencing 3 06-12-2012 05:37 PM
Can paired-end mapping produce more reads than single-end ? warrenemmett Bioinformatics 13 03-20-2012 11:10 PM
Paired-end Bam from single-end aligned sam ramouz87 Bioinformatics 4 08-17-2011 12:55 PM
RNA-seq: Replicates, single-end, paired-end story pasta Bioinformatics 2 07-04-2011 11:51 PM
Does Cufflinks support single-end and paired end data together ? ersenkavak Bioinformatics 1 10-22-2010 07:26 AM

Reply
 
Thread Tools
Old 02-16-2011, 07:42 AM   #1
m_elena_bioinfo
Member
 
Location: Ospedali Riuniti di Bergamo, ITALY

Join Date: Oct 2009
Posts: 99
Default duplicate PCR Single and Paired End

Hi NGS user,
I'm developing a script to remove PCR duplicates. I observed that samtools rmdup function doesn't always work well.
I'm using bioperl object and Bio:B::Sam utilities.
I have to distinguish if a bam is paired-end or single-read. Is there anyway to do this by Bio:B::Sam?

To remove duplicate,
if the NGS experiment is a single read, I considered duplicates read with the same start/end and with the same sequence.
If the NGS experiment is a paired end, a read is a PCR duplicates if has the same start of the first mate and the end of the second mate. Is it right?

Thanx a lot!
ME
m_elena_bioinfo is offline   Reply With Quote
Old 02-21-2011, 10:08 AM   #2
Thomas Doktor
Senior Member
 
Location: University of Southern Denmark (SDU), Denmark

Join Date: Apr 2009
Posts: 105
Default

I would say that PCR duplicates in SE sequencing can appear as having different end coordinates and different sequences because the PCR process itself can induce errors and the sequencing can introduce errors as well. That is why it's really difficult to remove actual PCR duplicates with single end reads because you really don't know if it's caused by PCR duplication, real enrichment (e.g. in an RNA-seq or CHIP-seq) or simply saturation of the sequencing depth. Optical duplicates in an illumina run can be identified by the actual cluster proximity and they will have identical sequence.

If you want to be really thorough and remove all duplicates (PCR or not), and you trust the quality of the alignment, you could define PCR duplicates in a single end experiment as reads with identical start point (5' of the read, so in the SAM/BAM format you would have to use the CIGAR string to determine the start position of the reverse strand reads) and then ignore any sequence information or indels in the read. That would remove all PCR duplicates, but also duplicate reads due to enrichment etc. It really depends on what type of analysis you want to make and the quality of the data.

PCR duplicates in PE experiments are much easier to identify since they would have equal total fragment length and fragment location, like you mention.

Last edited by Thomas Doktor; 02-21-2011 at 10:11 AM.
Thomas Doktor is offline   Reply With Quote
Old 02-21-2011, 10:14 AM   #3
m_elena_bioinfo
Member
 
Location: Ospedali Riuniti di Bergamo, ITALY

Join Date: Oct 2009
Posts: 99
Default

Thanx a lot Thomas for this clear answer!
m_elena_bioinfo is offline   Reply With Quote
Old 02-21-2011, 11:36 PM   #4
mapper
Member
 
Location: India

Join Date: Nov 2010
Posts: 17
Default

you can use picard tools to remove PCR duplicates...
mapper is offline   Reply With Quote
Old 02-22-2011, 03:49 AM   #5
Thomas Doktor
Senior Member
 
Location: University of Southern Denmark (SDU), Denmark

Join Date: Apr 2009
Posts: 105
Default

What if you wanted to collapse all PCR duplicates into a single read? Doesn't picard tools simply flag them as PCR duplicates and if you remove them you will lose that read position completely?
Or does picard tools only flag subsequent duplicates and not the first one it encountered?
Thomas Doktor is offline   Reply With Quote
Old 02-22-2011, 10:24 PM   #6
mapper
Member
 
Location: India

Join Date: Nov 2010
Posts: 17
Default

Quote:
Originally Posted by Thomas Doktor View Post
What if you wanted to collapse all PCR duplicates into a single read? Doesn't picard tools simply flag them as PCR duplicates and if you remove them you will lose that read position completely?
Or does picard tools only flag subsequent duplicates and not the first one it encountered?
I am not sure what "collapse all PCR reads mean"...Can you explain in more detail...

Picard mark reads as duplicates and remove it only when REMOVE_DUPLICATES is set to 'YES'...If you keep this option 'NO' it will just make duplicate and keep reads in file...Most of the downstream analysis tools considers uniquely mapped reads so one read per start position still gives all information....
mapper is offline   Reply With Quote
Old 02-23-2011, 03:57 AM   #7
Thomas Doktor
Senior Member
 
Location: University of Southern Denmark (SDU), Denmark

Join Date: Apr 2009
Posts: 105
Default

By collapsing I mean removing all but a single of the duplicate reads. My question is if picard tools marks all the duplicate reads as duplicates or keeps one read as the "original"?
Thomas Doktor is offline   Reply With Quote
Old 02-23-2011, 07:20 AM   #8
mapper
Member
 
Location: India

Join Date: Nov 2010
Posts: 17
Default

I believe it keeps one read per start position and discards other reads. I have seen these kind of alignments IGV...
mapper is offline   Reply With Quote
Old 02-23-2011, 09:57 AM   #9
Thomas Doktor
Senior Member
 
Location: University of Southern Denmark (SDU), Denmark

Join Date: Apr 2009
Posts: 105
Default

Thanks for clearing that up.
Thomas Doktor is offline   Reply With Quote
Old 02-23-2011, 12:11 PM   #10
Jon_Keats
Senior Member
 
Location: Phoenix, AZ

Join Date: Mar 2010
Posts: 279
Default

picard keeps the duplicate with best quality values and removes/marks the rest
Jon_Keats is offline   Reply With Quote
Old 03-02-2011, 01:14 PM   #11
Robin
Member
 
Location: US

Join Date: Nov 2009
Posts: 10
Default

I used the MarkDuplicate of Picard tools, and I got the very strange output. There are 0 count for UNPAIRED_READS_EXAMINED and UNPAIRED_READ_DUPLICATES. I found it is very strange.
Here is the detail output:
LIBRARY Unknown Library
UNPAIRED_READS_EXAMINED 0
READ_PAIRS_EXAMINED 25574519
UNMAPPED_READS 35323240
UNPAIRED_READ_DUPLICATES 0
READ_PAIR_DUPLICATES 5614734
READ_PAIR_OPTICAL_DUPLICATES 0
PERCENT_DUPLICATION 0.219544
ESTIMATED_LIBRARY_SIZE 49364927

Can someone explained it to me? Why there is no unpaired mapped reads found in the SAM file?
Thanks
R
Robin is offline   Reply With Quote
Old 03-03-2011, 04:04 AM   #12
Thomas Doktor
Senior Member
 
Location: University of Southern Denmark (SDU), Denmark

Join Date: Apr 2009
Posts: 105
Default

My guess is that it's not unpaired mapped reads, but simply unpaired reads. Both reads are present in the SAM file, even if only one of them is actually mapped. Picard is just telling you that it found all read pairs in the SAM file, not that all read pairs were mapped (since you have a lot of unmapped reads). If you remove all the unmapped reads, my guess is that Picard would report finding some unpaired reads.
Thomas Doktor is offline   Reply With Quote
Old 03-03-2011, 05:48 AM   #13
Robin
Member
 
Location: US

Join Date: Nov 2009
Posts: 10
Default

Please check the MarkDuplidate's column definition.

LIBRARY: The library on which the duplicate marking was performed.
UNPAIRED_READS_EXAMINED: The number of mapped reads examined which did not have a mapped mate pair.
READ_PAIRS_EXAMINED: The number of mapped read pairs examined.
UNMAPPED_READS: The total number of unmapped reads examined.
UNPAIRED_READ_DUPLICATES: The number of fragments that were marked as duplicates.
READ_PAIR_DUPLICATES: The number of read pairs that were marked as duplicates.
READ_PAIR_OPTICAL_DUPLICATES: The number of read pairs duplicates that were caused by optical duplication. Value is always < READ_PAIR_DUPLICATES, which counts all duplicates regardless of source.
PERCENT_DUPLICATION: The percentage of mapped sequence that is marked as duplicate.
ESTIMATED_LIBRARY_SIZE: The estimated number of unique molecules in the library based on PE duplication.


I don't believed that UNPAIRED_READS_EXAMINED = 0 count in my SAM file. From the Picard reports that I only have the unmapped reads, the mapped read pairs, and the duplicated mapped read pairs.
Robin is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 02:02 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO