[Brian Bushnell]

You can remove all the reads with any 1-count kmers using BBNorm if you really want to, by messing with some of the parameters, but that's not really what it's designed for. I should probably make that more straightforward. It's easier to do like this:

kcompress.sh in=file.fq out=kmers.fa minprob=0 mincount=1 maxcount=1
bbduk.sh in=file.fq out=filtered.fq ref=kmers.fa k=31 mm=f

However, there will still be singleton kmers. There will always be singleton kmers, because when you remove reads with singleton kmers, that creates new singleton kmers. For example, if a read has 50 1-count kmers and 50 2-count kmers, then after you remove it, you eliminate 50 1-count kmers but generate 50 new 1-count kmers. You can reduce them, but they never go away.

[horvathdp]

incomplete removal of singletons by BBNORM

Brian (or others),

I am still a bit disturbed by the inability to remove reads that contain singleton kmers. I understand that BBNORM generates an output file that has all reads that contain kmers between specific values, and that in doing so the output file still may contain reads with singleton kmers- thus the khist runs of output BBNORM files still may have a large number of singletons present. However, it seems that it should be possible to generate a list of all reads that contain singleton kmers and mask or delete them from the output - thus eliminating any fragment that contains any low-quality kmers. Note being a programer, is there a reason why this can not (or should not) be done? It seems to me that eliminating such reads might greatly simplify the graphs during any subsequent assembly.

[lac302]

Thanks for the quick reply Brain. I'll download the updated version.

[Brian Bushnell]

Perhaps I'm using an older version of bbmap but usage shows:
[CODE]Usage: dedupe.sh in=<file or stdin> out=<file or stdout>

An example of running Dedupe for clustering short reads:
dedupe.sh in=x.fq am ac fo c rnc=f mcs=4 mo=100 s=1 pto cc qin=33 csf=stats.txt pattern=cluster_%.fq dot=graph.dot

Input may be stdin or a fasta, fastq, or sam file, compressed or uncompressed.

Yep, that would be an older version, before I fixed that comment. Dedupe has never been able to process sam files, but I didn't realize it when I first wrote it, since I had never tried it. The problem is that BBTools Read objects have a single special field that can hold attachments; processing sam input uses that special field to hold the original sam line, but Dedupe needs it for storing other information.

WGS, 2x275 PE reads from a Nextera XT library. Looking at novel non-syn SNPS found in treatment group vs WT.

You can dedupe the fastqs before mapping or the sam files (by mapping) after mapping, then. Probably best to deduplicate by mapping in this case because 2x275 bp reads will have a high error rate so deduplication by sequence will be less effective and concentrate reads with errors.

[lac302]

Dedupe can't use sam files.

Perhaps I'm using an older version of bbmap but usage shows:

Usage:	dedupe.sh in=<file or stdin> out=<file or stdout>

An example of running Dedupe for clustering short reads:
dedupe.sh in=x.fq am ac fo c rnc=f mcs=4 mo=100 s=1 pto cc qin=33 csf=stats.txt pattern=cluster_%.fq dot=graph.dot

Input may be stdin or a fasta, fastq, or sam file, compressed or uncompressed.
Output may be stdout or a file.  With no output parameter, data will be written to stdout.
If 'out=null', there will be no output, but statistics will still be printed.
You can also use 'dedupe <infile> <outfile>' without the 'in=' and 'out='.

If you want to deduplicate based on mapping, I recommend converting to bam, sorting, and using some other program (Picard, for example). Or you can use dedupebymapping.sh which DOES accept sam files. It's much faster than Picard but does not fully take advantage of sorting, or write temp files, so it uses a lot of memory (no more than Dedupe, though).

Dedupe is for deduplicating by sequence. So, whether you deduplicate the fastq then map, or map and deduplicate the sam, depends on your goal - the results will differ slightly. In practice, it won't usually matter much for paired reads, but single-ended reads are much more impacted by deduplication methodology.

What kind of data are you using, and what's your experiment?

WGS, 2x275 PE reads from a Nextera XT library. Looking at novel non-syn SNPS found in treatment group vs WT.

[Brian Bushnell]

If you want to deduplicate based on mapping, I recommend converting to bam, sorting, and using some other program (Picard, for example). Or you can use dedupebymapping.sh which DOES accept sam files. It's much faster than Picard but does not fully take advantage of sorting, or write temp files, so it uses a lot of memory (no more than Dedupe, though).

