[Brian Bushnell]

OK, you can kill it - it won't finish. This is great feedback, by the way - I really appreciate it.

1) What version of BBMap are you using?
2) "unpair=f" is the default and that flag should always be ignored except when you are doing error-correction; it is exclusively for error-correcting paired reads. I don't see why that would have anything to do with this problem, but I'm not really sure.
3) This dataset is pretty small. Would it be possible for you to send it to me so I can try to replicate the issue?
4) I will note that your second version was run differently than the first one - it was using only 16GB RAM. The behavior of Clumpify is different when running with enough memory to hold all reads, and with not enough memory to hold all reads (in which case it needs to write temp files). It works in both cases, but when diagnosing errors, it's easiest to run with the same -Xmx parameter in all cases.

Thanks,
Brian

[chiayi]

Oh... in rare cases, "reorder" can cause it to run very slowly, stuck at 100% CPU utilization. That might be the issue here... try removing that flag. During normal execution, pigz is also using CPU-time and java is usually substantially higher than 100%. What did the screen output look like when top was showing this?

It was at the 'Deduping' step. I removed reorder tag and it still stuck at the same step:

openjdk version "1.8.0_121"
OpenJDK Runtime Environment (build 1.8.0_121-b13)
OpenJDK 64-Bit Server VM (build 25.121-b13, mixed mode)
java -ea -Xmx48g -Xms48g -cp ~/package/bbmap/current/ clump.Clumpify 
-Xmx48g in=in.fastq.gz out=out.fq.bz2 dedupe=t addcount=t dupedist=40 
optical=t
Executing clump.Clumpify [-Xmx48g, in=in.fastq.gz,
out=out.clumped.fq.bz2, dedupe=t, addcount=t, dupedist=40, optical=t]

Clumpify version 37.17
Read Estimate:          30447286
Memory Estimate:        13555 MB
Memory Available:       38586 MB
Set groups to 1
Executing clump.KmerSort [in1=in.fas in2=null, out1=out.clumped.fq.bz2, 
out2=null, groups=1, ecco=false, rename=false, shortname=f, 
[COLOR="Red"]unpair=false[/COLOR], repair=false, namesort=false, 
ow=true, -Xmx48g, dedupe=t, addcount=t, dupedist=40, optical=t]

Making comparator.
Made a comparator with k=31, seed=1, border=1, hashes=4
Starting cris 0.
Fetching reads.
Making fetch threads.
Starting threads.
Waiting for threads.
Fetch time:     95.473 seconds.
Closing input stream.
Combining thread output.
Combine time:   0.133 seconds.
Sorting.
Sort time:      36.426 seconds.
Making clumps.
Clump time:     6.788 seconds.
Deduping.

At the time the top looked like this:

PR  NI    VIRT    RES    SHR S  %CPU %MEM     TIME+ COMMAND
20   0 56.205g 0.023t  13008 S  99.7 18.7  10:00.83 java
20   0  188364   3932   2448 R   0.3  0.0   0:00.75 top
20   0   81452   1256   1012 S   0.0  0.0   0:00.00 pbzip2

I noticed the unpair=f and wasn't sure it was appropriate with my SE reads so I added unpair=t. Then the speed looked better yet it stuck at a different place:

Executing clump.KmerSplit [in1=in.fastq.gz, in2=null, 
out=out.fq.bz2, out2=null, groups=11, ecco=false, addname=f, 
shortname=f, [COLOR="Blue"]unpair=true[/COLOR], repair=f, namesort=f, ow=true, 
-Xmx16g, dedupe=t, addcount=t, dupedist=40, optical=t]

