While mean error is properly calculated within reads, is it properly calculated across reads?

Question

While mean error is properly calculated within reads, is it properly calculated across reads?

schorlton opened this issue 2 years ago · 9 comments

As per title.

cat test.fastq

@test_1
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
0,-/(/--/133451659;6875662/''$.110%&)(./1137540.+((&(##%&$)***,-30..*-,/20/0.,/124466640(-.',+)*,/.-,,+-&$$#%*(+,+/2/%'((&&$*+'*++'(()$%+-..1374,--'-&$$$$'.0/00340.--*++0/003513.0&&'-00,./03240/-,-.,,014'%%%0211+04753./,+1(*/222.&-,(()$$%$$/32300./2&,0%)00/6//00-+%&%&&'&/+/22321+-,'%
@test_2
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
AA=96802,2-:9539;9897,+56771<:;<<;;4:;8008.<>;92@BBBA?(&54<=;4233334/$&&+2257./7421/4;B?;A5699@<:@:-)8145.335000

Manual implementation:

#! /usr/bin/env python

from statistics import mean
from Bio.SeqIO.QualityIO import FastqGeneralIterator
import math


def decode(c):
    return ord(c) - 33

def phred_to_probability(phred_score):
    return 10**(-phred_score/10)

def probability_to_phred(error_probability):
    return -10 * math.log10(error_probability)

def raw_q_to_probability(q):
    return phred_to_probability(decode(q))


def calc_stats(fastq):
    mean_read_errors = []
    mean_read_quals = []
    with open(fastq) as in_handle:
        for _title, _seq, qual in FastqGeneralIterator(in_handle):

            error_probabilities = list(map(raw_q_to_probability, qual))
            mean_read_error = mean(error_probabilities)
            mean_read_errors.append(mean_read_error)
            mean_read_quals.append(probability_to_phred(mean_read_error))
    
    print(f"Averaging errors: {probability_to_phred(mean(mean_read_errors))}")
    print(f"Averaging quals: {mean(mean_read_quals)}")

calc_stats("test.fastq")

(base) mambauser@1:/phred_calc$ python3 mean_across_reads.py 
Averaging errors: 11.13979284230169
Averaging quals: 12.12473614554866

NanoStat:

(base) mambauser@1:/phred_calc$ NanoStat --fastq test.fastq -n nanostat.txt
(base) mambauser@1:/phred_calc$ cat nanostat.txt 
General summary:         
Mean read length:            198.0
Mean read quality:            12.1
Median read length:          198.0
Median read quality:          12.1
Number of reads:               2.0
Read length N50:             284.0
STDEV read length:            86.0
Total bases:                 396.0
Number, percentage and megabases of reads above quality cutoffs
>Q5:	2 (100.0%) 0.0Mb
>Q7:	2 (100.0%) 0.0Mb
>Q10:	1 (50.0%) 0.0Mb
>Q12:	1 (50.0%) 0.0Mb
>Q15:	1 (50.0%) 0.0Mb
Top 5 highest mean basecall quality scores and their read lengths
1:	15.2 (112)
2:	9.1 (284)
3:	NA
4:	NA
5:	NA
Top 5 longest reads and their mean basecall quality score
1:	284 (9.1)
2:	112 (15.2)
3:	NA
4:	NA
5:	NA

Thank you!!

Answer 1 · 2023-02-09T19:25:09.000Z

Hi,

I think this is more a comment than a question, as you have already demonstrated that it is only calculated adequately within reads, and not across reads. Is this then a suggestion to change things?

Cheers,
Wouter

Answer 2 · 2023-02-09T22:27:08.000Z

If possible? I recognize you're potentially working on cramino as a successor but as far as I can tell, that is not compatible with FASTQs, and this is a moderately confusing bug given the previous issues and blog post on how to correctly average Phred scores. Thank you for all of your tools and effort!

Answer 3 · 2023-03-21T14:51:13.000Z

It has taken a while, but nanomath v1.3.0 now has the requested change.

Cheers,
Wouter

Answer 4 · 2023-03-21T16:11:24.000Z

Thank you!! Really appreciate it.

Out of curiosity, how do you push updates to PyPI/Bioconda? I don't see it tagged in GitHub, so do you have to trigger a manual push to PyPI? I assume Bioconda tracks the PyPI version but don't see it pushed yet.

Thanks again!

Answer 5 · 2023-03-21T20:18:27.000Z

Ah, pushing to PyPI is on the command line, using twine (https://gigabaseorgigabyte.wordpress.com/2017/08/03/a-bash-alias-for-easier-uploading-a-project-to-pypi/). Usually, the biocondabot picks up that PyPI is updated, and then opens a PR on bioconda with the change. Let's give it some time :)
I will do it manually if necessary, and feel free to remind me if it looks like I'm forgetting this.

Answer 6 · 2023-04-03T02:36:33.000Z

This has made its way into bioconda today. Thank you again for the fix - closing!

Answer 7 · 2023-04-29T05:01:16.000Z

Hi @wdecoster , not sure if this is actually fixed? I'm looking at the median. My thinking is median read quality should be a median of medians. Does this make sense? See python logic.

docker run --rm -it -v $(pwd):$(pwd) -w $(pwd) mambaorg/micromamba
micromamba install -y nanostat 'nanomath>=1.3.0' -c conda-forge -c bioconda -c defaults

cat test.fastq

@03c90146-89b4-4e0f-aacb-29f6434c78a9
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
FDFGFEGHGHFDB@@AA72@/IHCBBBBEDEE??RJFGC@@;;;=GHDED<829--89HHIKKFGFHIKIIFFFFHN00000CFFJMJ@;JDEAAABCCBBBDFJHIHGEEEGLN{PLDDCCDHHDEEEOKEB752/-,+++++,./.----..%%&&&'&&&''(--544IHMGAAA@AEGGGPNHHGF@6=/,,,,2/))&%%%%%&)'('''(-78EG???=;=:=><8==;?><96/,)(&&'*'%$%'&%$$&%'&''++*****))(.-.022346709:;<1)((,1155?B@@@:-,++****+++++,9556;=210*&&&&&('(()*0///0+))),=?=:;443207-4555:732222:;;:?@:57<<A?>=,,,,1**'''()+*,--0:8763-,,//8AAA?A=(''''''('''''9A977*)))+77767:;>AA??B7666685-././;;753-+,-0111.--,(*56730*)+((((,*06-76443532...+),+++,6978,+(('''''2005.34481.***//.**()*&%%'+:=8543324.29IHB@?<32234=11-)$##$&'0=>=A852200000ACMRROQ?89('((()(&'200.--*)))('''')'&'''%&&%%&&%&(')**))('&&''&
@391f2d96-ea70-4665-8c4a-21b185c41525
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
3===RONT{{QKLIFGHEH=::::CAAA@BBA@CCC....00002::;IBCHKJIIHHHEGGGHHGIJGFFFFMKKPILMLGGGGGIMNMEB@CBCDEKNMQQSJJMT{UKOKLLEEDEELGKIJKPEBA==>??@@@HD@?@AH?>?>?IIKIPFFEEFFKIEDDCEHFFHIHIKMLPOHQOH
@7c333852-7b92-44d0-a532-f241c478ed97
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
"#%$$$$$$')++**'%%%%''&%$$$$&+++-.54542/-)&&&$$$$%&%'0-)&'''(&%%&&,,..,+--))+***+.*+)('&$$%&'(&$$$$%%$'))746?CEFDB=<;BDCCDCB@BAACCCDDEFHHJJHGGHJR{KIHHBBBBEGHLLEQHFCSLNQNOPCBAAAGFEFEFGEGIFLLJJLLQKLJIJOIIDBAAAAHHHEGQM<7640/-.?@EJ?;;<;IKKHF<9763)(&')****++BAA>>:7/4.+)&&'''(,*%%%%%57AC;AABB33333<<=C-,,,,IHGIFFFFELIKHIGHC@A@ADDDDEGHJHGGIIPEMMHGHFR{NJNHFGFEFFHIGHIJHIKHIKPQKLIRKOPOLNQU{IIGGFFGFFEFDEFDCEDEEEGFGFFDBIJD>====C000030/2('''1DEFMJGCDDEEDDFEKPPOKIIK<;;;;GHGGKQ{MF@6455?EDDFHJLIIGGHH32226BCDIJNLDJA=<<=?DFLI@=;@DJHHGGDB42*((),)**,,,,,;>CDGBCKRKMGEBA<<3,+.:C><<<<DKFHEEDDDEG;321113;NPTOLNLOMILJIKLIIJIHIJHJJIHGGIIJOKLLROFDEEDHHDCDDDHGDGEGHHCC@@960,,+&'&''(:;?GHHCILMKKHGHECCDDEJJIHEDBB,+++**&$%%+78=BBCCDDIIHJJKJGFGDEEDEEFIIKQLRL@@@@>;999:IHLMFDDCDD@G{KJGFFEDFHFEEDEFACBEFEHOIAAFGFDCCCECEEEFFEEC:956?E
@f52f5b20-8ca0-4d70-b034-7fdc9e4a4781
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
===>?HIGJMXOKLIGFFIMDDDCDIHH76666666JKKJKKJKIKLRLIKILLQBN{KJIHKLO??>;;;<=AAAA@@HJGGBAAABIEKKJGHHGKOIKMJKOMIKLKHIHFJJKJLJKGIHKKIGFEEE***=@@A{KMNNQQNKIQ{KQKQNLIIJKNOKJKHIIHHGFJHGIGFDDDAB@DDECDCCB<:1--2223CDGFHHLNNKLJJKLLNJHJJNGGECDDEJKIIIFEJOM{8881****75:OIGSJP{{NLP{UMJL{TQ{SMJHIJIDEDCEOIMOR{H:::::7,(''(.//222;=CLIKKHGGHJIIHFEFGEGIHGEDDDFHIKH
@e8b38fb1-f33e-4f85-9eb6-6f2544ac1a88
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
%$&&&&%%(((*4++))'%$$$(.+*('(&%%%''(',******%%%%&&'%%$###$##$$$%%$$%%$$$$'))%%%$$$&%%(((/((((((%$$$%%%%%$%)..,++(((%&(&'))++)%%%%'&&%%)'&%%%%%&$%&&'',*%%%$$)'((%%&&%'&%'()(%*('&&'%(&&)&&$##$$#$$##&%'&&(''''''&'&%%%%$%'*'&&%%&%&&&&$$&&&))()'''(&%$###%&'&&$%&&(,,86556779;?>AJGECDDHLLPLKHGEFEGIGKBAAABHHGHINII===:)('$$$&(('''''.48??DBA:5666683333?BDDCBDCBBBCEBBABB)(77:29=<<=A:54421444HGFEEDFECCBCCCCCGJGGJ?B0///4899BBAABCBBD3221000()*??=AEG????AF<;:99=<@@@CDFGHGC@??>ABFHK;:;;:::;666CCDCFGFEEFHHRMIFGFDBA@>>>??BEEEGIMC????FFHGDCCCBB=;;;:<<=OPKJGCDDEKHYLC>>==?BCEDBBA?ADFEADDHH>;;;;@BABCEDB@??=<<89<?@?==<<<:9850568?FFHFEFEHGFDDDDDCAA@=0+--/0=<;;;<KLMIKDDINQQKHDCDABCDLJIFECDDC7JICBBBBHFHEED****+0.&&'+../1122225BBDDCDCCDDBBD,++*)(()?;76:>44225;<CAA@@AD32222221600334=DFGFFGFAI9;;;:997;;:=)(((1/-----/DCCCEEEDEIGIKHGFFIL{GJNNLOKJE?989?<ABHCAAAAFGHGFJHIEGGGFCFHDIHKGFFEHG733@@?/......----99:;AADQQLJKJHJJKJIKNLOMNGNNLL;;;;;NJI+)))*.,+,-110+-0.0001BCCBCKPJKKP{JKH---,,++129;;;FGLJGIED<;..*(''',-=ACBGKIJUVWNOKNLNKJKJJJIHIM5544<;977767;KLNIFEGGBBCBCKID:99844>CLJMKLNL{LONQPKTMLKOLRMLLKJLNJLMNN9766:CEKHIGHJKMIJD97,,,,)&&%(.()../@DEDFEH>A::;HDC...**+,-8BKJmRILKNSHIRQSRQMKIJLIQOMQHLGIJNA7,,GJIKHCQUPNOGMJOJKIIL?@@@MUNMIVKIIO{{QKOORL{POMMQMMGCBCCEKJJQKNWN{QMQKNMMOJNMAAB@?2110000/5557<;;433*+AFJQMKJHGFDD@;;3.....'$$$$&####$)()(''',=ACCAA?899:IJMOMPJIAAA@A<65311666A=56ECDBA@><<<<>A?@GFHFB

NanoStat --fastq test.fastq 
General summary:         
Mean read length:             677.4
Mean read quality:             12.6
Median read length:           674.0
Median read quality:           12.7
Number of reads:                5.0
Read length N50:              805.0
STDEV read length:            416.8
Total bases:                3,387.0
Number, percentage and megabases of reads above quality cutoffs
>Q5:	5 (100.0%) 0.0Mb
>Q7:	5 (100.0%) 0.0Mb
>Q10:	5 (100.0%) 0.0Mb
>Q12:	3 (60.0%) 0.0Mb
>Q15:	2 (40.0%) 0.0Mb
Top 5 highest mean basecall quality scores and their read lengths
1:	26.3 (184)
2:	21.4 (342)
3:	12.7 (805)
4:	11.0 (1382)
5:	10.8 (674)
Top 5 longest reads and their mean basecall quality score
1:	1382 (11.0)
2:	805 (12.7)
3:	674 (10.8)
4:	342 (21.4)
5:	184 (26.3

#! /usr/bin/env python

from statistics import median
from Bio.SeqIO.QualityIO import FastqGeneralIterator
import math


def decode(c):
    return ord(c) - 33

def phred_to_probability(phred_score):
    return 10**(-phred_score/10)

def probability_to_phred(error_probability):
    return -10 * math.log10(error_probability)

def raw_q_to_probability(q):
    return phred_to_probability(decode(q))


def calc_stats(fastq):
    median_read_errors = []
    median_read_quals = []
    with open(fastq) as in_handle:
        for _title, _seq, qual in FastqGeneralIterator(in_handle):

            error_probabilities = list(map(raw_q_to_probability, qual))
            median_read_error = median(error_probabilities)
            median_read_errors.append(median_read_error)
    
    print(f"Median of median read errors: {probability_to_phred(median(median_read_errors))}")

calc_stats("test.fastq")

python calc_stats.py
Median of median read errors: 35.0

Answer 8 · 2023-04-29T06:08:45.000Z

Hmmm, I expect it to do what you expect, but I have to confirm. Median Q of 35 seems unreasonably high, though.

Answer 9 · 2023-05-01T13:22:14.000Z

You are correct. Taking median of the base error probabilities does not reflect the overall error probability for the read; average is more appropriate. I believe NanoStat performs this correctly and closing with an example:

Thanks for your help!