--- /dev/null
+# Copyright 2009 by Peter Cock. All rights reserved.
+# This code is part of the Biopython distribution and governed by its
+# license. Please see the LICENSE file that should have been included
+# as part of this package.
+#
+# This module is for reading and writing FASTQ and QUAL format files as
+# SeqRecord objects, and is expected to be used via the Bio.SeqIO API.
+
+"""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 quality value between 0 and 90 to each called base using a logged
+transformation of the error probability, Q = -10 log10( Pe ), for example::
+
+ Pe = 0.0, Q = 0
+ Pe = 0.1, Q = 10
+ Pe = 0.01, Q = 20
+ ...
+ Pe = 0.00000001, Q = 80
+ Pe = 0.000000001, Q = 90
+
+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". 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. Reasonable
+documentation exists at: http://maq.sourceforge.net/fastq.shtml
+
+Solexa/Illumina quality scores use Q = - 10 log10 ( Pe / (1-Pe) ), which can
+be negative or easily exceed 90. PHRED scores and Solexa scores are NOT
+interchangeable (but a reasonable mapping can be achieved between them).
+Confusingly Solexa produces a FASTQ like file but using their own score
+mapping instead.
+
+Also note that Roche 454 sequencers can output files in the QUAL format, and
+thankfully they use PHREP style scores like Sanger. 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.
+
+You are expected to use this module via the Bio.SeqIO functions, with the
+following format names:
+ - "fastq" means Sanger style FASTQ files using PHRED scores.
+ - "fastq-solexa" means Solexa/Illumina style FASTQ files.
+ - "qual" means simple quality files using PHRED scores.
+
+For example, consider the following short FASTQ file (extracted from a real
+NCBI dataset)::
+
+ @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 NCBI have
+followed the Sanger FASTQ convention and this actually uses PHRED style
+qualities. 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(open("Quality/example.fastq"), "fastq") :
+ ... print 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', SingleLetterAlphabet())
+ >>> 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 you can 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,
+
+ >>> print record.format("fastq")
+ @EAS54_6_R1_2_1_443_348
+ GTTGCTTCTGGCGTGGGTGGGGGGG
+ +
+ ;;;;;;;;;;;9;7;;.7;393333
+ <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>
+
+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', SingleLetterAlphabet())
+ >>> 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(open("Quality/example.fastq"), "fastq")
+ >>> out_handle = open("Quality/temp.qual", "w")
+ >>> SeqIO.write(record_iterator, out_handle, "qual")
+ 3
+ >>> out_handle.close()
+
+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(open("Quality/temp.qual"), "qual") :
+ ... print record.id, record.seq
+ EAS54_6_R1_2_1_413_324 ?????????????????????????
+ EAS54_6_R1_2_1_540_792 ?????????????????????????
+ EAS54_6_R1_2_1_443_348 ?????????????????????????
+
+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
+ UnknownSeq(25, alphabet = SingleLetterAlphabet(), character = '?')
+ >>> 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(open("Quality/example.fasta"), "fasta"))
+ >>> for rec in SeqIO.parse(open("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" for the Sanger standard using PHRED values, or "fastq-solexa"
+for the Solexa/Illumina variant), as this cannot be detected reliably
+automatically.
+"""
+__docformat__ = "epytext en" #Don't just use plain text in epydoc API pages!
+
+#See also http://blog.malde.org/index.php/2008/09/09/the-fastq-file-format-for-sequences/
+
+from Bio.Alphabet import single_letter_alphabet
+from Bio.Seq import Seq, UnknownSeq
+from Bio.SeqRecord import SeqRecord
+from Interfaces import SequentialSequenceWriter
+from math import log
+
+# 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) :
+ """Covert a PHRED quality (range 0 to about 90) to a Solexa quality.
+
+ 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(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(1),2)
+ -5.87
+ """
+ return 10*log(10**(phred_quality/10.0) - 1, 10)
+
+def phred_quality_from_solexa(solexa_quality) :
+ """Convert a Solexa quality (which can be negative) to a PHRED quality.
+
+ 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(-10),2)
+ 0.41
+ """
+ return 10*log(10**(solexa_quality/10.0) + 1, 10)
+
+def _get_phred_quality(record) :
+ """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)
+
+def _get_solexa_quality(record) :
+ """Extract Solexa qualities from a SeqRecord's letter_annotations (PRIVATE).
+
+ If there are no Solexa qualities, but there are PHRED qualities, those are
+ used instead after conversion.
+ """
+ try :
+ return record.letter_annotations["solexa_quality"]
+ except KeyError :
+ pass
+ try :
+ return [solexa_quality_from_phred(q) for \
+ q in record.letter_annotations["phred_quality"]]
+ except KeyError :
+ raise ValueError("No suitable quality scores found in letter_annotation "
+ "of SeqRecord (id=%s)." % record.id)
+
+
+#TODO - Default to nucleotide or even DNA?
+def FastqGeneralIterator(handle) :
+ """Iterate over Fastq records as string tuples (not as SeqRecord objects).
+
+ 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 Illumna 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 ommited. 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:
+
+ >>> handle = open("Quality/tricky.fastq", "rU")
+ >>> for (title, sequence, quality) in FastqGeneralIterator(handle) :
+ ... print title
+ ... print 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(++(%
+ >>> handle.close()
+
+ 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.
+ """
+ #Skip any text before the first record (e.g. blank lines, comments?)
+ while True :
+ line = handle.readline()
+ if line == "" : return #Premature end of file, or just empty?
+ if line[0] == "@" :
+ break
+
+ while True :
+ if line[0]!="@" :
+ raise ValueError("Records in Fastq files should start with '@' character")
+ title_line = line[1:].rstrip()
+
+ seq_lines = []
+ line = handle.readline()
+ while True:
+ if not line :
+ raise ValueError("End of file without quality information.")
+ if line[0] == "+":
+ #The title here is optional, but if present must match!
+ if line[1:].rstrip() and line[1:].rstrip() != title_line :
+ raise ValueError("Sequence and quality captions differ.")
+ break
+ seq_lines.extend(line.split()) #removes any whitespace
+ line = handle.readline()
+
+ seq_string = "".join(seq_lines)
+ del seq_lines
+
+ quality_lines = []
+ line = handle.readline()
+ while True:
+ if not line : break
+ 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("".join(quality_lines)) >= len(seq_string) :
+ #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) sequence data.
+
+ quality_lines.extend(line.split()) #removes any whitespace
+ line = handle.readline()
+
+ quality_string = "".join(quality_lines)
+ del quality_lines
+
+ if len(seq_string) != len(quality_string) :
+ raise ValueError("Lengths of sequence and quality values differs "
+ " for %s (%i and %i)." \
+ % (title_line, len(seq_string), len(quality_string)))
+
+ #Return the record and then continue...
+ yield (title_line, seq_string, quality_string)
+ if not line : return #StopIteration at end of file
+ assert False, "Should not reach this line"
+
+#This is a generator function!
+def FastqPhredIterator(handle, alphabet = single_letter_alphabet, title2ids = None) :
+ """Generator function to iterate over FASTQ records (as SeqRecord objects).
+
+ - handle - input file
+ - alphabet - optional alphabet
+ - title2ids - A function that, when given the title line from the FASTQ
+ file (without the beginning >), will return the id, name and
+ description (in that order) for the record as a tuple of
+ strings. If this is not given, then the entire title line
+ will be used as the description, and the first word as the
+ id and name.
+
+ Note that use of title2ids matches that of Bio.SeqIO.FastaIO.
+
+ For each sequence in a (Sanger style) FASTQ file there is a matching string
+ encoding the PHRED qualities (integers between 0 and about 90) using ASCII
+ values with an offset of 33.
+
+ For example, consider a file containing three short reads::
+
+ @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
+
+ For each sequence (e.g. "CCCTTCTTGTCTTCAGCGTTTCTCC") there is a matching
+ string encoding the PHRED qualities using a ASCI values with an offset of
+ 33 (e.g. ";;3;;;;;;;;;;;;7;;;;;;;88").
+
+ Using this module directly you might run:
+
+ >>> handle = open("Quality/example.fastq", "rU")
+ >>> for record in FastqPhredIterator(handle) :
+ ... print 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
+ >>> handle.close()
+
+ Typically however, you would call this via Bio.SeqIO instead with "fastq" as
+ the format:
+
+ >>> from Bio import SeqIO
+ >>> handle = open("Quality/example.fastq", "rU")
+ >>> for record in SeqIO.parse(handle, "fastq") :
+ ... print 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
+ >>> handle.close()
+
+ If you want to look at the qualities, they are record in each record's
+ per-letter-annotation dictionary as a simple list of integers:
+
+ >>> 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]
+ """
+ for title_line, seq_string, quality_string in FastqGeneralIterator(handle) :
+ if title2ids :
+ id, name, descr = title2ids(title_line)
+ else :
+ descr = title_line
+ id = descr.split()[0]
+ name = id
+ record = SeqRecord(Seq(seq_string, alphabet),
+ id=id, name=name, description=descr)
+
+ assert SANGER_SCORE_OFFSET == ord("!")
+ #According to BioPerl documentation at least, the first character should
+ #be an "!" (and therefore quality zero). This seems crazy - what if the
+ #sequence has been trimmed to remove any poor quality sequence? In any
+ #case real examples from the NCBI don't follow this practice, so we
+ #won't enforce it here.
+ #e.g. ftp://ftp.ncbi.nih.gov/pub/TraceDB/ShortRead/SRA000271/fastq/200x36x36-071113_EAS56_0053-s_1_1.fastq.gz
+ #
+ #if quality_string[0] != "!" :
+ # raise ValueError("The quality string should always start with a ! character.")
+ qualities = [ord(letter)-SANGER_SCORE_OFFSET for letter in quality_string]
+ if qualities :
+ if min(qualities) < 0 or max(qualities) > 90 :
+ raise ValueError("Quality score outside 0 to 90 found - these are perhaps "
+ "in a Solexa/Illumina format, not the Sanger FASTQ format "
+ "which uses PHRED scores.")
+ record.letter_annotations["phred_quality"] = qualities
+ yield record
+
+#This is a generator function!
+def FastqSolexaIterator(handle, alphabet = single_letter_alphabet, title2ids = None) :
+ """Parsing the Solexa/Illumina FASTQ like files (which differ in the quality mapping).
+
+ The optional arguments are the same as those for the FastqPhredIterator.
+
+ For each sequence in Solexa/Illumina FASTQ files there is a matching string
+ encoding the Solexa integer qualities using ASCII values with an offset
+ of 64. Solexa scores are scaled differently to PHRED scores, and Biopython
+ will NOT perform any automatic conversion when loading.
+
+ For example, consider a file containing these five records::
+
+ @SLXA-B3_649_FC8437_R1_1_1_610_79
+ GATGTGCAATACCTTTGTAGAGGAA
+ +SLXA-B3_649_FC8437_R1_1_1_610_79
+ YYYYYYYYYYYYYYYYYYWYWYYSU
+ @SLXA-B3_649_FC8437_R1_1_1_397_389
+ GGTTTGAGAAAGAGAAATGAGATAA
+ +SLXA-B3_649_FC8437_R1_1_1_397_389
+ YYYYYYYYYWYYYYWWYYYWYWYWW
+ @SLXA-B3_649_FC8437_R1_1_1_850_123
+ GAGGGTGTTGATCATGATGATGGCG
+ +SLXA-B3_649_FC8437_R1_1_1_850_123
+ YYYYYYYYYYYYYWYYWYYSYYYSY
+ @SLXA-B3_649_FC8437_R1_1_1_362_549
+ GGAAACAAAGTTTTTCTCAACATAG
+ +SLXA-B3_649_FC8437_R1_1_1_362_549
+ YYYYYYYYYYYYYYYYYYWWWWYWY
+ @SLXA-B3_649_FC8437_R1_1_1_183_714
+ GTATTATTTAATGGCATACACTCAA
+ +SLXA-B3_649_FC8437_R1_1_1_183_714
+ YYYYYYYYYYWYYYYWYWWUWWWQQ
+
+ Using this module directly you might run:
+
+ >>> handle = open("Quality/solexa_example.fastq", "rU")
+ >>> for record in FastqSolexaIterator(handle) :
+ ... print record.id, record.seq
+ SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA
+ SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA
+ SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG
+ SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG
+ SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA
+ >>> handle.close()
+
+ Typically however, you would call this via Bio.SeqIO instead with "fastq" as
+ the format:
+
+ >>> from Bio import SeqIO
+ >>> handle = open("Quality/solexa_example.fastq", "rU")
+ >>> for record in SeqIO.parse(handle, "fastq-solexa") :
+ ... print record.id, record.seq
+ SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA
+ SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA
+ SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG
+ SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG
+ SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA
+ >>> handle.close()
+
+ If you want to look at the qualities, they are recorded in each record's
+ per-letter-annotation dictionary as a simple list of integers:
+
+ >>> print record.letter_annotations["solexa_quality"]
+ [25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 23, 25, 25, 25, 25, 23, 25, 23, 23, 21, 23, 23, 23, 17, 17]
+
+ These scores aren't very good, but they are high enough that they map
+ almost exactly onto PHRED scores:
+
+ >>> print "%0.2f" % phred_quality_from_solexa(25)
+ 25.01
+
+ Let's look at another example read which is even worse, where there are
+ more noticeable differences between the Solexa and PHRED scores::
+
+ @slxa_0013_1_0001_24
+ ACAAAAATCACAAGCATTCTTATACACC
+ +slxa_0013_1_0001_24
+ ??????????????????:??<?<-6%.
+
+ Again, you would typically use Bio.SeqIO to read this file in (rather than
+ calling the Bio.SeqIO.QualtityIO module directly). Most FASTQ files will
+ contain thousands of reads, so you would normally use Bio.SeqIO.parse()
+ as shown above. This example has only as one entry, so instead we can
+ use the Bio.SeqIO.read() function:
+
+ >>> from Bio import SeqIO
+ >>> handle = open("Quality/solexa.fastq", "rU")
+ >>> record = SeqIO.read(handle, "fastq-solexa")
+ >>> handle.close()
+ >>> print record.id, record.seq
+ slxa_0013_1_0001_24 ACAAAAATCACAAGCATTCTTATACACC
+ >>> print record.letter_annotations["solexa_quality"]
+ [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -6, -1, -1, -4, -1, -4, -19, -10, -27, -18]
+
+ These quality scores are so low that when converted from the Solexa scheme
+ into PHRED scores they look quite different:
+
+ >>> print "%0.2f" % phred_quality_from_solexa(-1)
+ 2.54
+
+ Note you can use the Bio.SeqIO.write() function or the SeqRecord's format
+ method to output the record(s):
+
+ >>> print record.format("fastq-solexa")
+ @slxa_0013_1_0001_24
+ ACAAAAATCACAAGCATTCTTATACACC
+ +
+ ??????????????????:??<?<-6%.
+ <BLANKLINE>
+
+ Note this output is slightly different from the input file as Biopython
+ has left out the optional repetition of the sequence identifier on the "+"
+ line. If you want the to use PHRED scores, use "fastq" or "qual" as the
+ output format instead, and Biopython will do the conversion for you:
+
+ >>> print record.format("fastq")
+ @slxa_0013_1_0001_24
+ ACAAAAATCACAAGCATTCTTATACACC
+ +
+ $$$$$$$$$$$$$$$$$$"$$"$"!!!!
+ <BLANKLINE>
+
+ >>> print record.format("qual")
+ >slxa_0013_1_0001_24
+ 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 1 3 1 0 0 0 0
+ <BLANKLINE>
+ """
+ for title_line, seq_string, quality_string in FastqGeneralIterator(handle) :
+ if title2ids :
+ id, name, descr = title_line
+ else :
+ descr = title_line
+ id = descr.split()[0]
+ name = id
+ record = SeqRecord(Seq(seq_string, alphabet),
+ id=id, name=name, description=descr)
+ qualities = [ord(letter)-SOLEXA_SCORE_OFFSET for letter in quality_string]
+ #DO NOT convert these into PHRED qualities automatically!
+ record.letter_annotations["solexa_quality"] = qualities
+ yield record
+
+def QualPhredIterator(handle, alphabet = single_letter_alphabet, title2ids = None) :
+ """For QUAL files which include PHRED quality scores, but no sequence.
+
+ For example, consider this short QUAL file::
+
+ >EAS54_6_R1_2_1_413_324
+ 26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26
+ 26 26 26 23 23
+ >EAS54_6_R1_2_1_540_792
+ 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26
+ 26 18 26 23 18
+ >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
+
+ Using this module directly you might run:
+
+ >>> handle = open("Quality/example.qual", "rU")
+ >>> for record in QualPhredIterator(handle) :
+ ... print record.id, record.seq
+ EAS54_6_R1_2_1_413_324 ?????????????????????????
+ EAS54_6_R1_2_1_540_792 ?????????????????????????
+ EAS54_6_R1_2_1_443_348 ?????????????????????????
+ >>> handle.close()
+
+ Typically however, you would call this via Bio.SeqIO instead with "qual"
+ as the format:
+
+ >>> from Bio import SeqIO
+ >>> handle = open("Quality/example.qual", "rU")
+ >>> for record in SeqIO.parse(handle, "qual") :
+ ... print record.id, record.seq
+ EAS54_6_R1_2_1_413_324 ?????????????????????????
+ EAS54_6_R1_2_1_540_792 ?????????????????????????
+ EAS54_6_R1_2_1_443_348 ?????????????????????????
+ >>> handle.close()
+
+ Becase QUAL files don't contain the sequence string itself, the seq
+ property is set to an UnknownSeq object. As no alphabet was given, this
+ has defaulted to a generic single letter alphabet and the character "?"
+ used.
+
+ By specifying a nucleotide alphabet, "N" is used instead:
+
+ >>> from Bio import SeqIO
+ >>> from Bio.Alphabet import generic_dna
+ >>> handle = open("Quality/example.qual", "rU")
+ >>> for record in SeqIO.parse(handle, "qual", alphabet=generic_dna) :
+ ... print record.id, record.seq
+ EAS54_6_R1_2_1_413_324 NNNNNNNNNNNNNNNNNNNNNNNNN
+ EAS54_6_R1_2_1_540_792 NNNNNNNNNNNNNNNNNNNNNNNNN
+ EAS54_6_R1_2_1_443_348 NNNNNNNNNNNNNNNNNNNNNNNNN
+ >>> handle.close()
+
+ However, the quality scores themselves are available as a list of integers
+ in each record's per-letter-annotation:
+
+ >>> 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 still slice one of these SeqRecord objects with an UnknownSeq:
+
+ >>> sub_record = record[5:10]
+ >>> print sub_record.id, sub_record.letter_annotations["phred_quality"]
+ EAS54_6_R1_2_1_443_348 [26, 26, 26, 26, 26]
+ """
+ #Skip any text before the first record (e.g. blank lines, comments)
+ while True :
+ line = handle.readline()
+ if line == "" : return #Premature end of file, or just empty?
+ if line[0] == ">" :
+ break
+
+ while True :
+ if line[0]!=">" :
+ raise ValueError("Records in Fasta files should start with '>' character")
+ if title2ids :
+ id, name, descr = title2ids(line[1:].rstrip())
+ else :
+ descr = line[1:].rstrip()
+ id = descr.split()[0]
+ name = id
+
+ qualities = []
+ line = handle.readline()
+ while True:
+ if not line : break
+ if line[0] == ">": break
+ qualities.extend([int(word) for word in line.split()])
+ line = handle.readline()
+
+ if qualities :
+ if min(qualities) < 0 or max(qualities) > 90 :
+ raise ValueError(("Quality score range for %s is %i to %i, outside the " \
+ +"expected 0 to 90. Perhaps these are Solexa/Illumina " \
+ +"scores, and not PHRED scores?") \
+ % (id, min(qualities), max(qualities)))
+
+ #Return the record and then continue...
+ record = SeqRecord(UnknownSeq(len(qualities), alphabet),
+ id = id, name = name, description = descr)
+ record.letter_annotations["phred_quality"] = qualities
+ yield record
+
+ if not line : return #StopIteration
+ assert False, "Should not reach this line"
+
+class FastqPhredWriter(SequentialSequenceWriter):
+ """Class to write FASTQ format files (using PHRED quality scores).
+
+ Although you can use this class directly, you are strongly encouraged
+ to use the Bio.SeqIO.write() function instead. For example, this code
+ reads in a FASTQ (PHRED) file and re-saves it as another FASTQ (PHRED)
+ file:
+
+ >>> from Bio import SeqIO
+ >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq")
+ >>> out_handle = open("Quality/temp.fastq", "w")
+ >>> SeqIO.write(record_iterator, out_handle, "fastq")
+ 3
+ >>> out_handle.close()
+
+ You might want to do this if the original file included extra line breaks,
+ which while valid may not be supported by all tools. The output file from
+ Biopython will have each sequence on a single line, and each quality
+ string on a single line (which is considered desirable for maximum
+ compatibility).
+
+ In this next example, a Solexa FASTQ file is converted into a standard
+ Sanger style FASTQ file using PHRED qualities:
+
+ >>> from Bio import SeqIO
+ >>> record_iterator = SeqIO.parse(open("Quality/solexa.fastq"), "fastq-solexa")
+ >>> out_handle = open("Quality/temp.fastq", "w")
+ >>> SeqIO.write(record_iterator, out_handle, "fastq")
+ 1
+ >>> out_handle.close()
+
+ This code is also called if you use the .format("fastq") method of a
+ SeqRecord.
+
+ P.S. To avoid cluttering up your working directory, you can delete this
+ temporary file now:
+
+ >>> import os
+ >>> os.remove("Quality/temp.fastq")
+
+ """
+ def write_record(self, record):
+ """Write a single FASTQ record to the file."""
+ assert self._header_written
+ assert not self._footer_written
+ self._record_written = True
+
+ #TODO - Is an empty sequence allowed in FASTQ format?
+ assert SANGER_SCORE_OFFSET == ord("!")
+ #This rounds to the nearest integer:
+ qualities = "".join([chr(int(round(q+SANGER_SCORE_OFFSET,0))) for q \
+ in _get_phred_quality(record)])
+ if record.seq is None:
+ raise ValueError("No sequence for record %s" % record.id)
+ if len(qualities) != len(record) :
+ raise ValueError("Record %s has sequence length %i but %i quality scores" \
+ % (record.id, len(record), len(qualities)))
+
+ title = self.clean(record.id) #TODO - add the description too? cf Fasta output
+ self.handle.write("@%s\n%s\n+\n%s\n" % (title, record.seq, qualities))
+
+class QualPhredWriter(SequentialSequenceWriter):
+ """Class to write QUAL format files (using PHRED quality scores).
+
+ Although you can use this class directly, you are strongly encouraged
+ to use the Bio.SeqIO.write() function instead. For example, this code
+ reads in a FASTQ file and saves the quality scores into a QUAL file:
+
+ >>> from Bio import SeqIO
+ >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq")
+ >>> out_handle = open("Quality/temp.qual", "w")
+ >>> SeqIO.write(record_iterator, out_handle, "qual")
+ 3
+ >>> out_handle.close()
+
+ This code is also called if you use the .format("qual") method of a
+ SeqRecord.
+
+ P.S. Don't forget to clean up the temp file if you don't need it anymore:
+
+ >>> import os
+ >>> os.remove("Quality/temp.qual")
+ """
+ def __init__(self, handle, wrap=60, record2title=None):
+ """Create a QUAL writer.
+
+ Arguments:
+ - handle - Handle to an output file, e.g. as returned
+ by open(filename, "w")
+ - wrap - Optional line length used to wrap sequence lines.
+ Defaults to wrapping the sequence at 60 characters
+ Use zero (or None) for no wrapping, giving a single
+ long line for the sequence.
+ - record2title - Optional function to return the text to be
+ used for the title line of each record. By default
+ a combination of the record.id and record.description
+ is used. If the record.description starts with the
+ record.id, then just the record.description is used.
+
+ The record2title argument is present for consistency with the
+ Bio.SeqIO.FastaIO writer class.
+ """
+ SequentialSequenceWriter.__init__(self, handle)
+ #self.handle = handle
+ self.wrap = None
+ if wrap :
+ if wrap < 1 :
+ raise ValueError
+ self.wrap = wrap
+ self.record2title = record2title
+
+ def write_record(self, record):
+ """Write a single QUAL record to the file."""
+ assert self._header_written
+ assert not self._footer_written
+ self._record_written = True
+
+ if self.record2title :
+ title=self.clean(self.record2title(record))
+ else :
+ id = self.clean(record.id)
+ description = self.clean(record.description)
+
+ #if description[:len(id)]==id :
+ if description and description.split(None,1)[0]==id :
+ #The description includes the id at the start
+ title = description
+ else :
+ title = "%s %s" % (id, description)
+
+ assert "\n" not in title
+ assert "\r" not in title
+ self.handle.write(">%s\n" % title)
+
+ #This rounds to the nearest integer.
+ #TODO - can we put a float in a qual file?
+ qualities = [("%i" % round(q,0)) for q in _get_phred_quality(record)]
+
+ if self.wrap :
+ while qualities :
+ line=qualities.pop(0)
+ while qualities \
+ and len(line) + 1 + len(qualities[0]) < self.wrap :
+ line += " " + qualities.pop(0)
+ self.handle.write(line + "\n")
+ else :
+ data = " ".join(qualities)
+ self.handle.write(data + "\n")
+
+class FastqSolexaWriter(SequentialSequenceWriter):
+ """Class to write FASTQ format files (using Solexa quality scores).
+
+ Although you can use this class directly, you are strongly encouraged
+ to use the Bio.SeqIO.write() function instead. For example, this code
+ reads in a FASTQ file and re-saves it as another FASTQ file:
+
+ >>> from Bio import SeqIO
+ >>> record_iterator = SeqIO.parse(open("Quality/solexa.fastq"), "fastq-solexa")
+ >>> out_handle = open("Quality/temp.fastq", "w")
+ >>> SeqIO.write(record_iterator, out_handle, "fastq-solexa")
+ 1
+ >>> out_handle.close()
+
+ You might want to do this if the original file included extra line
+ breaks, which (while valid) may not be supported by all tools. The
+ output file from Biopython will have each sequence on a single line, and
+ each quality string on a single line (which is considered desirable for
+ maximum compatibility).
+
+ This code is also called if you use the .format("fastq-solexa") method of
+ a SeqRecord.
+
+ P.S. Don't forget to delete the temp file if you don't need it anymore:
+
+ >>> import os
+ >>> os.remove("Quality/temp.fastq")
+ """
+ def write_record(self, record):
+ """Write a single FASTQ record to the file."""
+ assert self._header_written
+ assert not self._footer_written
+ self._record_written = True
+
+ #TODO - Is an empty sequence allowed in FASTQ format?
+ qualities = "".join([chr(int(round(q+SOLEXA_SCORE_OFFSET,0))) for q \
+ in _get_solexa_quality(record)])
+ if record.seq is None:
+ raise ValueError("No sequence for record %s" % record.id)
+ if len(qualities) != len(record) :
+ raise ValueError("Record %s has sequence length %i but %i quality scores" \
+ % (record.id, len(record), len(qualities)))
+
+ title = self.clean(record.id) #TODO - add the description too? cf Fasta output
+ self.handle.write("@%s\n%s\n+\n%s\n" % (title, record.seq, qualities))
+
+def PairedFastaQualIterator(fasta_handle, qual_handle, alphabet = single_letter_alphabet, title2ids = None) :
+ """Iterate over matched FASTA and QUAL files as SeqRecord objects.
+
+ For example, consider this short QUAL file::
+
+ >EAS54_6_R1_2_1_413_324
+ 26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26
+ 26 26 26 23 23
+ >EAS54_6_R1_2_1_540_792
+ 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26
+ 26 18 26 23 18
+ >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
+
+ And a matching FASTA file::
+
+ >EAS54_6_R1_2_1_413_324
+ CCCTTCTTGTCTTCAGCGTTTCTCC
+ >EAS54_6_R1_2_1_540_792
+ TTGGCAGGCCAAGGCCGATGGATCA
+ >EAS54_6_R1_2_1_443_348
+ GTTGCTTCTGGCGTGGGTGGGGGGG
+
+ You can parse these separately using Bio.SeqIO with the "qual" and
+ "fasta" formats, but then you'll get a group of SeqRecord objects with
+ no sequence, and a matching group with the sequence but not the
+ qualities. Because it only deals with one input file handle, Bio.SeqIO
+ can't be used to read the two files together - but this function can!
+ For example,
+
+ >>> rec_iter = PairedFastaQualIterator(open("Quality/example.fasta", "rU"),
+ ... open("Quality/example.qual", "rU"))
+ >>> for record in rec_iter :
+ ... print 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
+
+ As with the FASTQ or QUAL parsers, if you want to look at the qualities,
+ they are in each record's per-letter-annotation dictionary as a simple
+ list of integers:
+
+ >>> 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]
+
+ If you have access to data as a FASTQ format file, using that directly
+ would be simpler and more straight forward. Note that you can easily use
+ this function to convert paired FASTA and QUAL files into FASTQ files:
+
+ >>> from Bio import SeqIO
+ >>> rec_iter = PairedFastaQualIterator(open("Quality/example.fasta", "rU"),
+ ... open("Quality/example.qual", "rU"))
+ >>> out_handle = open("Quality/temp.fastq", "w")
+ >>> SeqIO.write(rec_iter, out_handle, "fastq")
+ 3
+ >>> out_handle.close()
+
+ And don't forget to clean up the temp file if you don't need it anymore:
+
+ >>> import os
+ >>> os.remove("Quality/temp.fastq")
+ """
+ from Bio.SeqIO.FastaIO import FastaIterator
+ fasta_iter = FastaIterator(fasta_handle, alphabet=alphabet, \
+ title2ids=title2ids)
+ qual_iter = QualPhredIterator(qual_handle, alphabet=alphabet, \
+ title2ids=title2ids)
+
+ #Using zip(...) would create a list loading everything into memory!
+ #It would also not catch any extra records found in only one file.
+ while True :
+ try :
+ f_rec = fasta_iter.next()
+ except StopIteration :
+ f_rec = None
+ try :
+ q_rec = qual_iter.next()
+ except StopIteration :
+ q_rec = None
+ if f_rec is None and q_rec is None :
+ #End of both files
+ break
+ if f_rec is None :
+ raise ValueError("FASTA file has more entries than the QUAL file.")
+ if q_rec is None :
+ raise ValueError("QUAL file has more entries than the FASTA file.")
+ if f_rec.id != q_rec.id :
+ raise ValueError("FASTA and QUAL entries do not match (%s vs %s)." \
+ % (f_rec.id, q_rec.id))
+ if len(f_rec) != len(q_rec.letter_annotations["phred_quality"]) :
+ raise ValueError("Sequence length and number of quality scores disagree for %s" \
+ % f_rec.id)
+ #Merge the data....
+ f_rec.letter_annotations["phred_quality"] = q_rec.letter_annotations["phred_quality"]
+ yield f_rec
+ #Done
+
+
+def _test():
+ """Run the Bio.SeqIO module's doctests.
+
+ This will try and locate the unit tests directory, and run the doctests
+ from there in order that the relative paths used in the examples work.
+ """
+ import doctest
+ import os
+ if os.path.isdir(os.path.join("..","..","Tests")) :
+ print "Runing doctests..."
+ cur_dir = os.path.abspath(os.curdir)
+ os.chdir(os.path.join("..","..","Tests"))
+ assert os.path.isfile("Quality/example.fastq")
+ assert os.path.isfile("Quality/example.fasta")
+ assert os.path.isfile("Quality/example.qual")
+ assert os.path.isfile("Quality/tricky.fastq")
+ assert os.path.isfile("Quality/solexa.fastq")
+ doctest.testmod()
+ os.chdir(cur_dir)
+ del cur_dir
+ print "Done"
+
+if __name__ == "__main__" :
+ _test()
+
git://source.jalview.org / jabaws.git / blobdiff