Input is being processed as unpaired
Made a comparator with k=31, seed=1, border=1, hashes=4
Time:                           54.743 seconds.
Reads Processed:      27914k    509.91k reads/sec
Bases Processed:       1423m    26.01m bases/sec
Executing clump.KmerSort3 [in1=in.clumped_clumpify_p1_temp%_2e77f340ad809301.fq.bz2, in2=null,
 out=out.clumped.fq.bz2, out2=null, groups=11, ecco=f, addname=false,
 shortname=f, [COLOR="red"]unpair=f[/COLOR], repair=false,
 namesort=false, ow=true, -Xmx16g, dedupe=t, addcount=t, dupedist=40, optical=t]

Making comparator.
Made a comparator with k=31, seed=1, border=1, hashes=4
Making 2 fetch threads.
Starting threads.
Fetching reads.
Exception in thread "Thread-57"
 *Control-C or similar caught [sig=15], quitting...
Exception in thread "Thread-58" Terminator thread: premature exit requested - quitting...
java.lang.RuntimeException: Duplicate process for file nutrinetat_vegb11_edi-0_r2y.clumped_clumpify_p1_temp0_2e77f340ad809301.fq.bz2
        at fileIO.ReadWrite.addProcess(ReadWrite.java:1599)
        at fileIO.ReadWrite.getInputStreamFromProcess(ReadWrite.java:1050)
        at fileIO.ReadWrite.getUnpbzip2Stream(ReadWrite.java:986)
        at fileIO.ReadWrite.getBZipInputStream2(ReadWrite.java:1086)
        at fileIO.ReadWrite.getBZipInputStream(ReadWrite.java:1066)
        at fileIO.ReadWrite.getInputStream(ReadWrite.java:802)
        at fileIO.ByteFile1.open(ByteFile1.java:261)
        at fileIO.ByteFile1.<init>(ByteFile1.java:96)
        at fileIO.ByteFile.makeByteFile(ByteFile.java:26)
        at stream.FastqReadInputStream.<init>(FastqReadInputStream.java:61)
        at stream.ConcurrentReadInputStream.getReadInputStream(ConcurrentReadInputStream.java:119)
        at stream.ConcurrentReadInputStream.getReadInputStream(ConcurrentReadInputStream.java:55)
        at clump.KmerSort3$FetchThread.fetchNext(KmerSort3.java:853)
        at clump.KmerSort3$FetchThread.run(KmerSort3.java:825)

When this happened, top looked like this,

PR  NI    VIRT    RES    SHR S  %CPU %MEM     TIME+ COMMAND
20   0  188332   3956   2452 R   0.3  0.0   0:04.02 top
20   0 23.961g 4.125g  12916 S   0.0  1.6   0:42.05 java
20   0   81452   1256   1012 S   0.0  0.0   0:00.01 pbzip2
20   0 1950332 173572   1056 S   0.0  0.1   0:03.36 pbzip2
20   0 1941568 173768   1052 S   0.0  0.1   0:03.49 pbzip2

Not sure if it's related but is the "unpair=f" expected during the clump, even though I added unpair=t in the command?

[Brian Bushnell]

1) pigz is in my path.
*Updated* When I monitored a running job, here's what I saw:

PR  NI    VIRT    RES    SHR S  %CPU %MEM     TIME+ COMMAND
20   0 56.006g 0.028t  13060 S 100.3 22.8 343:27.03 java
20   0   48424   1864   1324 S   0.0  0.0   0:00.00 pigz

It looked java is using all the requested CPUs. Also, the VM reached the limit of requested memory while the physical did not. I also confirmed that there's no limit on virtual memory. Any thought what might have caused the job stuck like this?

Oh... in rare cases, "reorder" can cause it to run very slowly, stuck at 100% CPU utilization. That might be the issue here... try removing that flag. During normal execution, pigz is also using CPU-time and java is usually substantially higher than 100%. What did the screen output look like when top was showing this?

2) I was looking for 'count=x' in the header thus didn't find anything. When I looked for 'copies=X' then I saw the duplicates.

Ooops, sorry, I updated my post to correct that!

[chiayi]

