potential bug on ONT uncorrected Q20 output

[colindaven]

Hi,

thanks for this nice tool which I have been enjoying using. The side by side QC bar charts and tables are very nice indeed.

I am getting strange results with current ONT Q20 reads (without herro/dorado correct read correction). These should have Q values of Q10-30 (peak at Q25 ?) or similar in a distribution. However we are seeing reports with the modal peak distribution around Q40, which is way too high.

When I map the long reads to a genome and use cramino or googles' Best to judge mapped QC from BAM, I get a distribution around Q20 as expected.

Has anyone else seen this ? I'll try to find a public dataset ONT R10.4.1 I can see this in from the last few months.

Thanks,
Colin

[colindaven]

To check this you can get a public fungal dataset here (120mb fastq.gz) and try it.

https://www.ebi.ac.uk/ena/browser/view/ERR14662109

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR146/009/ERR14662109/ERR14662109.fastq.gz

I get 86% of bases over Q30 after standard filtering, which is way, way too high.

Q20 bases:	117.891109 M (95.808300%)
Q30 bases:	106.061401 M (86.194477%)
Q40 bases:	16.867597 M (13.708038%)

Importantly, we also see this problem with the tool fastp, not just fastplong, even though it is not so granular with its Qscore output bins, and it does not provide graphics.

Thanks for looking at this.

Colin

[colindaven]

There is also another public file, here, which was likely called with methylation headers ( abstract here https://www.ncbi.nlm.nih.gov/sra/?term=ERR14126976 )

Download - 25 mb.

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR141/076/ERR14126976/ERR14126976.fastq.gz

QV seen in fastplong - peak at QV50

fastplong -i ERR14126976.fastq.gz 
Trying to detect adapter sequence at read start
Detected: TTTAACCGTGTTACCGTGGGAATGAATCCTTGTTTTCGCATTTATCGTGAAACGCTTTCGCGTT
Trying to detect adapter sequence at read end
Detected: GTTACCGTGGGAATGAATCCTTGTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTC

Before filtering:
total reads: 7930
total bases: 34337668
Q20 bases: 31990855(93.1655%)
Q30 bases: 30040552(87.4857%)

After filtering:
total reads: 7700
total bases: 34223032
Q20 bases: 31943656(93.3396%)
Q30 bases: 30006871(87.6803%)

Filtering result:
reads passed filter: 7700
reads failed due to low quality: 217
reads failed due to too many N: 0
reads failed due to too short: 43
reads with adapter trimmed: 802
bases trimmed due to adapters: 96277

JSON report: fastplong.json
HTML report: fastplong.html

fastplong -i ERR14126976.fastq.gz 
fastplong v0.2.2, time used: 2 seconds

[s-j-mc]

I could be wrong, but I think this is due to how fastplong is calculating quality scores (#10). I think fastplong is taking an average/median of the quality scores, as opposed to taking an average/median of the error rate and then converting this to a quality score. The latter will lead the the Qscores in the range we usually expect for ONT data. How does your data look when you using something like nanoplot?

[colindaven]

Hi @s-j-mc thanks for your suggestion. I think you may very well be correct since the data works as expected with Nanoplot etc.

Checking the code, it seems qualities are just added to a total and the average is taken.

See

fastplong/src/read.cpp

Line 136 in aad10f3

total += mQuality->length();

and

fastplong/src/filter.cpp

Line 44 in aad10f3

else if(mOptions->qualfilter.avgQualReq > 0 && (totalQual / rlen)<mOptions->qualfilter.avgQualReq)

@sfchen I think this is a bug and needs to be changed. What do you think ?

We should keep in mind that modern variants of the dorado basecaller (0.9.1 and 0.9.5) are starting to output a qs tag, eg qs:22.6324. This is only for ONT nanopore reads so would not be a general solution for all reads (pacbio, older nanopore etc, nanopore from public repos like SRR etc).

[SwiftSeal]

I'm also seeing this - attached is a fastplong summary of bacterial data basecalled with dna_r10.4.1_e8.2_400bps_sup@v5.2.0_5mC_5hmC@v1, where the majority of reads have been classified as Q50. Also seen similar in plant reads.