bbduk.sh masking shorter than k-mer size

[cstack]

I've come across a situation using bbduk to mask the k-mers from one genome assembly in another where the resulting number of masked bases is some sequences is shorter than the k value.

I am running the command like this:

# version 37.something.  
bbduk.sh in=mygenome.fa out=mygenome_masked.fa ref=ecoli.fasta k=15 
qkmask=X maskfullycovered=t maskmiddle=f

In mygenome_masked.fa (a multisequence fasta file), there are a sizable number of sequences with a total of 0 < bases_masked < k(=15). It seems strange to have 3 nucleotides masked when k=15, and I am wondering if anyone can point out what options I should be using to prevent this from happening

[Brian Bushnell]

That would be the "maskfullycovered" flag. That means only bases covered entirely by reference kmers will be masked.

For example, take this sequence you want to mask:

ACGTTGCA

And this reference:

CGTTG

The ref kmers (ignoring reverse complement, at K=3) are CGT, GTT, and TTG. They line up like:

   TTG
  GTT
 CGT
ACGTTGCA

Every based is covered by 3 kmers, but only the first T is "fully covered" - covered by 3 ref kmers. So it's the only one masked. Whereas without "maskfullycovered", the entire CGTTG would be masked.

Incidentally, it depends on what your goal is, but normally I find K=15 to be very short for masking... typically I use K=31.

[cstack]

Thanks! I had been thinking about maskfullycovered incorrectly.

Incidentally, it depends on what your goal is, but normally I find K=15 to be very short for masking... typically I use K=31.

The sequences I was trying to mask were from a very rough (initial) assembly of 10-20x worth of uncorrected 1d Minion reads. The assembly probably has a fair few errors, so I tried a range of k values to see if I could find any obvious contamination. K=15 turned out to be way too sensitive, like you say.

Are there parameter values you'd recommend (for either bbduk or seal) for identifying contaminants in high-error-rate (~15%) reads?