1) pigz is in my path.
*Updated* When I monitored a running job, here's what I saw:

PR  NI    VIRT    RES    SHR S  %CPU %MEM     TIME+ COMMAND
20   0 56.006g 0.028t  13060 S 100.3 22.8 343:27.03 java
20   0   48424   1864   1324 S   0.0  0.0   0:00.00 pigz

It looked java is using all the requested CPUs. Also, the VM reached the limit of requested memory while the physical did not. I also confirmed that there's no limit on virtual memory. Any thought what might have caused the job stuck like this?

2) I was looking for 'count=x' in the header thus didn't find anything. When I looked for 'copies=X' then I saw the duplicates.

Thank you!

[Brian Bushnell]

1) No, that's not expected, nor do I see a reason why that might be happening. Also, the performance is amazingly slow on your system in both cases. Here's what I get:

clumpify.sh in=reads.fq.gz out=clumped.fq.gz reorder

Time:                           40.496 seconds.
Reads Processed:      18059k    445.97k reads/sec
Bases Processed:       1824m    45.04m bases/sec

Reads In:           18059968
Clumps Formed:        809674
Total time:     40.680 seconds.

clumpify.sh in=reads.fq.gz out=clumped2.fq.gz reorder dedupe optical dupedist=40 addcount

Time:                           39.087 seconds.
Reads Processed:      18059k    462.04k reads/sec
Bases Processed:       1824m    46.67m bases/sec

Reads In:           18059968
Clumps Formed:        405473
Duplicates Found:        864
Total time:     39.135 seconds.

These are respectively about 400 and 4000 times faster than what you're seeing... I can only speculate that it might have something to do with network congestion... also, do you have pigz in your path? When you look at top while the process is running, what is the reported CPU utilization of java and pigz?

2) The count is only present for reads that had duplicates; duplicate-free reads won't have it. For example, it looks like this:

reformat.sh in=reads.fq.gz out=1pair.fq reads=1
cat 1pair.fq 1pair.fq 1pair.fq > 3pair.fq
clumpify.sh in=3pair.fq out=deduped_3pair.fq dedupe addcount
cat deduped_3pair.fq

@HISEQ13:296:CA8E8ANXX:7:1101:1432:1903 1:N:0:ACGGAACTGTTCCG copies=3
NTTGGATTGCATTAAGGAGCGAGGATGAGCATTCCATTCACCCGCTGGCCGGAAGAGTTTGCCCGTCGCTATCGGGAAAAAGGCTACTGGCAGGATTTGCC
+
!<=@=BFDGGGGGGBDDFGGGBGDGGFGGGGGGGGGGGGGGGGGGGGGBG@C/>EGGGGGGCG>DGGGGGGGCGGGGGGGGB>GGGGGGGGGGGGGCDGED
@HISEQ13:296:CA8E8ANXX:7:1101:1432:1903 2:N:0:ACGGAACTGTTCCG copies=3
TCGTTGAGCAGTTGCACCACGCGAATGGAGGAATGTTCTGTGACGAAAGTATTGAGGAAATCATCCCCGCTAAACAGCGCATGTTGGCGATCGGCAATCAG
+
BBBBBF0=1;@FGGDDGEGGDGD/EGGG:>E/BGGFDGGGGG@FGBEFGGEGGGGGDCGGCDGGCCCGGBGGDDGGEGGD>GGGGGGFD8?D<FGDG/EGC

[GenoMax]

2) For the output, where is the 'count=x' supposed to be stored? I did not see it in the out fq file.

It is stored in the fastq header in the "clumped" file. Only where clustered reads were found.

3) I would like to get an idea about the ratio of optical duplicates vs PCR duplicates. To do that, should I run the first time with optical=t and the second time with optical=f, and then compare the difference between these two runs to get the PCR duplicates?

Thanks a lot!

dedupe=f optical=f (default)
Nothing happens with regards to duplicates.

