I use cutadapt to remove low quality bases from my Illumina reads. The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed. Why? Why doesn't the algorithm remove all low quality bases? It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?
I am commenting on this part:
The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed.
According to its user guide, cutadap is designed this way: it trims off bases from the 3'-end until it sees a base with quality higher than a threshold. This is not a good algorithm. For example, when we see a quality string
30,30,30,3,3,3,3,3,20,3,3, we would prefer to trim it down to
30,30,30 because a
3,3,3,3,3,20 tail is still not useful.
Some tools such as trimmomatic uses a sliding window to avoid this issue. In my view, a better algorithm is the so-called the "modified Mott algorithm" used by phred 20 years ago. seqtk among others implements this algorithm. Note that the Mott algorithm always trims from both ends, which is sometimes not preferred. BWA implements a variant of this algorithm, trimming from the 3'-end only.
In practice, though, different quality trimming algorithms probably work equally well because it is relatively rare to see a high-quality base in the middle of a low-quality tail. In addition, modern aligners can well handle low-quality tails; assemblers and error correctors can correct through such tails, too. It is usually not necessary to apply quality trimming.
It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?
With a low-quality base, you might have an error/mismatch, but with an N, you always have a mismatch. Masking low-quality bases to N is not as good.
EDIT: responding to the following comment from OP:
What if the base is low-quality and then it matches the reference by mistake? Wouldn't it be better if it did not match and was penalized? I guess this should be modeled mathematically.
If a read could be mismapped due to a sequencing error, its true location is often in the genome. In this case, a competent mapper will give the alignment a low mapping quality (mapQ). A quality-aware mapper such as novoalign will penalize mapQ even more. In comparison, if you hard mask this sequencing error to "N", you will get a mapQ=0. You can see the difference between the two approaches mostly comes from mapQ.
In modern Illumina data, it is quite often to see a Q8 base in the middle of high-quality bases. >80% of them (in theory) are still correct. My hunch is that hard masking all of them would lead to considerable data loss.
$\begingroup$ Is it possible to expand on the last part? 'With a low-quality base, you might have an error/mismatch, but with an N, you always have a mismatch. Masking low-quality bases to N is not as good.' What if the base is low-quality and then it matches the reference by mistake? Wouldn't it be better if it did not match and was penalised? I guess this should be modelled mathematically. $\endgroup$ Feb 26, 2018 at 16:34
If you check the read QC statistics of an Illumina run in e.g. fastQC, you will see that at the end of the read the quality decreases. This is because of exhaustion of chemicals at the end of the run. This is a general trend seen in all runs, therefore you can remove these low quality bases from the end of your run. If you have incidentally a bad quality base in the middle of a read this is not the same general trend but is sporadic. Which is hard to automatically remove with tools such as cutadapt. If you really want an N for each low quality bases within the read, you can make a simple program for that, but why would you need that? Most aligners can handle a few bad quality bases within a read, and for variant calling you'll need many reads (so the incidental error base will be overruled by the good ones).
2$\begingroup$ It's also worth noting that many aligners penalize an N by default, so it's not like replacing super-low quality bases with an N is going to greatly boost alignment scores. $\endgroup$ Feb 26, 2018 at 12:47
For Illumina (and 454) reads, the quality decreases with read length. It's not linear and is run/library-dependent. It has less to do with exhaustion of reagents and more to do with the strands on a spot being out of phase due to incomplete/missed base incorporation during sequencing.
It is common practice to trim off the low quality 3' ends as the first step in QC. Most analyses and most can deal with a few low errors in the middle of the read and what you do about such sequencing errors may vary depending on what type and use of the data. But in general, nobody wants the garbage at the end of the reads.
P.S. Even for the same library prep, the fwd/rev runs will have different quality curves. P.P.S. Sometimes the first couple bases are also poor too.
$\begingroup$ "in general, nobody wants the garbage at the end of the reads." For RNA-seq, quality trimming may be ok, but for de novo assembly, trimming may reduce N50. $\endgroup$ Mar 1, 2018 at 14:24
$\begingroup$ If you don't trim prior to assembly then a DBG assembler will have a lot of extra kmers, blowing up RAM. It's been my experience that trimming with increase N50, not decrease it. Although I don't think N50 is a great metric to evaluate assembly quality. I suppose as an alternative to end-trimming, one could do error correction, e.g. with bfc. $\endgroup$ Mar 1, 2018 at 23:24
1$\begingroup$ Modern assemblers (e.g. allpaths, sga and spades) often come with an error correction step. They only trim an end when they can't correct it. Only a small fraction of rubbish k-mers can go into the graph. $\endgroup$ Mar 1, 2018 at 23:50
$\begingroup$ Thanks, user172818, you're probably right. But we have found it beneficial to use bfc prior to metaspades. $\endgroup$ Mar 3, 2018 at 0:36