Why does cutadapt remove low quality bases from the ends of reads only?

Question

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?

Answer 1

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.

Answer 2

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).

Answer 3

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.