4
$\begingroup$

In the fastq format every 4k+4-th line contains the positionwise qualityscore (ascii encoded):

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 
...

which program uses this information besides trimming tools (cutadapt) or quality control (fastqc)?

I heard that segemehl (mapping) abandons this quality string in favor of running time (maybe other mapping tools too?) But what about abyss, trinity, ...?

UPDATE:

So i atrifically generated a fastq file containing 2 identical reads that perfectly match a sequence in some genome (C.fna). The two reads have the most extreme phred scores (only '!!'=highest sanger phread and only 'II'=lowest sanger phred)

Now using segemehl the mapping quality is identical (60)

@HD VN:1.0
@SQ SN:C_fna    LN:15440
@RG ID:A1   SM:sample1  LB:library1 PU:unit1    PL:illumina
@PG ID:segemehl VN:0.3.4    CL:segemehl -i C2.fna.idx -d C2.fna -q test.fq
A   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
B   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
[SEGEMEHL] Thu Feb  6 08:35:15 2020: matching w/ suffixarray has taken 0.000000 seconds.

So segemehl does NOT use the read quality phred (basically only using the read sequence for mapping)

Improve this question
$\endgroup$
1
  • $\begingroup$ Assemblers will generally ignore base quality scores since it's assumed that the input has been heavily QCed already. $\endgroup$
    – Devon Ryan
    Feb 6, 2020 at 7:52

2 Answers 2

Reset to default
4
$\begingroup$

I think almost all alignment programs would use base quality as a factor to calculate mapping quality (See here). These would also be used in Variant Calling to pick reads where high quality bases confidently align as either the reference base or one of the alternate bases.

Improve this answer
$\endgroup$
3
  • $\begingroup$ see update above: segemehl does not use phred $\endgroup$
    – Paul
    Feb 6, 2020 at 7:49
  • $\begingroup$ @Paul Some aligners use phred (e.g., bowtie2), others largely ignore it. Phred scores have changed a bit over time and these days there are only 5 possible values, so their utility has decreased as sequencers have gotten more accurate. $\endgroup$
    – Devon Ryan
    Feb 6, 2020 at 7:50
  • $\begingroup$ But it is kinda weird that in this benchmark (ecseq.com/support/benchmark.html) segemehl is more accurate than bowtie2,bwa,... but slower and more memory hungry. If bowtie2 uses phred, why does it not gain a significant advantage over segemehl? $\endgroup$
    – Paul
    Feb 6, 2020 at 7:56
1
$\begingroup$

LoFreq is a variant calling tool that specifically leverages per-base quality scores.

For long read technologies (PacBio SMRT, ONT), I think quality score usage is mostly limited to qc steps as well, such as filtering out poor reads either by mean quality score of the whole read (as in ONT's basecalling softwares and Filtlong) or by the lowest quality score within a window in a read (Filtlong). I think most long read mappers ignore these values (for example minimap2 and ngmlr do and can take fasta for read inputs).

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.