[Brian Bushnell]

Interesting, I guess a good cut-off there is 20000.

I'm not really sure. The graph is strange. Genomax estimated that the tiles appear to be roughly 33k by 37k, so dist=20000 would span a third of the tile. I'm guessing that maybe this was a highly amplified library with lots of PCR duplicates, and the reason the slope decreases and eventually levels is due to a high error rate when using "subs=0". I will investigate further. But thanks @Genomax for posting it!

Gotta love needing different settings per machine.

Yes, it makes deciding on default parameters very difficult, especially when the company never really publicizes any of this so you need to discover it empirically...

[GenoMax]

This sample had

Reads In:         1606392082
Clumps Formed:     356710875

Based on the name of the file it must have come from a S1 flowcell.

[dpryan]

Interesting, I guess a good cut-off there is 20000. Gotta love needing different settings per machine.

Edit: Out of curiosity, what sort of percentage would 100 million reads represent? I also wonder if the setting will depend on which of S1/S2/S3/S4 flow cells are being used.

[GenoMax]

I looked at one sample from NovaSeq data that Illumina made available recently. Here is a scan for dupedist= (spantiles=f dupesubs=0) values similar to one @Devon posted above.

Temporarily removing the image until @Brian reconfirms things.

[Brian Bushnell]

Keep in mind that 10-100 on the HiSeq 2500 is equivalent to 100-1000 on HiSeq 3000/4000/X, apparently Illumina multiplied everything by 10 (because hey, why not?).

Oh, thanks for the heads-up! Now Genomax's earlier observations make a lot more sense; on earlier platforms you see the asymptote well before 100.

Also, it looks like maybe spantiles should default to false since it is really only applicable to the NextSeq. The reason it shows an immediate increase at 1 is because cross-tile coordinate-based duplicate detection does not use Euclidean distances like the intra-tile duplicates. Instead, since only one coordinate is expected to be comparable (either X or Y), it uses the minimum of the X difference and Y difference. So, at dist=1, you are basically comparing a cluster to every other cluster in the same row or column. This does not have much affect on HiSeq 2500 or NextSeq (other than for tile-edge duplicates) but it looks like it has a huge impact on patterned flowcells because everything is already aligned in rows or columns.

[GenoMax]

Will try dupedist=2500. Thanks for the heads-up.

Not many had been systematically looking at duplicates (optical or otherwise) since Picard requires alignments for marking dups. We may find more duplicates than we expect/imagine once we start looking routinely.

[dpryan]

Keep in mind that 10-100 on the HiSeq 2500 is equivalent to 100-1000 on HiSeq 3000/4000/X, apparently Illumina multiplied everything by 10 (because hey, why not?).

Anyway, yes, it's essentially a linear increase starting at 20 (or 0 if spantiles is enabled).

Edit: Originally wrote 8 instead of 20. I had something else in mind apparently...

[GenoMax]

Interesting. I was only looking at distances of up to 100 (does a plot of 10-100 show consistent increase, it should). It does look like 2500 should be a reasonable setting.

[dpryan]

I'm attaching two plots from a sample with ~152 million reads, one truncated at a distance of 10000 and another going all the way to 50000. For what it's worth, I noticed that the Broad is using a distance of 2500 for patterned flow cell, which seems pretty reasonable. If one enables tile spanning, then you don't see saturation until ~20000, which seems a bit over the top.

The tile spanning results seem a bit over the top, though interestingly a distance of 1 is sufficient to find stuff with that enabled. I'll post the density of duplicates according to X/Y coordinate to ensure there's no NextSeq-like tile edge effect.

Update: Yup, no tile-edge effect, so not spanning tiles makes sense.

[Brian Bushnell]

Feature request: It'd be quite nice to be able to write marked duplicates to a different file or files. At the moment, I have to mark duplicates and write everything to a temporary file, which is then processed. Granted, one CAN use "out=stderr.fastq" and send that to a pipe, but then one needs to deal with all of the normal stuff that's written to stderr.

The impetus behind this is removing optical duplicates before delivery to the our labs but still writing them to a separate file or files in case they need them for some reason.

I will plan to add a new output stream for duplicate files as well, though I might not get to it this week.

[Brian Bushnell]

Additionally, is there any way to make clumpify itself respect the "threads=" setting? pigz seems to, but clumpify itself seems to use as many as it can get regardless of what I specify. This is in version 36.86.

Oh, hmmm... that will be very tricky. When running with one group, Clumpify should respect threads correctly. But when writing temp files (when happens whenever the reads won't all fit in memory), it uses at least one thread per temp file, and the default is a minimum of 11 temp files. Your best bet, unfortunately, would be to bind the process to a certain number of cores. You can also manually set the number of groups which indirectly affect the number of threads used.

Clumpify also uses multithreaded sorting, which uses all available cores, but normally that only happens for a small fraction of the runtime. However, I will add a flag to disable it.

[GenoMax]

Thanks, I'm running this now on a single sample, I'll post in image when I have a worthwhile sweep range.

Do the sweep with spanfiles=f and t. I was only interested in optical duplicates when I did mine.

[dpryan]

For now use the following workaround provided by @Brian.

clumpify.sh in=x.fq out=y.fq markduplicates [optical allduplicates subs=0]
filterbyname.sh in=y.fq out=dupes.fq names=duplicate substring include
filterbyname.sh in=y.fq out=unique.fq names=duplicate substring include=f

Thanks, in the interim I just wrote something in C that I can just call once to do this (it also strips "duplicate" from the read names).

This is a bit murky. I have done the sweeps with 4000 data I have access to. If I keep the spantiles=f then I don't see any optical dups until dupedist=20. Note: The edge duplicates problem seen with NextSeq (which has @Brian setting spantiles=t by default) is not present in HiSeq 4000/MiSeq (again based on data I have seen).

Thanks, I'm running this now on a single sample, I'll post in image when I have a worthwhile sweep range.

[dpryan]

Additionally, is there any way to make clumpify itself respect the "threads=" setting? pigz seems to, but clumpify itself seems to use as many as it can get regardless of what I specify. This is in version 36.86.

[GenoMax]

Feature request: It'd be quite nice to be able to write marked duplicates to a different file or files. At the moment, I have to mark duplicates and write everything to a temporary file, which is then processed. Granted, one CAN use "out=stderr.fastq" and send that to a pipe, but then one needs to deal with all of the normal stuff that's written to stderr.

That request has been in for some time I also wanted to see counts (with associated sequence) to see how acute of a problem the duplicates may be.

For now use the following workaround provided by @Brian.

clumpify.sh in=x.fq out=y.fq markduplicates [optical allduplicates subs=0]
filterbyname.sh in=y.fq out=dupes.fq names=duplicate substring include
filterbyname.sh in=y.fq out=unique.fq names=duplicate substring include=f

BTW, do you have any recommendations for the "dist" parameter on a HiSeq 4000? I was planning to just do a parameter sweep, but if that's already been done by someone else...

This is a bit murky. I have done the sweeps with 4000 data I have access to. If I keep the spantiles=f then I don't see any optical dups until dupedist=20. Note: The edge duplicates problem seen with NextSeq (which has @Brian setting spantiles=t by default) is not present in HiSeq 4000/MiSeq (again based on data I have seen).

I have not pulled out the reads using the method above to look at the co-ordinates/sequence as yet.

It may be good to see what you get.