Dedupe is for deduplicating by sequence. So, whether you deduplicate the fastq then map, or map and deduplicate the sam, depends on your goal - the results will differ slightly. In practice, it won't usually matter much for paired reads, but single-ended reads are much more impacted by deduplication methodology.

What kind of data are you using, and what's your experiment?

[GenoMax]

sam files are an accepted input file type for dedupe.sh correct? Or would you recommend removing duplicates from the fastq files?

Dedupe can't use sam files.

Input may be fasta or fastq, compressed or uncompressed.

[lac302]

sam files are an accepted input file type for dedupe.sh correct? Or would you recommend removing duplicates from the fastq files?

[root@g300-149-b0 bbmap]# ./dedupe.sh Um57.sam Um57_nodups.sam 
Max memory cannot be determined.  Attempting to use 3200 MB.
If this fails, please add the argument -Xmx29g (adjusted to roughly 85 percent of physical RAM).
java -Djava.library.path=/app/bbmap/jni/ -ea -Xmx3200m -Xms3200m -cp /app/bbmap/current/ jgi.Dedupe Um57.sam Um57_nodups.sam
Executing jgi.Dedupe [Um57.sam, Um57_nodups.sam]

Initial:
Memory: free=3182m, used=34m

Exception in thread "Thread-3" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-14" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-13" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-11" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-12" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-10" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-9" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-8" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-7" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-5" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-4" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Exception in thread "Thread-6" java.lang.ClassCastException: stream.SamLine cannot be cast to jgi.Dedupe$Unit
	at jgi.Dedupe$HashThread.processRead(Dedupe.java:3067)
	at jgi.Dedupe$HashThread.processReadOuter(Dedupe.java:3045)
	at jgi.Dedupe$HashThread.run(Dedupe.java:2980)
Found 0 duplicates.
Finished exact matches.    Time: 0.072 seconds.
Memory: free=2930m, used=286m

Found 0 contained sequences.
Finished containment.      Time: 0.006 seconds.
Memory: free=2712m, used=504m

Removed 0 invalid entries.
Finished invalid removal.  Time: 0.001 seconds.
Memory: free=2712m, used=504m

Input:                  	12 reads 		3295 bases.
Duplicates:             	0 reads (0.00%) 	0 bases (0.00%)     	0 collisions.
Containments:           	0 reads (0.00%) 	0 bases (0.00%)    	0 collisions.
Result:                 	0 reads (0.00%) 	3295 bases (100.00%)

Printed output.            Time: 0.003 seconds.
Memory: free=2695m, used=521m

Time:   			0.089 seconds.
Reads Processed:          12 	0.14k reads/sec
Bases Processed:        3295 	0.04m bases/sec
[root@g300-149-b0 bbmap]

[Brian Bushnell]

BBMap has no SNP-calling functionality, so you'll need to use a variant-caller (samtools, GATK, FreeBayes, etc).

[lac302]

Thanks for the quick reply Brian...Does bbmap have any SNP calling functionality? Or should I just feed these .sam files into the samtools pipeline?

[lac302]

Is the insert size measured or inferred?

measured via fragment analyzer

[GenoMax]

Is the insert size measured or inferred?

[Brian Bushnell]

Oh, that's probably fine. You'd have a much higher error rate if this was a Nextera LMP library that had not been correctly preprocesed. As you can see, there is a 0.4% per-base substitution rate for read 1, which is pretty good for such long reads. Read 2 is higher, as expected. It's normal for long reads to have a high fraction containing at least one substitution.

[lac302]

Ouput below. This is a nextera XT fragment library that I made and sequenced. Average insert size is around 900bp. I used the defualut bbmap settings adding qtrim=r and pairlen=2000

BBMap version 34.94
Retaining first best site only for ambiguous mappings.
Executing dna.FastaToChromArrays2 [/data/U.maydis/GCA_000328475.2_Umaydis521_2.0_genomic.fna, 1, writeinthread=false, genscaffoldinfo=true, retain, waitforwriting=false, gz=true, maxlen=536670912, writechroms=false, minscaf=1, midpad=300, startpad=8000, stoppad=8000, nodisk=true]

Set genScaffoldInfo=true
Set genome to 1

