interpretation of pcr duplication

[saturatedfunk]

We have high rates (50%) of duplication flags identified by picard using proton's rnaseq protocol. So I conducted an in-silico experiment below.

The data comes from 2 almost-technical replicates of human breast cancer. They differ in starting RNA input, 50 vs 25 ng of RNA.

0% optical reads were detected, as I wanted, but:

Relevant flags:
1024 1040
25ng 8615686 6934223
50ng 11481086 8776399

These reads represent roughly 50% of the total, and that has been consistent regardless of the starting input material or specimen.

Since the mean read length is short and these are single end reads, I thought the duplication rate might be caused by short, highly expressed transcripts.

Therefor, I merged the bam files from the replicates, set RG the same with AddOrReplaceReadGroups, then markedDuplicates on the merged bam

Results:
1024 1040
5025 observed 20508673 16097011

and expected results:
1024 1040
5025expectedminimum 20096772 15710622
5025 observed 20508673 16097011

where I expected the minimum number of duplicates to be the sum of the duplicates from the individual runs.

Since the observed duplication rate in the merged sample is only slightly higher, I conclude that the majority of original reads marked as duplicate really are pcr duplicates. And that the 'false pcr duplicates' rate is only about 3%, given this library preperation.

Is this interpretation correct?

[jwfoley]

The tables are hard to read and this logic is hard to follow:

Since the observed duplication rate in the merged sample is only slightly higher, I conclude that the majority of original reads marked as duplicate really are pcr duplicates. And that the 'false pcr duplicates' rate is only about 3%, given this library preperation.

Can you elaborate?

[dpryan]

I'm not sure that Picard's markDuplicates command respects read groups and, if not, I wouldn't conclude anything from your test. Having said that, high alleged duplication rates are expected in RNAseq due to highly expressed species (this is also why you don't bother marking duplicates in RNAseq unless you're looking for edit sites).

[saturatedfunk]

Sorry for the tables, they looked ok before I posted
I will attempt to simplify the question and describe my logic.

Replicate 1 had 15.5 million duplicates, which represented roughly 50% of total reads
Replicate 2 had 20.3 million duplicates, which represented roughly 50% of total reads.

Assuming,
1. These are truly PCR duplicates
2. There is 0% pcr duplication between the two runs ( since they were pcr'd seperately)

The assumption in #1 is a bit naive and has been discussed exhaustively, I do not intend to rehash it here. The assumption in #2 is true by definition.

I expect to find (at least) 35.8 million duplicates in the merged file.
I found 36.6 million duplicates in the merged file.

Therefor, (1- (35.8/36.6))=2% of the reads from replicate 1 share a read start site with replicate 2, but are not truly pcr duplicates. I refer to these 2% as "false-pcr-duplicate"

My question is, can I assume that the false pcr duplication rate in the individual runs is only 2%?