-
Notifications
You must be signed in to change notification settings - Fork 1.7k
/
QualityIO.py
2428 lines (2044 loc) · 91.9 KB
/
QualityIO.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Copyright 2009-2020 by Peter Cock. All rights reserved.
# Copyright 2020 by Michael R. Crusoe
#
# This file is part of the Biopython distribution and governed by your
# choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
# Please see the LICENSE file that should have been included as part of this
# package.
"""Bio.SeqIO support for the FASTQ and QUAL file formats.
Note that you are expected to use this code via the Bio.SeqIO interface, as
shown below.
The FASTQ file format is used frequently at the Wellcome Trust Sanger Institute
to bundle a FASTA sequence and its PHRED quality data (integers between 0 and
90). Rather than using a single FASTQ file, often paired FASTA and QUAL files
are used containing the sequence and the quality information separately.
The PHRED software reads DNA sequencing trace files, calls bases, and
assigns a non-negative quality value to each called base using a logged
transformation of the error probability, Q = -10 log10( Pe ), for example::
Pe = 1.0, Q = 0
Pe = 0.1, Q = 10
Pe = 0.01, Q = 20
...
Pe = 0.00000001, Q = 80
Pe = 0.000000001, Q = 90
In typical raw sequence reads, the PHRED quality valuea will be from 0 to 40.
In the QUAL format these quality values are held as space separated text in
a FASTA like file format. In the FASTQ format, each quality values is encoded
with a single ASCI character using chr(Q+33), meaning zero maps to the
character "!" and for example 80 maps to "q". For the Sanger FASTQ standard
the allowed range of PHRED scores is 0 to 93 inclusive. The sequences and
quality are then stored in pairs in a FASTA like format.
Unfortunately there is no official document describing the FASTQ file format,
and worse, several related but different variants exist. For more details,
please read this open access publication::
The Sanger FASTQ file format for sequences with quality scores, and the
Solexa/Illumina FASTQ variants.
P.J.A.Cock (Biopython), C.J.Fields (BioPerl), N.Goto (BioRuby),
M.L.Heuer (BioJava) and P.M. Rice (EMBOSS).
Nucleic Acids Research 2010 38(6):1767-1771
https://doi.org/10.1093/nar/gkp1137
The good news is that Roche 454 sequencers can output files in the QUAL format,
and sensibly they use PHREP style scores like Sanger. Converting a pair of
FASTA and QUAL files into a Sanger style FASTQ file is easy. To extract QUAL
files from a Roche 454 SFF binary file, use the Roche off instrument command
line tool "sffinfo" with the -q or -qual argument. You can extract a matching
FASTA file using the -s or -seq argument instead.
The bad news is that Solexa/Illumina did things differently - they have their
own scoring system AND their own incompatible versions of the FASTQ format.
Solexa/Illumina quality scores use Q = - 10 log10 ( Pe / (1-Pe) ), which can
be negative. PHRED scores and Solexa scores are NOT interchangeable (but a
reasonable mapping can be achieved between them, and they are approximately
equal for higher quality reads).
Confusingly early Solexa pipelines produced a FASTQ like file but using their
own score mapping and an ASCII offset of 64. To make things worse, for the
Solexa/Illumina pipeline 1.3 onwards, they introduced a third variant of the
FASTQ file format, this time using PHRED scores (which is more consistent) but
with an ASCII offset of 64.
i.e. There are at least THREE different and INCOMPATIBLE variants of the FASTQ
file format: The original Sanger PHRED standard, and two from Solexa/Illumina.
The good news is that as of CASAVA version 1.8, Illumina sequencers will
produce FASTQ files using the standard Sanger encoding.
You are expected to use this module via the Bio.SeqIO functions, with the
following format names:
- "qual" means simple quality files using PHRED scores (e.g. from Roche 454)
- "fastq" means Sanger style FASTQ files using PHRED scores and an ASCII
offset of 33 (e.g. from the NCBI Short Read Archive and Illumina 1.8+).
These can potentially hold PHRED scores from 0 to 93.
- "fastq-sanger" is an alias for "fastq".
- "fastq-solexa" means old Solexa (and also very early Illumina) style FASTQ
files, using Solexa scores with an ASCII offset 64. These can hold Solexa
scores from -5 to 62.
- "fastq-illumina" means newer Illumina 1.3 to 1.7 style FASTQ files, using
PHRED scores but with an ASCII offset 64, allowing PHRED scores from 0
to 62.
We could potentially add support for "qual-solexa" meaning QUAL files which
contain Solexa scores, but thus far there isn't any reason to use such files.
For example, consider the following short FASTQ file::
@EAS54_6_R1_2_1_413_324
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
;;;;;;;;;;;9;7;;.7;393333
This contains three reads of length 25. From the read length these were
probably originally from an early Solexa/Illumina sequencer but this file
follows the Sanger FASTQ convention (PHRED style qualities with an ASCII
offset of 33). This means we can parse this file using Bio.SeqIO using
"fastq" as the format name:
>>> from Bio import SeqIO
>>> for record in SeqIO.parse("Quality/example.fastq", "fastq"):
... print("%s %s" % (record.id, record.seq))
EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC
EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA
EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
The qualities are held as a list of integers in each record's annotation:
>>> print(record)
ID: EAS54_6_R1_2_1_443_348
Name: EAS54_6_R1_2_1_443_348
Description: EAS54_6_R1_2_1_443_348
Number of features: 0
Per letter annotation for: phred_quality
Seq('GTTGCTTCTGGCGTGGGTGGGGGGG')
>>> print(record.letter_annotations["phred_quality"])
[26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
You can use the SeqRecord format method to show this in the QUAL format:
>>> print(record.format("qual"))
>EAS54_6_R1_2_1_443_348
26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18
24 18 18 18 18
<BLANKLINE>
Or go back to the FASTQ format, use "fastq" (or "fastq-sanger"):
>>> print(record.format("fastq"))
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
;;;;;;;;;;;9;7;;.7;393333
<BLANKLINE>
Or, using the Illumina 1.3+ FASTQ encoding (PHRED values with an ASCII offset
of 64):
>>> print(record.format("fastq-illumina"))
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
ZZZZZZZZZZZXZVZZMVZRXRRRR
<BLANKLINE>
You can also get Biopython to convert the scores and show a Solexa style
FASTQ file:
>>> print(record.format("fastq-solexa"))
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
ZZZZZZZZZZZXZVZZMVZRXRRRR
<BLANKLINE>
Notice that this is actually the same output as above using "fastq-illumina"
as the format! The reason for this is all these scores are high enough that
the PHRED and Solexa scores are almost equal. The differences become apparent
for poor quality reads. See the functions solexa_quality_from_phred and
phred_quality_from_solexa for more details.
If you wanted to trim your sequences (perhaps to remove low quality regions,
or to remove a primer sequence), try slicing the SeqRecord objects. e.g.
>>> sub_rec = record[5:15]
>>> print(sub_rec)
ID: EAS54_6_R1_2_1_443_348
Name: EAS54_6_R1_2_1_443_348
Description: EAS54_6_R1_2_1_443_348
Number of features: 0
Per letter annotation for: phred_quality
Seq('TTCTGGCGTG')
>>> print(sub_rec.letter_annotations["phred_quality"])
[26, 26, 26, 26, 26, 26, 24, 26, 22, 26]
>>> print(sub_rec.format("fastq"))
@EAS54_6_R1_2_1_443_348
TTCTGGCGTG
+
;;;;;;9;7;
<BLANKLINE>
If you wanted to, you could read in this FASTQ file, and save it as a QUAL file:
>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("Quality/example.fastq", "fastq")
>>> with open("Quality/temp.qual", "w") as out_handle:
... SeqIO.write(record_iterator, out_handle, "qual")
3
You can of course read in a QUAL file, such as the one we just created:
>>> from Bio import SeqIO
>>> for record in SeqIO.parse("Quality/temp.qual", "qual"):
... print("%s read of length %d" % (record.id, len(record.seq)))
EAS54_6_R1_2_1_413_324 read of length 25
EAS54_6_R1_2_1_540_792 read of length 25
EAS54_6_R1_2_1_443_348 read of length 25
Notice that QUAL files don't have a proper sequence present! But the quality
information is there:
>>> print(record)
ID: EAS54_6_R1_2_1_443_348
Name: EAS54_6_R1_2_1_443_348
Description: EAS54_6_R1_2_1_443_348
Number of features: 0
Per letter annotation for: phred_quality
Undefined sequence of length 25
>>> print(record.letter_annotations["phred_quality"])
[26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
Just to keep things tidy, if you are following this example yourself, you can
delete this temporary file now:
>>> import os
>>> os.remove("Quality/temp.qual")
Sometimes you won't have a FASTQ file, but rather just a pair of FASTA and QUAL
files. Because the Bio.SeqIO system is designed for reading single files, you
would have to read the two in separately and then combine the data. However,
since this is such a common thing to want to do, there is a helper iterator
defined in this module that does this for you - PairedFastaQualIterator.
Alternatively, if you have enough RAM to hold all the records in memory at once,
then a simple dictionary approach would work:
>>> from Bio import SeqIO
>>> reads = SeqIO.to_dict(SeqIO.parse("Quality/example.fasta", "fasta"))
>>> for rec in SeqIO.parse("Quality/example.qual", "qual"):
... reads[rec.id].letter_annotations["phred_quality"]=rec.letter_annotations["phred_quality"]
You can then access any record by its key, and get both the sequence and the
quality scores.
>>> print(reads["EAS54_6_R1_2_1_540_792"].format("fastq"))
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
<BLANKLINE>
It is important that you explicitly tell Bio.SeqIO which FASTQ variant you are
using ("fastq" or "fastq-sanger" for the Sanger standard using PHRED values,
"fastq-solexa" for the original Solexa/Illumina variant, or "fastq-illumina"
for the more recent variant), as this cannot be detected reliably
automatically.
To illustrate this problem, let's consider an artificial example:
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> test = SeqRecord(Seq("NACGTACGTA"), id="Test", description="Made up!")
>>> print(test.format("fasta"))
>Test Made up!
NACGTACGTA
<BLANKLINE>
>>> print(test.format("fastq"))
Traceback (most recent call last):
...
ValueError: No suitable quality scores found in letter_annotations of SeqRecord (id=Test).
We created a sample SeqRecord, and can show it in FASTA format - but for QUAL
or FASTQ format we need to provide some quality scores. These are held as a
list of integers (one for each base) in the letter_annotations dictionary:
>>> test.letter_annotations["phred_quality"] = [0, 1, 2, 3, 4, 5, 10, 20, 30, 40]
>>> print(test.format("qual"))
>Test Made up!
0 1 2 3 4 5 10 20 30 40
<BLANKLINE>
>>> print(test.format("fastq"))
@Test Made up!
NACGTACGTA
+
!"#$%&+5?I
<BLANKLINE>
We can check this FASTQ encoding - the first PHRED quality was zero, and this
mapped to a exclamation mark, while the final score was 40 and this mapped to
the letter "I":
>>> ord('!') - 33
0
>>> ord('I') - 33
40
>>> [ord(letter)-33 for letter in '!"#$%&+5?I']
[0, 1, 2, 3, 4, 5, 10, 20, 30, 40]
Similarly, we could produce an Illumina 1.3 to 1.7 style FASTQ file using PHRED
scores with an offset of 64:
>>> print(test.format("fastq-illumina"))
@Test Made up!
NACGTACGTA
+
@ABCDEJT^h
<BLANKLINE>
And we can check this too - the first PHRED score was zero, and this mapped to
"@", while the final score was 40 and this mapped to "h":
>>> ord("@") - 64
0
>>> ord("h") - 64
40
>>> [ord(letter)-64 for letter in "@ABCDEJT^h"]
[0, 1, 2, 3, 4, 5, 10, 20, 30, 40]
Notice how different the standard Sanger FASTQ and the Illumina 1.3 to 1.7 style
FASTQ files look for the same data! Then we have the older Solexa/Illumina
format to consider which encodes Solexa scores instead of PHRED scores.
First let's see what Biopython says if we convert the PHRED scores into Solexa
scores (rounding to one decimal place):
>>> for q in [0, 1, 2, 3, 4, 5, 10, 20, 30, 40]:
... print("PHRED %i maps to Solexa %0.1f" % (q, solexa_quality_from_phred(q)))
PHRED 0 maps to Solexa -5.0
PHRED 1 maps to Solexa -5.0
PHRED 2 maps to Solexa -2.3
PHRED 3 maps to Solexa -0.0
PHRED 4 maps to Solexa 1.8
PHRED 5 maps to Solexa 3.3
PHRED 10 maps to Solexa 9.5
PHRED 20 maps to Solexa 20.0
PHRED 30 maps to Solexa 30.0
PHRED 40 maps to Solexa 40.0
Now here is the record using the old Solexa style FASTQ file:
>>> print(test.format("fastq-solexa"))
@Test Made up!
NACGTACGTA
+
;;>@BCJT^h
<BLANKLINE>
Again, this is using an ASCII offset of 64, so we can check the Solexa scores:
>>> [ord(letter)-64 for letter in ";;>@BCJT^h"]
[-5, -5, -2, 0, 2, 3, 10, 20, 30, 40]
This explains why the last few letters of this FASTQ output matched that using
the Illumina 1.3 to 1.7 format - high quality PHRED scores and Solexa scores
are approximately equal.
"""
import warnings
from math import log
from abc import abstractmethod
from typing import Callable
from typing import IO
from collections.abc import Iterator
from collections.abc import Mapping
from typing import Optional
from collections.abc import Sequence
from typing import Union
from Bio import BiopythonParserWarning
from Bio import BiopythonWarning
from Bio import StreamModeError
from Bio.File import as_handle
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from .Interfaces import _clean
from .Interfaces import _get_seq_string
from .Interfaces import _TextIOSource
from .Interfaces import SequenceIterator
from .Interfaces import SequenceWriter
# define score offsets. See discussion for differences between Sanger and
# Solexa offsets.
SANGER_SCORE_OFFSET = 33
SOLEXA_SCORE_OFFSET = 64
def solexa_quality_from_phred(phred_quality: float) -> float:
"""Convert a PHRED quality (range 0 to about 90) to a Solexa quality.
PHRED and Solexa quality scores are both log transformations of a
probality of error (high score = low probability of error). This function
takes a PHRED score, transforms it back to a probability of error, and
then re-expresses it as a Solexa score. This assumes the error estimates
are equivalent.
How does this work exactly? Well the PHRED quality is minus ten times the
base ten logarithm of the probability of error::
phred_quality = -10*log(error,10)
Therefore, turning this round::
error = 10 ** (- phred_quality / 10)
Now, Solexa qualities use a different log transformation::
solexa_quality = -10*log(error/(1-error),10)
After substitution and a little manipulation we get::
solexa_quality = 10*log(10**(phred_quality/10.0) - 1, 10)
However, real Solexa files use a minimum quality of -5. This does have a
good reason - a random base call would be correct 25% of the time,
and thus have a probability of error of 0.75, which gives 1.25 as the PHRED
quality, or -4.77 as the Solexa quality. Thus (after rounding), a random
nucleotide read would have a PHRED quality of 1, or a Solexa quality of -5.
Taken literally, this logarithic formula would map a PHRED quality of zero
to a Solexa quality of minus infinity. Of course, taken literally, a PHRED
score of zero means a probability of error of one (i.e. the base call is
definitely wrong), which is worse than random! In practice, a PHRED quality
of zero usually means a default value, or perhaps random - and therefore
mapping it to the minimum Solexa score of -5 is reasonable.
In conclusion, we follow EMBOSS, and take this logarithmic formula but also
apply a minimum value of -5.0 for the Solexa quality, and also map a PHRED
quality of zero to -5.0 as well.
Note this function will return a floating point number, it is up to you to
round this to the nearest integer if appropriate. e.g.
>>> print("%0.2f" % round(solexa_quality_from_phred(80), 2))
80.00
>>> print("%0.2f" % round(solexa_quality_from_phred(50), 2))
50.00
>>> print("%0.2f" % round(solexa_quality_from_phred(20), 2))
19.96
>>> print("%0.2f" % round(solexa_quality_from_phred(10), 2))
9.54
>>> print("%0.2f" % round(solexa_quality_from_phred(5), 2))
3.35
>>> print("%0.2f" % round(solexa_quality_from_phred(4), 2))
1.80
>>> print("%0.2f" % round(solexa_quality_from_phred(3), 2))
-0.02
>>> print("%0.2f" % round(solexa_quality_from_phred(2), 2))
-2.33
>>> print("%0.2f" % round(solexa_quality_from_phred(1), 2))
-5.00
>>> print("%0.2f" % round(solexa_quality_from_phred(0), 2))
-5.00
Notice that for high quality reads PHRED and Solexa scores are numerically
equal. The differences are important for poor quality reads, where PHRED
has a minimum of zero but Solexa scores can be negative.
Finally, as a special case where None is used for a "missing value", None
is returned:
>>> print(solexa_quality_from_phred(None))
None
"""
if phred_quality is None:
# Assume None is used as some kind of NULL or NA value; return None
# e.g. Bio.SeqIO gives Ace contig gaps a quality of None.
return None
elif phred_quality > 0:
# Solexa uses a minimum value of -5, which after rounding matches a
# random nucleotide base call.
return max(-5.0, 10 * log(10 ** (phred_quality / 10.0) - 1, 10))
elif phred_quality == 0:
# Special case, map to -5 as discussed in the docstring
return -5.0
else:
raise ValueError(
f"PHRED qualities must be positive (or zero), not {phred_quality!r}"
)
def phred_quality_from_solexa(solexa_quality: float) -> float:
"""Convert a Solexa quality (which can be negative) to a PHRED quality.
PHRED and Solexa quality scores are both log transformations of a
probality of error (high score = low probability of error). This function
takes a Solexa score, transforms it back to a probability of error, and
then re-expresses it as a PHRED score. This assumes the error estimates
are equivalent.
The underlying formulas are given in the documentation for the sister
function solexa_quality_from_phred, in this case the operation is::
phred_quality = 10*log(10**(solexa_quality/10.0) + 1, 10)
This will return a floating point number, it is up to you to round this to
the nearest integer if appropriate. e.g.
>>> print("%0.2f" % round(phred_quality_from_solexa(80), 2))
80.00
>>> print("%0.2f" % round(phred_quality_from_solexa(20), 2))
20.04
>>> print("%0.2f" % round(phred_quality_from_solexa(10), 2))
10.41
>>> print("%0.2f" % round(phred_quality_from_solexa(0), 2))
3.01
>>> print("%0.2f" % round(phred_quality_from_solexa(-5), 2))
1.19
Note that a solexa_quality less then -5 is not expected, will trigger a
warning, but will still be converted as per the logarithmic mapping
(giving a number between 0 and 1.19 back).
As a special case where None is used for a "missing value", None is
returned:
>>> print(phred_quality_from_solexa(None))
None
"""
if solexa_quality is None:
# Assume None is used as some kind of NULL or NA value; return None
return None
if solexa_quality < -5:
warnings.warn(
f"Solexa quality less than -5 passed, {solexa_quality!r}", BiopythonWarning
)
return 10 * log(10 ** (solexa_quality / 10.0) + 1, 10)
def _get_phred_quality(record: SeqRecord) -> Union[list[float], list[int]]:
"""Extract PHRED qualities from a SeqRecord's letter_annotations (PRIVATE).
If there are no PHRED qualities, but there are Solexa qualities, those are
used instead after conversion.
"""
try:
return record.letter_annotations["phred_quality"]
except KeyError:
pass
try:
return [
phred_quality_from_solexa(q)
for q in record.letter_annotations["solexa_quality"]
]
except KeyError:
raise ValueError(
"No suitable quality scores found in "
"letter_annotations of SeqRecord (id=%s)." % record.id
) from None
# Only map 0 to 93, we need to give a warning on truncating at 93
_phred_to_sanger_quality_str = {
qp: chr(min(126, qp + SANGER_SCORE_OFFSET)) for qp in range(93 + 1)
}
# Only map -5 to 93, we need to give a warning on truncating at 93
_solexa_to_sanger_quality_str = {
qs: chr(min(126, int(round(phred_quality_from_solexa(qs)) + SANGER_SCORE_OFFSET)))
for qs in range(-5, 93 + 1)
}
def _get_sanger_quality_str(record: SeqRecord) -> str:
"""Return a Sanger FASTQ encoded quality string (PRIVATE).
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> r = SeqRecord(Seq("ACGTAN"), id="Test",
... letter_annotations = {"phred_quality":[50, 40, 30, 20, 10, 0]})
>>> _get_sanger_quality_str(r)
'SI?5+!'
If as in the above example (or indeed a SeqRecord parser with Bio.SeqIO),
the PHRED qualities are integers, this function is able to use a very fast
pre-cached mapping. However, if they are floats which differ slightly, then
it has to do the appropriate rounding - which is slower:
>>> r2 = SeqRecord(Seq("ACGTAN"), id="Test2",
... letter_annotations = {"phred_quality":[50.0, 40.05, 29.99, 20, 9.55, 0.01]})
>>> _get_sanger_quality_str(r2)
'SI?5+!'
If your scores include a None value, this raises an exception:
>>> r3 = SeqRecord(Seq("ACGTAN"), id="Test3",
... letter_annotations = {"phred_quality":[50, 40, 30, 20, 10, None]})
>>> _get_sanger_quality_str(r3)
Traceback (most recent call last):
...
TypeError: A quality value of None was found
If (strangely) your record has both PHRED and Solexa scores, then the PHRED
scores are used in preference:
>>> r4 = SeqRecord(Seq("ACGTAN"), id="Test4",
... letter_annotations = {"phred_quality":[50, 40, 30, 20, 10, 0],
... "solexa_quality":[-5, -4, 0, None, 0, 40]})
>>> _get_sanger_quality_str(r4)
'SI?5+!'
If there are no PHRED scores, but there are Solexa scores, these are used
instead (after the appropriate conversion):
>>> r5 = SeqRecord(Seq("ACGTAN"), id="Test5",
... letter_annotations = {"solexa_quality":[40, 30, 20, 10, 0, -5]})
>>> _get_sanger_quality_str(r5)
'I?5+$"'
Again, integer Solexa scores can be looked up in a pre-cached mapping making
this very fast. You can still use approximate floating point scores:
>>> r6 = SeqRecord(Seq("ACGTAN"), id="Test6",
... letter_annotations = {"solexa_quality":[40.1, 29.7, 20.01, 10, 0.0, -4.9]})
>>> _get_sanger_quality_str(r6)
'I?5+$"'
Notice that due to the limited range of printable ASCII characters, a
PHRED quality of 93 is the maximum that can be held in an Illumina FASTQ
file (using ASCII 126, the tilde). This function will issue a warning
in this situation.
"""
# TODO - This functions works and is fast, but it is also ugly
# and there is considerable repetition of code for the other
# two FASTQ variants.
try:
# These take priority (in case both Solexa and PHRED scores found)
qualities = record.letter_annotations["phred_quality"]
except KeyError:
# Fall back on solexa scores...
pass
else:
# Try and use the precomputed mapping:
try:
return "".join(_phred_to_sanger_quality_str[qp] for qp in qualities)
except KeyError:
# Could be a float, or a None in the list, or a high value.
pass
if None in qualities:
raise TypeError("A quality value of None was found")
if max(qualities) >= 93.5:
warnings.warn(
"Data loss - max PHRED quality 93 in Sanger FASTQ", BiopythonWarning
)
# This will apply the truncation at 93, giving max ASCII 126
return "".join(
chr(min(126, int(round(qp)) + SANGER_SCORE_OFFSET)) for qp in qualities
)
# Fall back on the Solexa scores...
try:
qualities = record.letter_annotations["solexa_quality"]
except KeyError:
raise ValueError(
"No suitable quality scores found in "
"letter_annotations of SeqRecord (id=%s)." % record.id
) from None
# Try and use the precomputed mapping:
try:
return "".join(_solexa_to_sanger_quality_str[qs] for qs in qualities)
except KeyError:
# Either no PHRED scores, or something odd like a float or None
pass
if None in qualities:
raise TypeError("A quality value of None was found")
# Must do this the slow way, first converting the PHRED scores into
# Solexa scores:
if max(qualities) >= 93.5:
warnings.warn(
"Data loss - max PHRED quality 93 in Sanger FASTQ", BiopythonWarning
)
# This will apply the truncation at 93, giving max ASCII 126
return "".join(
chr(min(126, int(round(phred_quality_from_solexa(qs))) + SANGER_SCORE_OFFSET))
for qs in qualities
)
# Only map 0 to 62, we need to give a warning on truncating at 62
assert 62 + SOLEXA_SCORE_OFFSET == 126
_phred_to_illumina_quality_str = {
qp: chr(qp + SOLEXA_SCORE_OFFSET) for qp in range(62 + 1)
}
# Only map -5 to 62, we need to give a warning on truncating at 62
_solexa_to_illumina_quality_str = {
qs: chr(int(round(phred_quality_from_solexa(qs))) + SOLEXA_SCORE_OFFSET)
for qs in range(-5, 62 + 1)
}
def _get_illumina_quality_str(record: SeqRecord) -> str:
"""Return an Illumina 1.3 to 1.7 FASTQ encoded quality string (PRIVATE).
Notice that due to the limited range of printable ASCII characters, a
PHRED quality of 62 is the maximum that can be held in an Illumina FASTQ
file (using ASCII 126, the tilde). This function will issue a warning
in this situation.
"""
# TODO - This functions works and is fast, but it is also ugly
# and there is considerable repetition of code for the other
# two FASTQ variants.
try:
# These take priority (in case both Solexa and PHRED scores found)
qualities = record.letter_annotations["phred_quality"]
except KeyError:
# Fall back on solexa scores...
pass
else:
# Try and use the precomputed mapping:
try:
return "".join(_phred_to_illumina_quality_str[qp] for qp in qualities)
except KeyError:
# Could be a float, or a None in the list, or a high value.
pass
if None in qualities:
raise TypeError("A quality value of None was found")
if max(qualities) >= 62.5:
warnings.warn(
"Data loss - max PHRED quality 62 in Illumina FASTQ", BiopythonWarning
)
# This will apply the truncation at 62, giving max ASCII 126
return "".join(
chr(min(126, int(round(qp)) + SOLEXA_SCORE_OFFSET)) for qp in qualities
)
# Fall back on the Solexa scores...
try:
qualities = record.letter_annotations["solexa_quality"]
except KeyError:
raise ValueError(
"No suitable quality scores found in "
"letter_annotations of SeqRecord (id=%s)." % record.id
) from None
# Try and use the precomputed mapping:
try:
return "".join(_solexa_to_illumina_quality_str[qs] for qs in qualities)
except KeyError:
# Either no PHRED scores, or something odd like a float or None
pass
if None in qualities:
raise TypeError("A quality value of None was found")
# Must do this the slow way, first converting the PHRED scores into
# Solexa scores:
if max(qualities) >= 62.5:
warnings.warn(
"Data loss - max PHRED quality 62 in Illumina FASTQ", BiopythonWarning
)
# This will apply the truncation at 62, giving max ASCII 126
return "".join(
chr(min(126, int(round(phred_quality_from_solexa(qs))) + SOLEXA_SCORE_OFFSET))
for qs in qualities
)
# Only map 0 to 62, we need to give a warning on truncating at 62
assert 62 + SOLEXA_SCORE_OFFSET == 126
_solexa_to_solexa_quality_str = {
qs: chr(min(126, qs + SOLEXA_SCORE_OFFSET)) for qs in range(-5, 62 + 1)
}
# Only map -5 to 62, we need to give a warning on truncating at 62
_phred_to_solexa_quality_str = {
qp: chr(min(126, int(round(solexa_quality_from_phred(qp))) + SOLEXA_SCORE_OFFSET))
for qp in range(62 + 1)
}
def _get_solexa_quality_str(record: SeqRecord) -> str:
"""Return a Solexa FASTQ encoded quality string (PRIVATE).
Notice that due to the limited range of printable ASCII characters, a
Solexa quality of 62 is the maximum that can be held in a Solexa FASTQ
file (using ASCII 126, the tilde). This function will issue a warning
in this situation.
"""
# TODO - This functions works and is fast, but it is also ugly
# and there is considerable repetition of code for the other
# two FASTQ variants.
try:
# These take priority (in case both Solexa and PHRED scores found)
qualities = record.letter_annotations["solexa_quality"]
except KeyError:
# Fall back on PHRED scores...
pass
else:
# Try and use the precomputed mapping:
try:
return "".join(_solexa_to_solexa_quality_str[qs] for qs in qualities)
except KeyError:
# Could be a float, or a None in the list, or a high value.
pass
if None in qualities:
raise TypeError("A quality value of None was found")
if max(qualities) >= 62.5:
warnings.warn(
"Data loss - max Solexa quality 62 in Solexa FASTQ", BiopythonWarning
)
# This will apply the truncation at 62, giving max ASCII 126
return "".join(
chr(min(126, int(round(qs)) + SOLEXA_SCORE_OFFSET)) for qs in qualities
)
# Fall back on the PHRED scores...
try:
qualities = record.letter_annotations["phred_quality"]
except KeyError:
raise ValueError(
"No suitable quality scores found in "
"letter_annotations of SeqRecord (id=%s)." % record.id
) from None
# Try and use the precomputed mapping:
try:
return "".join(_phred_to_solexa_quality_str[qp] for qp in qualities)
except KeyError:
# Either no PHRED scores, or something odd like a float or None
# or too big to be in the cache
pass
if None in qualities:
raise TypeError("A quality value of None was found")
# Must do this the slow way, first converting the PHRED scores into
# Solexa scores:
if max(qualities) >= 62.5:
warnings.warn(
"Data loss - max Solexa quality 62 in Solexa FASTQ", BiopythonWarning
)
return "".join(
chr(min(126, int(round(solexa_quality_from_phred(qp))) + SOLEXA_SCORE_OFFSET))
for qp in qualities
)
# TODO - Default to nucleotide or even DNA?
def FastqGeneralIterator(source: _TextIOSource) -> Iterator[tuple[str, str, str]]:
"""Iterate over Fastq records as string tuples (not as SeqRecord objects).
Arguments:
- source - input stream opened in text mode, or a path to a file
This code does not try to interpret the quality string numerically. It
just returns tuples of the title, sequence and quality as strings. For
the sequence and quality, any whitespace (such as new lines) is removed.
Our SeqRecord based FASTQ iterators call this function internally, and then
turn the strings into a SeqRecord objects, mapping the quality string into
a list of numerical scores. If you want to do a custom quality mapping,
then you might consider calling this function directly.
For parsing FASTQ files, the title string from the "@" line at the start
of each record can optionally be omitted on the "+" lines. If it is
repeated, it must be identical.
The sequence string and the quality string can optionally be split over
multiple lines, although several sources discourage this. In comparison,
for the FASTA file format line breaks between 60 and 80 characters are
the norm.
**WARNING** - Because the "@" character can appear in the quality string,
this can cause problems as this is also the marker for the start of
a new sequence. In fact, the "+" sign can also appear as well. Some
sources recommended having no line breaks in the quality to avoid this,
but even that is not enough, consider this example::
@071113_EAS56_0053:1:1:998:236
TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA
+071113_EAS56_0053:1:1:998:236
IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III
@071113_EAS56_0053:1:1:182:712
ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG
+
@IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+
@071113_EAS56_0053:1:1:153:10
TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT
+
IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6
@071113_EAS56_0053:1:3:990:501
TGGGAGGTTTTATGTGGA
AAGCAGCAATGTACAAGA
+
IIIIIII.IIIIII1@44
@-7.%<&+/$/%4(++(%
This is four PHRED encoded FASTQ entries originally from an NCBI source
(given the read length of 36, these are probably Solexa Illumina reads where
the quality has been mapped onto the PHRED values).
This example has been edited to illustrate some of the nasty things allowed
in the FASTQ format. Firstly, on the "+" lines most but not all of the
(redundant) identifiers are omitted. In real files it is likely that all or
none of these extra identifiers will be present.
Secondly, while the first three sequences have been shown without line
breaks, the last has been split over multiple lines. In real files any line
breaks are likely to be consistent.
Thirdly, some of the quality string lines start with an "@" character. For
the second record this is unavoidable. However for the fourth sequence this
only happens because its quality string is split over two lines. A naive
parser could wrongly treat any line starting with an "@" as the beginning of
a new sequence! This code copes with this possible ambiguity by keeping
track of the length of the sequence which gives the expected length of the
quality string.
Using this tricky example file as input, this short bit of code demonstrates
what this parsing function would return:
>>> with open("Quality/tricky.fastq") as handle:
... for (title, sequence, quality) in FastqGeneralIterator(handle):
... print(title)
... print("%s %s" % (sequence, quality))
...
071113_EAS56_0053:1:1:998:236
TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III
071113_EAS56_0053:1:1:182:712
ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG @IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+
071113_EAS56_0053:1:1:153:10
TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6
071113_EAS56_0053:1:3:990:501
TGGGAGGTTTTATGTGGAAAGCAGCAATGTACAAGA IIIIIII.IIIIII1@44@-7.%<&+/$/%4(++(%
Finally we note that some sources state that the quality string should
start with "!" (which using the PHRED mapping means the first letter always
has a quality score of zero). This rather restrictive rule is not widely
observed, so is therefore ignored here. One plus point about this "!" rule
is that (provided there are no line breaks in the quality sequence) it
would prevent the above problem with the "@" character.
"""
with as_handle(source) as handle:
if handle.read(0) != "":
raise StreamModeError("Fastq files must be opened in text mode") from None
try:
line = next(handle)
except StopIteration:
return # Premature end of file, or just empty?
while True:
if line[0] != "@":
raise ValueError(
"Records in Fastq files should start with '@' character"
)
title_line = line[1:].rstrip()
seq_string = ""
# There will now be one or more sequence lines; keep going until we
# find the "+" marking the quality line:
for line in handle:
if line[0] == "+":
break
seq_string += line.rstrip()
else:
if seq_string:
raise ValueError("End of file without quality information.")
else:
raise ValueError("Unexpected end of file")
# The title here is optional, but if present must match!
second_title = line[1:].rstrip()
if second_title and second_title != title_line:
raise ValueError("Sequence and quality captions differ.")
# This is going to slow things down a little, but assuming
# this isn't allowed we should try and catch it here:
if " " in seq_string or "\t" in seq_string:
raise ValueError("Whitespace is not allowed in the sequence.")
seq_len = len(seq_string)
# There will now be at least one line of quality data, followed by
# another sequence, or EOF
line = None
quality_string = ""
for line in handle:
if line[0] == "@":
# This COULD be the start of a new sequence. However, it MAY just
# be a line of quality data which starts with a "@" character. We
# should be able to check this by looking at the sequence length
# and the amount of quality data found so far.
if len(quality_string) >= seq_len:
# We expect it to be equal if this is the start of a new record.
# If the quality data is longer, we'll raise an error below.
break
# Continue - its just some (more) quality data.
quality_string += line.rstrip()
else:
if line is None:
raise ValueError("Unexpected end of file")
line = None
if seq_len != len(quality_string):
raise ValueError(
"Lengths of sequence and quality values differs for %s (%i and %i)."
% (title_line, seq_len, len(quality_string))
)
# Return the record and then continue...
yield (title_line, seq_string, quality_string)
if line is None:
break
class FastqIteratorAbstractBaseClass(SequenceIterator[str]):
"""Abstract base class for FASTQ file parsers."""
modes = "t"