Loaded Reference:	0.005 seconds.
Loading index for chunk 1-1, build 1
Indexing threads started for block 0-1
Indexing threads finished for block 0-1
Generated Index:	1.753 seconds.
Analyzed Index:   	3.196 seconds.
Started output stream:	0.018 seconds.
Cleared Memory:    	0.180 seconds.
Processing reads in paired-ended mode.
Started read stream.
Started 12 mapping threads.
Detecting finished threads: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11

   ------------------   Results   ------------------   

Genome:                	1
Key Length:            	13
Max Indel:             	16000
Minimum Score Ratio:  	0.56
Mapping Mode:         	normal
Reads Used:           	10712672	(2767819477 bases)

Mapping:          	2375.634 seconds.
Reads/sec:       	4509.39
kBases/sec:      	1165.09


Pairing data:   	pct reads	num reads 	pct bases	   num bases

mated pairs:     	 95.9222% 	  5137917 	 95.9216% 	  2654935852
bad pairs:       	  0.3883% 	    20798 	  0.4112% 	    11380974
insert size avg: 	  457.71


Read 1 data:      	pct reads	num reads 	pct bases	   num bases

mapped:          	 97.7238% 	  5234416 	 97.8595% 	  1353706554
unambiguous:     	 96.7062% 	  5179910 	 96.8763% 	  1340106299
ambiguous:       	  1.0176% 	    54506 	  0.9832% 	    13600255
low-Q discards:  	  0.0474% 	     2537 	  0.0064% 	       88795

perfect best site:	 62.2963% 	  3336797 	 61.0063% 	   843909803
semiperfect site:	 62.3484% 	  3339588 	 61.0594% 	   844644584
rescued:         	  0.0469% 	     2511

Match Rate:      	      NA 	       NA 	 97.9361% 	  1345904929
Error Rate:      	 36.1452% 	  1892020 	  2.0011% 	    27499816
Sub Rate:        	 35.0788% 	  1836199 	  0.4172% 	     5733414
Del Rate:        	  2.4959% 	   130648 	  1.4962% 	    20562060
Ins Rate:        	  1.7790% 	    93120 	  0.0876% 	     1204342
N Rate:          	  0.4722% 	    24718 	  0.0629% 	      863869


Read 2 data:      	pct reads	num reads 	pct bases	   num bases

mapped:          	 96.5932% 	  5173855 	 96.6576% 	  1338226877
unambiguous:     	 95.5860% 	  5119909 	 95.6852% 	  1324764021
ambiguous:       	  1.0071% 	    53946 	  0.9724% 	    13462856
low-Q discards:  	  0.0478% 	     2561 	  0.0067% 	       92659

perfect best site:	 28.1750% 	  1509145 	 24.9411% 	   345310812
semiperfect site:	 28.2006% 	  1510517 	 24.9662% 	   345657705
rescued:         	  0.0547% 	     2932

Match Rate:      	      NA 	       NA 	 96.0340% 	  1307794049
Error Rate:      	 70.7813% 	  3662191 	  3.9015% 	    53130664
Sub Rate:        	 70.3060% 	  3637596 	  2.0494% 	    27908811
Del Rate:        	  2.8907% 	   149564 	  1.7313% 	    23576787
Ins Rate:        	  2.7860% 	   144147 	  0.1208% 	     1645066
N Rate:          	  0.5726% 	    29628 	  0.0645% 	      878951

Total time:     	2381.122 seconds.

[Brian Bushnell]

BBMap reports error rates in two pairs of columns; the first is the total number and fraction of reads with any errors, the second is the total number and fraction of bases that are errors. Typically, for Illumina data, the most important number is be the per-base substitution rate, which is hopefully under 3% or so, but it can get pretty high with such long reads (I assume this was a MiSeq 2x300bp kit). Posting the entire screen output here might be helpful.

Also, note that there are (at least?) two completely different Nextera protocols, one for fragment libraries and one for long-mate pair libraries. Since you have an insert range of 400-2000bp, I assume you are using Nextera long-mate pairs. Is that correct? These cannot be mapped directly, but need to be preprocessed first, as they have a chimeric junction somewhere. To do that, first adapter-trim using BBDuk, then run splitnextera.sh on them:

splitnextera.sh in=<file> out=<file> outf=<file> outu=<file> outs=<file> mask=t

This will remove the junction adapters and split the reads into multiple output files which should be mapped independently. out is for long-mate pairs; outf is for short fragments; outu is for unknown orientation pairs; and outs is for singletons. If nothing comes out as long mate pairs, then you don't have Nextera LMP data and should not perform this step. It's best to check with the people who prepared the library, though.