De-duplicate to remove PCR artifacts in BS-seq data

[Jerry_Zhao]

Hi, All,

I have a question for de-duplicate for the bisulfite sequencing data to remove PCR atfifacts.

There are mainly two methods to do this.

1. The Lister method. This method consider mapped reads with identical 5' end (start site) as "clonal" reads, only keep the one with highest sequencing quality scores, and exclude others.

For example:
ATATATCGTAGTGGACCGTAACTGACGTTTTCAGC
-------------------
------------------
-----------------
----------------

2. The method used in deduplicate_bismark_alignment_output.pl of the Bismark package.
In the 163 line of this perl file, my $composite = join (":",$strand,$chr,$start,$end);
If I understand correctly, this mean that reads with identical start and end sites were considered as duplicates.

For example,
ATATATCGTAGTGGACCGTAACTGACGTTTTCAGC
----------------
----------------
----------------
----------------


However, both of these two methods do not consider the mapping status.
For example, can we consider the four reads bellow as duplicates?

ATATATCGTAGTGGACCGTAACTGACGTTTTCAGC
-----C----------
-----T----------
-----T----------
-----C----------


Thanks.
Jerry

[fkrueger]

Hi Jerry,

The best way of deduplicating is probably a combination of start position (and end-position for paired-end alignments), orientation and the actual sequence. It is also a matter of how much energy you want to spend on this and how many system resources you have. You would need to decide whether a sequence with e.g. a single sequencing error would count as a unique sequence or whether you might just classify it a duplicate regardless. A similar question arises in cases where a C position is observed as C in one read or T in another (like in your last example), would you treat this the same way as a mismatch at another position?

It is probably also a matter of feasibility, since keeping for example all chromosome, start and end postions, all read orientations as well as all sequences including quality scores in memory may require a LOT of RAM for lets say a 2x100bp HiSeq lane with > 200 M alignments.

For most WGSBS cases it is probably sufficient to just discard everything that looks like a PCR duplicate, and rely on getting enough coverage by authentic sequences spanning a certain position (like the one marked in bold in the following example) while starting at different positions, like so:

ATATAT[B]C[/B]GTAGTGGACCGTAACTGACGTTTTCAGC...
 TATAT[B]C[/B]GTAGTGGACCGTAACTGACGTTTTCAGC... 
  ATAT[B]C[/B]GTAGTGGACCGTAACTGACGTTTTCAGC...
   TAT[B]C[/B]GTAGTGGACCGTAACTGACGTTTTCAGC...
    AT[B]C[/B]GTAGTGGACCGTAACTGACGTTTTCAGC...
     T[B]C[/B]GTAGTGGACCGTAACTGACGTTTTCAGC...
      [B]C[/B]GTAGTGGACCGTAACTGACGTTTTCAGC...

This should theoretically still give you a maximum possible coverage of each position in the genome equivalent to the read length (or even higher for paired-end alignments), so you don't even have to start thinking about the "what is actually a PCR duplicate?" outlined above.

I realise that the solution in the deduplicate_bismark method is also not flawless, since if you quality trim your data prior to alignments you might end up with some reads that could be PCR duplicates but have a different read length because of quality trimming, which are thus not removed. If you want to be stricter you could change the line:

my $composite = join (":",$strand,$chr,$start,$end);

to

my $composite = join (":",$strand,$chr,$start);

I suppose if you are looking at sequencing small genomes where you expect a huge coverage or want to filter for really high read coverage you might indeed want to think about alternative ways of de-duplication.

[Jerry_Zhao]

Hi, Felix,

It is really helpful for me. I don't know how to define PCR duplicates before.
As you mentioned, I will think more deeply to imagine what the duplicates might be after the trimming and removing adapter sequences.

Actually, my initial idea is change
my $composite = join (":",$strand,$chr,$start,$end);

to

my $composite = join (":",$strand,$chr,$start,$end,$cigar);
and the $cigar will be the sixth column of the SAM file (CIGAR).
Nevertheless, I am realized that this is not appropriate.

It is really hard to define PCR duplicates.
I think I will try both of the two method, and to see whether there is big difference between them or not.


Many thanks and best regards,
Jerry