e-fuhrmann commented 1 month ago

Issue Report

Please describe the issue:

Hey there, I am currently working on some analyses in which I am particularly interested in the basecalling qscores (qs:f). However, I have a problem with reproducing the mean basecall qscores dorado generates. Whenever I calculate the mean basecall qscore manually for any given read (from bam-column QUAL), the result varies (sometimes drastically) from the qs:f-value dorado outputs. I don't quite understand how dorado arrives at these values.

This is based on my understanding (please correct me if I'm wrong here) that

the QUAL String encodes the Phred QS in ASCII (base 33)
the 'mean basecall qscore' (qs:f) is the arithmetic mean of all individual qscores in QUAL in a given read (corresponding to the given nucleotides in SEQ)

Steps to reproduce the issue:

Use dorado basecaller

choose any read and calculate the mean qscore from QUAL. Here a short piece of R-code I used to do so:


# R-env
qual_test <- "[Insert QUAL String here]" # add backslashes before quotations marks first
qual_test <- str_split_1(qual_test, "")
qual_num_test <- list()
for (i in 1:length(qual_test )) {
qual_num_test <- append(qual_num_test, (as.integer(charToRaw(qual_test[i])) - 33))
}
qual_num_test <- unlist(qual_num_test)

mean

mean(qual_num_test)

median

median(qual_num_test)

mode

uniq <- unique(qual_num_test) tab <- tabulate(match(qual_num_test, uniq)) uniq[tab == max(tab)]



## Run environment:
- Dorado version: 0.6.0 (but also seems to happen in 0.5.3, 0.7.0)
- Dorado command: basecaller (with/without --trim adapters, --estimate-poly-a); rna004_130bps_hac@v3.0.1 (but also seems to happen with rna004_130bps_sup@v3.0.1)
- Operating system: Linux
- Hardware: NVIDIA A100 SXM4, 40GB
- Source data type: pod5
- Source data location:  Networked drive
- Details about data:
   - MinION Mk1B (MIN-101B)
   - FLO-MIN004RA
   - seems to work on different datasets, independet of read lengths / number of reads / total dataset size

## Example
`/path/to/0.6.0/dorado basecaller --trim adapters -r /path/to/rna004_130bps_hac@v3.0.1 /path/to/pod5-files > /path/to/out.bam`

Some reads from above output and manual calculation as example:

`d21cfef9-9882-4029-bff0-6c2eb37769ca    4       *       0       0       *       *       0       0       GCAATTCCCAGTATGACCAGGTGTTAGATGAGGCCATCAAAGCTTGCAAAACCACGCAGGTGAACAACAAGGGCATCCAGATCATCTACACCCGCAATCGAGGTGAAGAGTGAG &&'.3;>@@;<@25B@DAB910***69,,+,8;?BAB=4BA.'.<>E?DAB<-*),//87557;???<887<243:9:9=>ADACB;;9-*.;?>6/11::<<<A<70'&&&&% qs:i:13 du:f:0.74125    ns:i:2965       ts:i:0  mx:i:1  ch:i:173   st:Z:2024-05-21T12:33:48.568+00:00      rn:i:-1 fn:Z:FAY52980_pass_020e2a4b_6b8e5b7d_0.pod5     sm:f:98.6667    sd:f:22.2867       sv:Z:quantile   dx:i:0  RG:Z:6b8e5b7dc343c35993f0c99c718711911603f68a_rna004_130bps_hac@v3.0.1  pi:Z:4d59b25b-47e7-4aae-b380-48a69ad4f386  
sp:i:21800`
dorado mean qs: **13**
manual mean/median/mode: **21.89**, 24, 27

`f6ad4e24-a99a-4e98-ae96-528d52b81fdc    4       *       0       0       *       *       0       0       CCTGTACTGAAAATGTGTAAATACCTTGTAACCTTGTCCCATGCGTCATGTACATAGATATATGTCATATCTATATTTATATCCATAGACAGTCAAGGATTTTAATGATCTAAATCATTTGCCACTGGTGCTATGGGAGTAGT    0.,++,/0??EE93'%&&'0/,(''(((***(%$$$334842.,))*()(&%$$-6?>;6;6+')(/,)('&''&')((('%%%&)+??<=>=DDDCA66334ABF,,''((''+()(('%#$&'(&'((()015AB-->>/5       qs:i:8  du:f:3.03625    ns:i:12145      ts:i:8600       mx:i:4  ch:i:424        st:Z:2024-05-22T01:45:00.377+00:00      rn:i:23630 fn:Z:FAY52980_pass_020e2a4b_6b8e5b7d_137.pod5   sm:f:79.5601    sd:f:16.2484    sv:Z:quantile   dx:i:0  RG:Z:6b8e5b7dc343c35993f0c99c718711911603f68a_rna004_130bps_hac@v3.0.1`
dorado mean qs: **8**
manual mean/median/mode: **13.52**, 9, 7

`910b3f2e-3ed5-4c9a-ae35-21f5faaa028e    4       *       0       0       *       *       0       0       CAGCCGGGATTTGTTCGGTTTTGTTTGTGGATCGCTGTGATCGCTCTTGACAATGCAGATCTTCGTGAAGACTCTGACTGGTAAGACCATTCTCCGAGGTTGAGCCCAGTGACACCATCGAGAATGTCAAGGCAAAGATCCAAGATAAGGAAGGCATCCCTCCTGACCAGCAGAGGCTGATCTTTGCTGGAAACAGCTGGAAGATGGGCGCACCCTGTCTGACTACAACATCCAGAAAGAGTCCACCTGCACCTGGTGCTCCGTCTCAGAGGTGGGATGCAAATCTTCGTGAAGACACTCTGGCAAGACCATCACCCTTGAGGTGGAGCCCAGTGACACCATCGAGAACGTCAAAGCAAAGATCCAGGACAAGGAAGGCATTCCTGACCAGCAGAGGTTGATCTTTGCCGGAAAGCAGCTAGAAGATGGGCGCACCTGTCTGACTACAACATCCAGAAAGAGTCTACCCTGCACCTGGTGCTCCGTCTCAGAGGTGGGATGCAGATCTTCGTGAAGACCCTGACTGGTAAGACCATCACCCTCGAGGTGGAGCCCAGTGACACCATCGAGAATGTCAAGGCAGATCCAAGATAAGGAAGGCATTCCTGATCAGCAGAGGTTGATCTTTGCCGGAAAACAGCTGGAAGATGGTCGTACCTGTCTGACTACAACATCCAGAAAGAGTCCTTCTCTGCACCTGGTACTCCGTCTCAGAGGTGGGATGCAAATCTTCGTGAAGACACTCTGCAAGACCATCACCCTTGAGGTCGAGCCCAGTGACACTATCGAGAACGTCAAAGCAAAGATCCAAGACAAGGAAGGCATTCCTGACCAGCAGAGGTTGATCTTTGCCGGAAAGCAGCTGGAAGATGGGCGCACCCTGTTCTGACTACAACATCCAGAAGGAGAGTCTACCCTGCACCTGGTGCTCCGTCTCAGAGGTGGGATGCAGATCTTCGTGAAGACCCTGACTGGTAAGACCATCACCCGAAGTGGAGCCGAGTGACACCATTGAGAATGTCAAGGCAAAGATCCAAGACAAGGAAGGCATCCCTCCTGACCAGCAGAGGTTGATCTTTGCCGGAAACAGCTGGAAGATGGTCGTACCCTGTCTGACTACAACATCCAGAAAGAGTCCTCTCTGCACCTGGTGCTCCGTCTCAGAGGTGGGATGCAGATCTTCGTGAAGACCCTGACTGGTAAGACCATCACTTCGAGGTGGAGCCGAGTGACACCATAGAGAATGTCAAGGCAAAGATCCAAGACAAGGAAGGCATCCCTCCTGACCAGCAGAGGTTGATTTGCTGGGAAACAGCTGGAAGATGGACGCCCTGTCTGACTACAACATCCAGAAAGAGTCCACCCTGCACCTGGTGCTCCGTCTTAGAGGTGGGATGCAGATCTTCGTGAAGACCCTGACTGGTAAGACCATCACACTGAGGAGCCGAGTGTTCCGCTGAGAATGTCAAGGCAAAGATCCAAGACAAGGAAGGCATCCCTCCTGACCAGCAGAGGTTGATCTTTGCTGGGAAACAGCTGGAAGATGGACGCACCCTGTCTGACTACAACATCCAGAAAGAGTCCACCCTGCACCTGGTGCTCCGTCTTAGAGGTGGGATGCAGATCTTCGTGAAGACCCTGACTGGTAAGACCATCACTCAAAGTGGAGCCGAGTGACCCATTGAGAATGTCAAGGCAAAGATCCAAGACAAGGAAGGCATCCCTCCTGACCAGCAGAGGTTGATCTTTGCTGGGAAACAGCTGGAAGATGGACGCACCCTGTCTGACTACAACCAGAAAGAGTCCACCCTGCACCTGGTGCTCCGTCTCAGAGGTGGGATGCAGATCTTCGTGAAGACCCTGACTGGTAAGACCATCACCCTCGAGGTGGAGCCCAGTGACACCATCGAGAATGTCAAGGCAAAGATCCAAGATAAGGAAGGCATCCCTCCTGATCAGCAGAGGTTGATCTTTGCTGGGAAACAGCTGGAAGATGGACGCACCCTGTCTGACTACAACATCCAGAAAGAGTCCACTTCTGCACTTTGGTCCTGCGCTTGAGGGGTGTCTGGTTTCCCCTTTTAAGGTTTCAACAAATTCATTGCACTTTCCTCTTCAATCAGTTGTTGCATTCCCAAAAAAAAAAAAAAAATCC        /1+++,;-,,()(,/37,,.3<=@BB;992000469=951**)*549?>:9:998?B=;92;;<=EFDC>;;<AB3>??56=:8769-*)%'(*068:>>EGOKKKBLHG@C?>@AHLLSF@?@BINSJEEHFGDFEJI?==?===?EEFEFHGEA..759@FFHIGNLLKJDFFLK3DBB:66457:B=;;=EE9885247:CFHFEFHKLOQJID;:ED>::876,+''(=;A:ACA6.))*,3<899453+,,-7;<<<66662*))*-/3<@DB@=;:,324@B<<@@?5122--97BFE;CDEDBA>DGNFFEMLBAACFKKIIHEHKK8>EGA;=>?IM?;999@1.<869;7DGBI=DFGILOC>DGFC91--4AGLKGMJRPQFJGD<9:=51/,)-?ABF?45EIKK?''''CHJ<<;<FC:95.>A@AD:KHM5EJLOMLJHGIAMKIGH@;112/2CDFFGDC?:;<A<3+%$&,391--/9::9:>EF9BEBDB?:;89BD@>>>CDHKG2/+(*.9=AB;LKDH<720159:;DGMSGSSSCB@<9<;B>5SSSLMDECDEEHA<3+*<*+1/)'',>@B4IHCC87/.+AFD9668F;7445)((,;>'32124../<<:9=?BB9776=>C8EE/)))11)9@@ABD@@?=DOPSHFGP@KLHGGHEF<33+.--43++++????D=<:6668E></**+57687?CEK99A<+++786:8:5&%&%$&'&;;:;?A>KFMHE9629CGFEC<>?PNMGJ?<;971891112;?>:833=>8>2/487;A@C@EEGIIMEFEGFBDC>/4?BDB>CBB@A?<3<<<==20,,+,989<?A7>8889?@?ABBEIKF@;99<?8:CA)'(22.01381013438:?=+0'((5;>@CDHDJKIE@A@CB>AB=89:>:..//0CDDAEBC?>C?ABG44339;;;<IIFG?;CA@@==98:<DFLL>@@JMD<;;<((&&%'/+***+9;:@>A.+++@KPM<<:9:@A@::;AAGB>=<ACD@@ABC105/8687699:?GNMIKKOMJD@<///0CE@;69:50007;@@AGF?64456;9GE<3336@?>>DCA@@@CDDHPNQMQMLBA<9:AFFKLC4356;;A>>>>=9::;:93+*)*3=@;:99:766?EJED9=<==>=??DEBIIB>GFJ?;:2''(*;<;60,++*+-8?BD>=5:<<CDINGFA???D1<<==DLF7:88654*8;?F?EE>?;767:INILLGI?::DCB@<<:AG<<<<KCCCC@;<<=J4343.---..++.0;><:<=<9<<:7-))(*<>@<==>E@>>@<C?@AA=@C>J2AA@A@DDDFHBBAA@<BCC?@?=?>BDA==>FEC??@AFLHSOPCF65)**+++//('',9??IL>=:998?BDG97.3*)*)'(()7<=6679=&&%%'&&')BEEB@60.034>=::==BA?<600099+.)(/20:;<923668CCHI>JMLLMJD>=;<9;72-+.7765676589BGEBA?>:50148AAAFPHGGDD?>7D@?>>0+.FDAC>?F@=@IKI<<<;58=<<9:8:535EIG<53-+(%&.<;<;??@>>BEIKDGF=<;>?@DFHJGCCDFLG;72113=?><7,))('&%#-8<===@?2114>97+,*4435EKPKA@328CD@???>?AB@E@EFGHIID3('(&()9?JH==@CCIBCB=;9867:767=4446/>??C@A<;:CACB@@>AADB>BDPDMMIRNQIJD//.-))(-)(&+9A=9:89<JCEHCABE>;;;7<?8FJL>DA@9.,*+,1?A@DE:;<LLOJE:::;8767@DABDC=AEKL8-))+/3>EGCHFGF<9**+,@@?HHHNKMNNGIHE?9></5INNRRIIHFDB;8753587==9/.,45<LNEEJNIO<;53/1++3;?.**-A==>?I;::..-,0+)++1668<76568;8898BBA@CF=;:9;AAB;=;;;7FLSOHGQSKO>DC@A?>?BCD=@<3,-./.-+++,-05598><,++&&'99623463(&/02:766:8350227:;<2/&&&&(*+(&*...5;>>A:966''*+11,)')/2+;>@BC>9840/19ADDC@=<;954,#,---00     qs:i:16 du:f:15.4008    ns:i:61603 ts:i:2150       mx:i:2  ch:i:337        st:Z:2024-05-21T13:14:32.272+00:00      rn:i:1219       fn:Z:FAY52980_pass_020e2a4b_6b8e5b7d_7.pod5        sm:f:88.0974    sd:f:21.0791    sv:Z:quantile   dx:i:0  RG:Z:6b8e5b7dc343c35993f0c99c718711911603f68a_rna004_130bps_hac@v3.0.1`
dorado mean qs: **16**
manual mean/median/mode: **26.46**, 27, 27

##
Am I missing something here? Does dorado calculate something different?

Many thanks in advance!

Psy-Fer commented 1 month ago

Hey, I'm not from nanopore, but I faced this problem when first looking at qscores way back

Mean qscore is the average of the probabilities associated with the quality scores of your read. It is calculated with the following equation:

in python that's

import numpy as np

def calculate_qscore(qstring):
    '''
    calculate a qscore from a qstring
    '''
    qs = (np.array(qstring, 'c').view(np.uint8) - 33)
    mean_err = np.exp(qs * (-np.log(10) / 10.)).mean()
    score = -10 * np.log10(max(mean_err, 1e-4))
    return score

If you run this code on your first qstring: &&'.3;>@@;<@25B@DAB910***69,,+,8;?BAB=4BA.'.<>E?DAB<-*),//87557;???<887<243:9:9=>ADACB;;9-*.;?>6/11::<<<A<70'&&&&%

you get 13.47433188066914

Dorado then floors this into an int, so you get 13

I hope that helps James

malton-ont commented 1 month ago

Hi @e-fuhrmann,

@Psy-Fer is correct regarding the calculation. The other thing to note is that dorado may skip some bases when calculating the mean qscore, for example to exclude the adapter for RNA reads. See here and here for previous discussions on this topic.

e-fuhrmann commented 1 month ago

Hey there,

I understand now, that does the trick!

Thanks a lot!

nanoporetech / dorado