Single base difference in overlapping R1 & R2 reads: Trust single higher quality basecall in a general lower quality stretch?

Question

[Carambakaracho]

Hi there,

I have some 16 S Amplicon seq data and sequenced roughly 460-480 bp amplicons using 2x300 MiSeq reads. Usually, I use flash2 to merge the reads after quality trimming. However, recently I stumbled across this protocol from Fred Mahé, which suggests to merge reads before trimming. This sort of made sense to me as trimming effects mostly the overlapping ends. However, the difference surprised me nevertheless, and in some cases doubled the number of reads that were merged.

Now I'm trying to understand some pieces of the puzzle in more detail - one is this example (end of R1 on top, R2 below):

AAAGAGGAATTGCTCCAGGAGGGGTCTGCGTCCGATTAGCTGGATGGTGGGGTAAGGGCCTACCATGGCGATGATCGGTAGCTGGTCGGAGAGGATGAGCAGCC
    AGGAATTGCTCCAGGAGGGGTCGGCGTCCGATTAGCGGGATGGTGGGGTAAAGGCCTACCATGGCGATGATCGGTAGCTGGTCTGAGAGGATGATCAGCCACA
                                                       ^

Where the respective quality string are (Illumina 1.8+ Phred+33)

EGGGGGGFFD?<FFGFGGGGGGGGGGFGGGCEDDCE5CFFGFF8FDFC?EEG@DF*6>3:)23>79FB?FFB:?4631?BFBFGF(-0435(.48:?0)-<C?0
    ########################################??>::4360*).02?C4@7:5(1?79FD83<D@=36GGGDGFF:;22)AGFGFF;EF>GAC=AE
                                                       ^

Based on typical Phred+33 scores, the . is better than * at the exact position, so in untrimmed sequences I get an A call at the position - when merging the trimmed reads, I removed the R2 end because of it's general low quality and hence retain the G. This sequence section covers a hypervariable region and single nucleotide differences in the right position can lead to drastic changes in the classification.

Would you trust the higher quality single base call in a lower quality stretch more than the low quality base call in an otherwise higher quality region?

[Bastien Hervé]

As I never done this kind of process, this will just be my personal advise :

I would take a look at the neighbour qualities of your base of interest (5 bases on each side)

R1 : GGTAA`G`GGCCT        Quality : EG@DF`*`6>3:)      Score : 36/38/31/35/37/`9`/21/29/18/25/8
R2 : GGTAA`A`GGCCT        Quality : 360*)`.`02?C4         Score : 18/21/15/9/8/`13`/15/17/30/34/19

R1 got a better quality mean overall but the machine is really not sure about this 'G' whereas it calls with high score the neighbouring bases. I would not trust this base for classification

R2 got heterogeous scores from good to very bad call, your 'A' base is 13 Q-score which is not that good 1 in 30 to be a bad call...

I won't believe either of the two bases for the classification. Anyhow I guess your classification is mostly based on base qualities so the soft will not used these bad bases.

[Carambakaracho]

I won't believe either of the two bases for the classification.

hm, matches my gut feeling.

Anyhow I guess your classification is mostly based on base qualities so the soft will not used these bad bases.

unfortunately not - it's based on the consensus sequence. I convert to fasta. Maybe I include a conversion to ambiguous nucleotides for low quality sequences before conversion to fasta

[Bastien Hervé]

Transform base to N if both quality below a specific threshold could be an idea