dedupe=t optical=f
All duplicates are detected, whether optical or not.  All copies except one are removed for each duplicate.

dedupe=f optical=t
Nothing happens.

dedupe=t optical=t

Only optical duplicates (those with an X or Y coordinate within dist) are detected.  All copies except one are removed for each duplicate.
The allduplicates flag makes all copies of duplicates removed, rather than leaving a single copy.  But like optical, it has no effect unless dedupe=t.

You should add "spantiles=f" for all non-NextSeq data. For optical duplicates you would need to specify "dupedist=". There are suggested values in this thread.

[chiayi]

Note that this does NOT work in "markduplicates" mode, only in "dedupe" mode. Reads that absorbed duplicates will have "count=X" appended to the end to indicate how many reads they represent (including themselves, so the minimum you would see is 2).

Hi Brian,

Thank you for adding these new features. I tried clumpify (37.17) to detect duplicates in the HiSeq 2500 single-ended strand-specific 50 bp reads. My pipeline is Clumpify-BBDuk-BBMap-featureCounte-a) edgeR for differential expression analysis or b) network construction. I got some qestuions about performance and output after several trials on dedup.

1) I noticed that when I added dedupe=t, the speed slowed down for ~10 time. Is this expected or did I set some options incorrectly?

No dedup:

clumpify.sh -Xmx44g in=$INPUT out=$OUTPUT.clumped.fq.gz reorder=t overwrite=t

Time:                           3540.405 seconds.
[COLOR="red"]Reads Processed:      27914k    7.88k reads/sec[/COLOR]
Bases Processed:       1423m    0.40m bases/sec

Reads In:           27914336
Clumps Formed:       2997579
Total time:     3540.575 seconds.

Remove optical duplicates:

clumpify.sh -Xmx48g in=$INPUT out=$OUTPUT.clumped.fq.gz reorder=t overwrite=t \
dedupe=t addcount=t dupedist=40 optical=t

Time:                           38030.384 seconds.
[COLOR="Red"]Reads Processed:      27914k    0.73k reads/sec[/COLOR]
Bases Processed:       1423m    0.04m bases/sec

Reads In:           27914336
Clumps Formed:       2997579
Duplicates Found:      27931
Total time:     41299.460 seconds.

2) For the output, where is the 'count=x' supposed to be stored? I did not see it in the out fq file.

3) I would like to get an idea about the ratio of optical duplicates vs PCR duplicates. To do that, should I run the first time with optical=t and the second time with optical=f, and then compare the difference between these two runs to get the PCR duplicates?

Thanks a lot!

[GenoMax]

Wonderful!

Will give it a try as soon as I can get bbmap upgraded to latest on the cluster.

[Brian Bushnell]

I also wanted to see counts (with associated sequence) to see how acute of a problem the duplicates may be.

Update! Clumpify can now report the number of duplicates per read. Usage:

clumpify.sh in=file.fq out=deduped.fq dedupe addcount

Note that this does NOT work in "markduplicates" mode, only in "dedupe" mode. Reads that absorbed duplicates will have "copies=X" appended to the end to indicate how many reads they represent (including themselves, so the minimum you would see is 2).

[dpryan]

Since inevitably others will run into this and wonder...

Please note that if you run FastQC on files run through clumpify, that you will see bias in the reported duplication rates between read #1 and read #2 in a pair. In other words, the duplication rate returned by FastQC for read #2 will be much much much (often 2-3x) higher than that of read #1. Note that this is a technical artifact due to FastQC's duplication module only looking at the first 100,000 reads in each file. Since the files can be reordered for increased compressibility, the first 100,000 reads are then expected to be similar and the results skewed.

In other words, don't worry about such observations.

[GenoMax]

A blog post about duplicates on HiSeq 4000 has been posted to QCFail site (with pretty graphics): https://sequencing.qcfail.com/articl...ated-sequences

[Stewart Russell]

