Edit to add info:

I thought maybe the loss of data at the Part A: quality recalibration ('Step 4'; 800 MB lost) might be due to some filtering of very low depth kmers. (e.g. a long tail of kmer depth=1 for plant rhizosphere data likely exists in the raw data; and this was being filtered out somehow).

To test, I used khist.sh:

khist.sh in="$IN" hist="$OUT"/"$f"_kmer-hist-zeroes.tsv peaks="$OUT"/"$f"_kmer-hist-peaks-zeroes.tsv minprob=0 minqual=0 mindepth=0 printzerocoverage=true minq=6 minprob=0.5

I am not that experienced at these analyses, so I was surprised to learn that the low-depth coverage was identical between raw and qual-recalib kmers for this sample. Instead, the quality recalibration seemed to have shaved off the higher depth kmers, if I'm interpreting the results correctly.

Specifically, the highest kmer depth values recorded in khist (tail) were about twice as deep for the RAW data:

         #Depth    Raw_Count    Unique_Kmers  
RAW:    2426           3703        2  
RECAL: 1228        2372       2

The loss of high-depth kmers reminds me of read normalization, which is commonly done before assembly, but I do not want read normalization to happen at this quality recalib step.

Here again is the code in the polyfilter pipeline that occurs during the quality recalibration step. Note that I did not include the usetiles flag in my pipeline because these simulated reads do not have tile info.

#Make recal matrices
calctruequality.sh "$HIGH" "$ARGS" in=clean.sam.gz usetiles ref=quick.fa callvars

#Recalibrate quality
bbduk.sh "$LOW" "$ARGS" in=clean.fq.gz out=clean_recal_tile.fq.gz recalibrate usetiles

Here's the kmer hist comparison, in case others have additional comments.

Tagging Brian Bushnell (author of BBTools).

Edit to add two more pieces of info:

Reading the bbduk.sh params again, I found this:

Note - if ktrim, kmask, and ksplit are unset, the default behavior is kfilter.
All kmer processing modes are mutually exclusive.

The defaults for these params are here:

qtrim=f             Trim read ends to remove bases with quality below trimq.
ksplit=f            For single-ended reads only.  Reads will be split into
kmask=              Replace bases matching ref kmers with another symbol.
                    Allows any non-whitespace character, and processes short
                    kmers on both ends if mink is set.  'kmask=lc' will
                    convert masked bases to lowercase.
    k=31

So, my assumption that bbduk.sh is running quality score recalibration only is wrong, it seems? It is also running a kfilter by default, for this quality score recalibration step.

To test this theory, I reran the entire polyfilter pipeline as I had done previously, but adding outm=recal-matched.fq.gz to the quality recalib step, as shown here:

bbduk.sh -Xmx146553m -Xms200m in=clean.fq.gz out=clean_recal.fq.gz outm=recal-matched.fq.gz recalibrate;

And re-checked the file sizes:

#SIZE     #STEP       #FILE NAME
4.4G    INFILE A    Anonymous-reads.fq.gz
    …   
3.6G    Step 4  clean_recal.fq.gz
28       Step 4 recal-matched.fq.gz

So, I am not able to explain this change from 4.4GB to 3.6GB with kfiltering, it seems, since the outmatched file is basically empty.

Second piece of info: checking difference in total number of reads, rather than file size. I used lhist for this and though I ran it twice, I am assuming I am making some errors because somehow I end up with slightly more reads after running the quality recalibration:

RAW anonymous_reads.fq.gz_readLen-hist.tsv
Length  Count
135 146
136 267
137 587
138 1262
139 2785
140 5557
141 10487
142 15919
143 16096
144 16049
150 33175769
SUM:    33,244,924

QUAL-RECALIB     clean_recal.fq.gz_readLen-hist.tsv
Length  Count
150 33,247,226

So, the 800 MB smaller file size for clean_recal.fq.gz may just be more efficient file compression? But confused with the change in read length distribution and read number from raw to recalibrated. Also, maybe the difference in kmer depth is just due to a few outlier reads?

Summary of read counts using read length histograms, lhist.

# Name                 File Size        read len range, bp    Sum read counts
anonymous_reads.fq.gz     4.4GB          135-150 bp           33,244,942
clean_recal.fq.gz         3.6GB            150 bp             33,247,226
histd2.fq.gz              3.6GB           135-150 bp         33,244,942
ecco.fq.gz                1.3GB           15 - 150 bp         10,590,400
merged_both.fq.gz         712MB           15- 291 bp         3,660,983 (7,321,966)

So, it seems like the only big chop of 2/3 of my reads is happening in the error correction step, to create ecco.fq.gz in the first step of Part B. I assume the adapters.fa file is the same as the one generated earlier in the pipeline.

The information for prefilter tells me that changing this setting will conserve my long tail of kmers at depth=1. I will check by re-running with prefilter=0 and compare. This adjustment to prefilter=0 seems like it will help us recover more rare taxa when we start doing co-assemblies, rather than individual sample assemblies, as long as we can wrangle enough RAM for the larger kmer tables.

RESULT: ecco.fq.gz is still 1.3GB when prefilter=0.

bbmerge.sh "$MAX" "$ARGS" in=hdist2.fq.gz out=ecco.fq.gz ecco mix adapters=adapters.fa kfilter=1 k=31 prefilter=1

prefilter=0          If set to a positive integer, use a countmin sketch to 
                     ignore kmers with depth of that value or lower, to
                     reduce memory usage.

The last parameter that caught my eye for this step is kfilter=1.

kfilter=0            Ban overlaps that create kmers with count below
                     this value (0 disables).  If this is used minprob should
                     probably be set to 0.  Requires good coverage.

minprob=0.5          Ignore kmers with overall probability of correctness 
                     below this, to reduce memory usage.

I reran with prefilter=0, kfilter=0 and minprob=0.

bbmerge-auto.sh -Xmx207642m -Xms200m t=2 in="$IN"/hdist2.fq.gz out="$OUT"/ecco.fq.gz ecco mix adapters="$IN"/adapters.fa kfilter=0 k=31 prefilter=0 minprob=0

RESULT: ecco.fq.gz is still 1.3GB.

I can see from B Bushnell et al 2017 paper announcing BBmerge, in Fig 1 panel B for example, when reads do not overlap, they are discarded. But, I thought with the ecco flag in place, there is error correction by overlap only. So, esp with filters set to zero, there should be no discarded reads at this step, regardless of whether they overlap or not, correct? The actual outfile of merged reads happens later in the pipeline.