by default BBMerge won't look for insert sizes shorter than 35, and by default Clumpify will not work on reads shorter than 31bp, so set those as appropriate. If you are not expecting insert sizes shorter than 35bp then you can remove those flags.

This explains why I was losing so many reads from BBmerge. I should have read the doc more carefully!

But, if you want to deduplicate the raw reads (though I'd really recommend adapter-trimming first), and ignore pairing, you can add the "unpair" flag to Clumpify

Given that I'm looking to merge/remove adapters, and then collapse based on sequence + 5' and 3' 4N, allowing for 3 non-templated bases, I think I need to do

bbmerge.sh in=r1.fq in2=r2.fq out=merged.fq outu=unmerged.fq mininsert=20
clumpify.sh in=merged.fq out=merged_clumped.fq dedupe dedupesubs=3 k=20 unpair

Note that after doing this pairing order will be lost so you'd need to run repair.sh to recover it.

At this point I don't think I need to re-pair because I have my data in a clean, single file for mapping.

Thanks for the help

[Brian Bushnell]

My question is whether clumpify respects the read names and will only remove paired duplicates, or if it can remove all duplicates in the data.

Clumpify does respect pairs and by default will only consider pairs duplicates if both reads match. In this situation it makes the most sense to me to merge reads first, then remove duplicates on the merged reads:

bbmerge.sh in=r1.fq in2=r2.fq out=merged.fq outu=unmerged.fq mininsert=20
clumpify.sh in=merged.fq out=merged_clumped.fq dedupe k=20

Merging will generally remove any adapter bases on the ends because those are not part of the overlapping sequence. Also, since merging removes bases off the ends of the adapters too, there will be fewer mismatches to confuse the issue. Note in this case I set "mininsert=20" and "k=20" because I don't know what your minimum expected length is for small RNAs, but by default BBMerge won't look for insert sizes shorter than 35, and by default Clumpify will not work on reads shorter than 31bp, so set those as appropriate. If you are not expecting insert sizes shorter than 35bp then you can remove those flags.

But, if you want to deduplicate the raw reads (though I'd really recommend adapter-trimming first), and ignore pairing, you can add the "unpair" flag to Clumpify:

clumpify.sh in1=r1.fq in2=r2.fq out=clumped.fq dedupe k=20 unpair

As Genomax mentioned, you will probably need to play with the "dupesubs" flag if there are expected mismatches due to nongenomic bases. Note that after doing this pairing order will be lost so you'd need to run repair.sh to recover it.

[GenoMax]

If you have extra bases at the 3' and 5'-ends then they may complicate duplicate identification (you can add dupesubs=4) when running clumpify.

Here are combinations of flags for clumpify that may be useful.

dedupe=f optical=f (default)
Nothing happens with regards to duplicates.

dedupe=t optical=f
All duplicates are detected, whether optical or not.  All copies except one are removed for each duplicate.

dedupe=f optical=t
Nothing happens.

dedupe=t optical=t

Only optical duplicates (those with an X or Y coordinate within dist) are detected.  All copies except one are removed for each duplicate.

The allduplicates flag makes all copies of duplicates removed, rather than leaving a single copy.  But like optical, it has no effect unless dedupe=t.

[Stewart Russell]

Clumpify - dedup fasta/fastq

Hey Brian, first off, thanks for designing these great tools and providing thorough explanation.

I'm trying to de-duplicate small RNA seq reads from a 2x150 bp MiSeq run. The library prep for ultra-low input has made the adapter sequences less than straight forward so I've removed those prior to attempting to merge and remove duplicates.

My question is whether clumpify respects the read names and will only remove paired duplicates, or if it can remove all duplicates in the data. My reads have 3' and 5' 4N sequences (to reduce bias in adapter ligation) that can be used to identify PCR duplicates if I can somehow collapse to unique reads before trimming them.

Thanks in advance for your suggestions,

Stewart