Package Bio :: Package SeqIO :: Module QualityIO
[hide private]
[frames] | no frames]

Source Code for Module Bio.SeqIO.QualityIO

   1  # Copyright 2009 by Peter Cock.  All rights reserved. 
   2  # This code is part of the Biopython distribution and governed by its 
   3  # license.  Please see the LICENSE file that should have been included 
   4  # as part of this package. 
   5  # 
   6  # This module is for reading and writing FASTQ and QUAL format files as 
   7  # SeqRecord objects, and is expected to be used via the Bio.SeqIO API. 
   8   
   9  """Bio.SeqIO support for the "fastq" and "qual" file formats. 
  10   
  11  Note that you are expected to use this code via the Bio.SeqIO interface, as 
  12  shown below. 
  13   
  14  The FASTQ file format is used frequently at the Wellcome Trust Sanger Institute 
  15  to bundle a FASTA sequence and its PHRED quality data (integers between 0 and 
  16  90).  Rather than using a single FASTQ file, often paired FASTA and QUAL files 
  17  are used containing the sequence and the quality information separately. 
  18   
  19  The PHRED software reads DNA sequencing trace files, calls bases, and 
  20  assigns a quality value between 0 and 90 to each called base using a logged 
  21  transformation of the error probability, Q = -10 log10( Pe ), for example:: 
  22   
  23      Pe = 0.0,         Q =  0 
  24      Pe = 0.1,         Q = 10 
  25      Pe = 0.01,        Q = 20 
  26      ... 
  27      Pe = 0.00000001,  Q = 80 
  28      Pe = 0.000000001, Q = 90 
  29   
  30  In the QUAL format these quality values are held as space separated text in 
  31  a FASTA like file format.  In the FASTQ format, each quality values is encoded 
  32  with a single ASCI character using chr(Q+33), meaning zero maps to the 
  33  character "!" and for example 80 maps to "q".  The sequences and quality are 
  34  then stored in pairs in a FASTA like format. 
  35   
  36  Unfortunately there is no official document describing the FASTQ file format, 
  37  and worse, several related but different variants exist.  Reasonable 
  38  documentation exists at: http://maq.sourceforge.net/fastq.shtml 
  39   
  40  Solexa/Illumina quality scores use Q = - 10 log10 ( Pe / (1-Pe) ), which can 
  41  be negative or easily exceed 90.  PHRED scores and Solexa scores are NOT 
  42  interchangeable (but a reasonable mapping can be achieved between them). 
  43  Confusingly Solexa produces a FASTQ like file but using their own score 
  44  mapping instead. 
  45   
  46  Also note that Roche 454 sequencers can output files in the QUAL format, and 
  47  thankfully they use PHREP style scores like Sanger.  To extract QUAL files from 
  48  a Roche 454 SFF binary file, use the Roche off instrument command line tool 
  49  "sffinfo" with the -q or -qual argument.  You can extract a matching FASTA file 
  50  using the -s or -seq argument instead. 
  51   
  52  You are expected to use this module via the Bio.SeqIO functions, with the 
  53  following format names: 
  54   - "fastq" means Sanger style FASTQ files using PHRED scores. 
  55   - "fastq-solexa" means Solexa/Illumina style FASTQ files. 
  56   - "qual" means simple quality files using PHRED scores. 
  57   
  58  For example, consider the following short FASTQ file (extracted from a real 
  59  NCBI dataset):: 
  60   
  61      @EAS54_6_R1_2_1_413_324 
  62      CCCTTCTTGTCTTCAGCGTTTCTCC 
  63      + 
  64      ;;3;;;;;;;;;;;;7;;;;;;;88 
  65      @EAS54_6_R1_2_1_540_792 
  66      TTGGCAGGCCAAGGCCGATGGATCA 
  67      + 
  68      ;;;;;;;;;;;7;;;;;-;;;3;83 
  69      @EAS54_6_R1_2_1_443_348 
  70      GTTGCTTCTGGCGTGGGTGGGGGGG 
  71      + 
  72      ;;;;;;;;;;;9;7;;.7;393333 
  73   
  74  This contains three reads of length 25.  From the read length these were 
  75  probably originally from an early Solexa/Illumina sequencer but NCBI have 
  76  followed the Sanger FASTQ convention and this actually uses PHRED style 
  77  qualities.  This means we can parse this file using Bio.SeqIO using "fastq" 
  78  as the format name: 
  79   
  80      >>> from Bio import SeqIO 
  81      >>> for record in SeqIO.parse(open("Quality/example.fastq"), "fastq") : 
  82      ...     print record.id, record.seq 
  83      EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC 
  84      EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA 
  85      EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG 
  86   
  87  The qualities are held as a list of integers in each record's annotation: 
  88   
  89      >>> print record 
  90      ID: EAS54_6_R1_2_1_443_348 
  91      Name: EAS54_6_R1_2_1_443_348 
  92      Description: EAS54_6_R1_2_1_443_348 
  93      Number of features: 0 
  94      Per letter annotation for: phred_quality 
  95      Seq('GTTGCTTCTGGCGTGGGTGGGGGGG', SingleLetterAlphabet()) 
  96      >>> print record.letter_annotations["phred_quality"] 
  97      [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] 
  98   
  99  You can use the SeqRecord format method you can show this in the QUAL format: 
 100   
 101      >>> print record.format("qual") 
 102      >EAS54_6_R1_2_1_443_348 
 103      26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18 
 104      24 18 18 18 18 
 105      <BLANKLINE> 
 106   
 107  Or go back to the FASTQ format, 
 108   
 109      >>> print record.format("fastq") 
 110      @EAS54_6_R1_2_1_443_348 
 111      GTTGCTTCTGGCGTGGGTGGGGGGG 
 112      + 
 113      ;;;;;;;;;;;9;7;;.7;393333 
 114      <BLANKLINE> 
 115   
 116  You can also get Biopython to convert the scores and show a Solexa style 
 117  FASTQ file: 
 118   
 119      >>> print record.format("fastq-solexa") 
 120      @EAS54_6_R1_2_1_443_348 
 121      GTTGCTTCTGGCGTGGGTGGGGGGG 
 122      + 
 123      ZZZZZZZZZZZXZVZZMVZRXRRRR 
 124      <BLANKLINE> 
 125   
 126  If you wanted to trim your sequences (perhaps to remove low quality regions, 
 127  or to remove a primer sequence), try slicing the SeqRecord objects.  e.g. 
 128   
 129      >>> sub_rec = record[5:15] 
 130      >>> print sub_rec 
 131      ID: EAS54_6_R1_2_1_443_348 
 132      Name: EAS54_6_R1_2_1_443_348 
 133      Description: EAS54_6_R1_2_1_443_348 
 134      Number of features: 0 
 135      Per letter annotation for: phred_quality 
 136      Seq('TTCTGGCGTG', SingleLetterAlphabet()) 
 137      >>> print sub_rec.letter_annotations["phred_quality"] 
 138      [26, 26, 26, 26, 26, 26, 24, 26, 22, 26] 
 139      >>> print sub_rec.format("fastq") 
 140      @EAS54_6_R1_2_1_443_348 
 141      TTCTGGCGTG 
 142      + 
 143      ;;;;;;9;7; 
 144      <BLANKLINE> 
 145       
 146  If you wanted to, you could read in this FASTQ file, and save it as a QUAL file: 
 147   
 148      >>> from Bio import SeqIO 
 149      >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq") 
 150      >>> out_handle = open("Quality/temp.qual", "w") 
 151      >>> SeqIO.write(record_iterator, out_handle, "qual") 
 152      3 
 153      >>> out_handle.close() 
 154   
 155  You can of course read in a QUAL file, such as the one we just created: 
 156   
 157      >>> from Bio import SeqIO 
 158      >>> for record in SeqIO.parse(open("Quality/temp.qual"), "qual") : 
 159      ...     print record.id, record.seq 
 160      EAS54_6_R1_2_1_413_324 ????????????????????????? 
 161      EAS54_6_R1_2_1_540_792 ????????????????????????? 
 162      EAS54_6_R1_2_1_443_348 ????????????????????????? 
 163   
 164  Notice that QUAL files don't have a proper sequence present!  But the quality 
 165  information is there: 
 166   
 167      >>> print record 
 168      ID: EAS54_6_R1_2_1_443_348 
 169      Name: EAS54_6_R1_2_1_443_348 
 170      Description: EAS54_6_R1_2_1_443_348 
 171      Number of features: 0 
 172      Per letter annotation for: phred_quality 
 173      UnknownSeq(25, alphabet = SingleLetterAlphabet(), character = '?') 
 174      >>> print record.letter_annotations["phred_quality"] 
 175      [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] 
 176   
 177  Just to keep things tidy, if you are following this example yourself, you can 
 178  delete this temporary file now: 
 179   
 180      >>> import os 
 181      >>> os.remove("Quality/temp.qual") 
 182   
 183  Sometimes you won't have a FASTQ file, but rather just a pair of FASTA and QUAL 
 184  files.  Because the Bio.SeqIO system is designed for reading single files, you 
 185  would have to read the two in separately and then combine the data.  However, 
 186  since this is such a common thing to want to do, there is a helper iterator 
 187  defined in this module that does this for you - PairedFastaQualIterator. 
 188   
 189  Alternatively, if you have enough RAM to hold all the records in memory at once, 
 190  then a simple dictionary approach would work: 
 191   
 192      >>> from Bio import SeqIO 
 193      >>> reads = SeqIO.to_dict(SeqIO.parse(open("Quality/example.fasta"), "fasta")) 
 194      >>> for rec in SeqIO.parse(open("Quality/example.qual"), "qual") : 
 195      ...     reads[rec.id].letter_annotations["phred_quality"]=rec.letter_annotations["phred_quality"] 
 196   
 197  You can then access any record by its key, and get both the sequence and the 
 198  quality scores. 
 199   
 200      >>> print reads["EAS54_6_R1_2_1_540_792"].format("fastq") 
 201      @EAS54_6_R1_2_1_540_792 
 202      TTGGCAGGCCAAGGCCGATGGATCA 
 203      + 
 204      ;;;;;;;;;;;7;;;;;-;;;3;83 
 205      <BLANKLINE> 
 206   
 207  It is important that you explicitly tell Bio.SeqIO which FASTQ variant you are 
 208  using ("fastq" for the Sanger standard using PHRED values, or "fastq-solexa" 
 209  for the Solexa/Illumina variant), as this cannot be detected reliably 
 210  automatically. 
 211  """ 
 212  __docformat__ = "epytext en" #Don't just use plain text in epydoc API pages! 
 213   
 214  #See also http://blog.malde.org/index.php/2008/09/09/the-fastq-file-format-for-sequences/ 
 215   
 216  from Bio.Alphabet import single_letter_alphabet 
 217  from Bio.Seq import Seq, UnknownSeq 
 218  from Bio.SeqRecord import SeqRecord 
 219  from Interfaces import SequentialSequenceWriter 
 220  from math import log 
 221   
 222  # define score offsets. See discussion for differences between Sanger and 
 223  # Solexa offsets. 
 224  SANGER_SCORE_OFFSET = 33 
 225  SOLEXA_SCORE_OFFSET = 64 
 226   
227 -def solexa_quality_from_phred(phred_quality) :
228 """Covert a PHRED quality (range 0 to about 90) to a Solexa quality. 229 230 This will return a floating point number, it is up to you to round this to 231 the nearest integer if appropriate. e.g. 232 233 >>> print "%0.2f" % round(solexa_quality_from_phred(80),2) 234 80.00 235 >>> print "%0.2f" % round(solexa_quality_from_phred(50),2) 236 50.00 237 >>> print "%0.2f" % round(solexa_quality_from_phred(20),2) 238 19.96 239 >>> print "%0.2f" % round(solexa_quality_from_phred(10),2) 240 9.54 241 >>> print "%0.2f" % round(solexa_quality_from_phred(1),2) 242 -5.87 243 """ 244 return 10*log(10**(phred_quality/10.0) - 1, 10)
245
246 -def phred_quality_from_solexa(solexa_quality) :
247 """Convert a Solexa quality (which can be negative) to a PHRED quality. 248 249 This will return a floating point number, it is up to you to round this to 250 the nearest integer if appropriate. e.g. 251 252 >>> print "%0.2f" % round(phred_quality_from_solexa(80),2) 253 80.00 254 >>> print "%0.2f" % round(phred_quality_from_solexa(20),2) 255 20.04 256 >>> print "%0.2f" % round(phred_quality_from_solexa(10),2) 257 10.41 258 >>> print "%0.2f" % round(phred_quality_from_solexa(0),2) 259 3.01 260 >>> print "%0.2f" % round(phred_quality_from_solexa(-10),2) 261 0.41 262 """ 263 return 10*log(10**(solexa_quality/10.0) + 1, 10)
264
265 -def _get_phred_quality(record) :
266 """Extract PHRED qualities from a SeqRecord's letter_annotations (PRIVATE). 267 268 If there are no PHRED qualities, but there are Solexa qualities, those are 269 used instead after conversion. 270 """ 271 try : 272 return record.letter_annotations["phred_quality"] 273 except KeyError : 274 pass 275 try : 276 return [phred_quality_from_solexa(q) for \ 277 q in record.letter_annotations["solexa_quality"]] 278 except KeyError : 279 raise ValueError("No suitable quality scores found in letter_annotations " 280 "of SeqRecord (id=%s)." % record.id)
281
282 -def _get_solexa_quality(record) :
283 """Extract Solexa qualities from a SeqRecord's letter_annotations (PRIVATE). 284 285 If there are no Solexa qualities, but there are PHRED qualities, those are 286 used instead after conversion. 287 """ 288 try : 289 return record.letter_annotations["solexa_quality"] 290 except KeyError : 291 pass 292 try : 293 return [solexa_quality_from_phred(q) for \ 294 q in record.letter_annotations["phred_quality"]] 295 except KeyError : 296 raise ValueError("No suitable quality scores found in letter_annotation " 297 "of SeqRecord (id=%s)." % record.id)
298 299 300 #TODO - Default to nucleotide or even DNA?
301 -def FastqGeneralIterator(handle) :
302 """Iterate over Fastq records as string tuples (not as SeqRecord objects). 303 304 This code does not try to interpret the quality string numerically. It 305 just returns tuples of the title, sequence and quality as strings. For 306 the sequence and quality, any whitespace (such as new lines) is removed. 307 308 Our SeqRecord based FASTQ iterators call this function internally, and then 309 turn the strings into a SeqRecord objects, mapping the quality string into 310 a list of numerical scores. If you want to do a custom quality mapping, 311 then you might consider calling this function directly. 312 313 For parsing FASTQ files, the title string from the "@" line at the start 314 of each record can optionally be omitted on the "+" lines. If it is 315 repeated, it must be identical. 316 317 The sequence string and the quality string can optionally be split over 318 multiple lines, although several sources discourage this. In comparison, 319 for the FASTA file format line breaks between 60 and 80 characters are 320 the norm. 321 322 WARNING - Because the "@" character can appear in the quality string, 323 this can cause problems as this is also the marker for the start of 324 a new sequence. In fact, the "+" sign can also appear as well. Some 325 sources recommended having no line breaks in the quality to avoid this, 326 but even that is not enough, consider this example:: 327 328 @071113_EAS56_0053:1:1:998:236 329 TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA 330 +071113_EAS56_0053:1:1:998:236 331 IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III 332 @071113_EAS56_0053:1:1:182:712 333 ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG 334 + 335 @IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+ 336 @071113_EAS56_0053:1:1:153:10 337 TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT 338 + 339 IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6 340 @071113_EAS56_0053:1:3:990:501 341 TGGGAGGTTTTATGTGGA 342 AAGCAGCAATGTACAAGA 343 + 344 IIIIIII.IIIIII1@44 345 @-7.%<&+/$/%4(++(% 346 347 This is four PHRED encoded FASTQ entries originally from an NCBI source 348 (given the read length of 36, these are probably Solexa Illumna reads where 349 the quality has been mapped onto the PHRED values). 350 351 This example has been edited to illustrate some of the nasty things allowed 352 in the FASTQ format. Firstly, on the "+" lines most but not all of the 353 (redundant) identifiers are ommited. In real files it is likely that all or 354 none of these extra identifiers will be present. 355 356 Secondly, while the first three sequences have been shown without line 357 breaks, the last has been split over multiple lines. In real files any line 358 breaks are likely to be consistent. 359 360 Thirdly, some of the quality string lines start with an "@" character. For 361 the second record this is unavoidable. However for the fourth sequence this 362 only happens because its quality string is split over two lines. A naive 363 parser could wrongly treat any line starting with an "@" as the beginning of 364 a new sequence! This code copes with this possible ambiguity by keeping track 365 of the length of the sequence which gives the expected length of the quality 366 string. 367 368 Using this tricky example file as input, this short bit of code demonstrates 369 what this parsing function would return: 370 371 >>> handle = open("Quality/tricky.fastq", "rU") 372 >>> for (title, sequence, quality) in FastqGeneralIterator(handle) : 373 ... print title 374 ... print sequence, quality 375 071113_EAS56_0053:1:1:998:236 376 TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III 377 071113_EAS56_0053:1:1:182:712 378 ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG @IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+ 379 071113_EAS56_0053:1:1:153:10 380 TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6 381 071113_EAS56_0053:1:3:990:501 382 TGGGAGGTTTTATGTGGAAAGCAGCAATGTACAAGA IIIIIII.IIIIII1@44@-7.%<&+/$/%4(++(% 383 >>> handle.close() 384 385 Finally we note that some sources state that the quality string should 386 start with "!" (which using the PHRED mapping means the first letter always 387 has a quality score of zero). This rather restrictive rule is not widely 388 observed, so is therefore ignored here. One plus point about this "!" rule 389 is that (provided there are no line breaks in the quality sequence) it 390 would prevent the above problem with the "@" character. 391 """ 392 #Skip any text before the first record (e.g. blank lines, comments?) 393 while True : 394 line = handle.readline() 395 if line == "" : return #Premature end of file, or just empty? 396 if line[0] == "@" : 397 break 398 399 while True : 400 if line[0]!="@" : 401 raise ValueError("Records in Fastq files should start with '@' character") 402 title_line = line[1:].rstrip() 403 404 seq_lines = [] 405 line = handle.readline() 406 while True: 407 if not line : 408 raise ValueError("End of file without quality information.") 409 if line[0] == "+": 410 #The title here is optional, but if present must match! 411 if line[1:].rstrip() and line[1:].rstrip() != title_line : 412 raise ValueError("Sequence and quality captions differ.") 413 break 414 seq_lines.extend(line.split()) #removes any whitespace 415 line = handle.readline() 416 417 seq_string = "".join(seq_lines) 418 del seq_lines 419 420 quality_lines = [] 421 line = handle.readline() 422 while True: 423 if not line : break 424 if line[0] == "@": 425 #This COULD be the start of a new sequence. However, it MAY just 426 #be a line of quality data which starts with a "@" character. We 427 #should be able to check this by looking at the sequence length 428 #and the amount of quality data found so far. 429 if len("".join(quality_lines)) >= len(seq_string) : 430 #We expect it to be equal if this is the start of a new record. 431 #If the quality data is longer, we'll raise an error below. 432 break 433 #Continue - its just some (more) sequence data. 434 435 quality_lines.extend(line.split()) #removes any whitespace 436 line = handle.readline() 437 438 quality_string = "".join(quality_lines) 439 del quality_lines 440 441 if len(seq_string) != len(quality_string) : 442 raise ValueError("Lengths of sequence and quality values differs " 443 " for %s (%i and %i)." \ 444 % (title_line, len(seq_string), len(quality_string))) 445 446 #Return the record and then continue... 447 yield (title_line, seq_string, quality_string) 448 if not line : return #StopIteration at end of file 449 assert False, "Should not reach this line"
450 451 #This is a generator function!
452 -def FastqPhredIterator(handle, alphabet = single_letter_alphabet, title2ids = None) :
453 """Generator function to iterate over FASTQ records (as SeqRecord objects). 454 455 - handle - input file 456 - alphabet - optional alphabet 457 - title2ids - A function that, when given the title line from the FASTQ 458 file (without the beginning >), will return the id, name and 459 description (in that order) for the record as a tuple of 460 strings. If this is not given, then the entire title line 461 will be used as the description, and the first word as the 462 id and name. 463 464 Note that use of title2ids matches that of Bio.SeqIO.FastaIO. 465 466 For each sequence in a (Sanger style) FASTQ file there is a matching string 467 encoding the PHRED qualities (integers between 0 and about 90) using ASCII 468 values with an offset of 33. 469 470 For example, consider a file containing three short reads:: 471 472 @EAS54_6_R1_2_1_413_324 473 CCCTTCTTGTCTTCAGCGTTTCTCC 474 + 475 ;;3;;;;;;;;;;;;7;;;;;;;88 476 @EAS54_6_R1_2_1_540_792 477 TTGGCAGGCCAAGGCCGATGGATCA 478 + 479 ;;;;;;;;;;;7;;;;;-;;;3;83 480 @EAS54_6_R1_2_1_443_348 481 GTTGCTTCTGGCGTGGGTGGGGGGG 482 + 483 ;;;;;;;;;;;9;7;;.7;393333 484 485 For each sequence (e.g. "CCCTTCTTGTCTTCAGCGTTTCTCC") there is a matching 486 string encoding the PHRED qualities using a ASCI values with an offset of 487 33 (e.g. ";;3;;;;;;;;;;;;7;;;;;;;88"). 488 489 Using this module directly you might run: 490 491 >>> handle = open("Quality/example.fastq", "rU") 492 >>> for record in FastqPhredIterator(handle) : 493 ... print record.id, record.seq 494 EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC 495 EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA 496 EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG 497 >>> handle.close() 498 499 Typically however, you would call this via Bio.SeqIO instead with "fastq" as 500 the format: 501 502 >>> from Bio import SeqIO 503 >>> handle = open("Quality/example.fastq", "rU") 504 >>> for record in SeqIO.parse(handle, "fastq") : 505 ... print record.id, record.seq 506 EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC 507 EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA 508 EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG 509 >>> handle.close() 510 511 If you want to look at the qualities, they are record in each record's 512 per-letter-annotation dictionary as a simple list of integers: 513 514 >>> print record.letter_annotations["phred_quality"] 515 [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] 516 """ 517 for title_line, seq_string, quality_string in FastqGeneralIterator(handle) : 518 if title2ids : 519 id, name, descr = title2ids(title_line) 520 else : 521 descr = title_line 522 id = descr.split()[0] 523 name = id 524 record = SeqRecord(Seq(seq_string, alphabet), 525 id=id, name=name, description=descr) 526 527 assert SANGER_SCORE_OFFSET == ord("!") 528 #According to BioPerl documentation at least, the first character should 529 #be an "!" (and therefore quality zero). This seems crazy - what if the 530 #sequence has been trimmed to remove any poor quality sequence? In any 531 #case real examples from the NCBI don't follow this practice, so we 532 #won't enforce it here. 533 #e.g. ftp://ftp.ncbi.nih.gov/pub/TraceDB/ShortRead/SRA000271/fastq/200x36x36-071113_EAS56_0053-s_1_1.fastq.gz 534 # 535 #if quality_string[0] != "!" : 536 # raise ValueError("The quality string should always start with a ! character.") 537 qualities = [ord(letter)-SANGER_SCORE_OFFSET for letter in quality_string] 538 if qualities : 539 if min(qualities) < 0 or max(qualities) > 90 : 540 raise ValueError("Quality score outside 0 to 90 found - these are perhaps " 541 "in a Solexa/Illumina format, not the Sanger FASTQ format " 542 "which uses PHRED scores.") 543 record.letter_annotations["phred_quality"] = qualities 544 yield record
545 546 #This is a generator function!
547 -def FastqSolexaIterator(handle, alphabet = single_letter_alphabet, title2ids = None) :
548 """Parsing the Solexa/Illumina FASTQ like files (which differ in the quality mapping). 549 550 The optional arguments are the same as those for the FastqPhredIterator. 551 552 For each sequence in Solexa/Illumina FASTQ files there is a matching string 553 encoding the Solexa integer qualities using ASCII values with an offset 554 of 64. Solexa scores are scaled differently to PHRED scores, and Biopython 555 will NOT perform any automatic conversion when loading. 556 557 For example, consider a file containing these five records:: 558 559 @SLXA-B3_649_FC8437_R1_1_1_610_79 560 GATGTGCAATACCTTTGTAGAGGAA 561 +SLXA-B3_649_FC8437_R1_1_1_610_79 562 YYYYYYYYYYYYYYYYYYWYWYYSU 563 @SLXA-B3_649_FC8437_R1_1_1_397_389 564 GGTTTGAGAAAGAGAAATGAGATAA 565 +SLXA-B3_649_FC8437_R1_1_1_397_389 566 YYYYYYYYYWYYYYWWYYYWYWYWW 567 @SLXA-B3_649_FC8437_R1_1_1_850_123 568 GAGGGTGTTGATCATGATGATGGCG 569 +SLXA-B3_649_FC8437_R1_1_1_850_123 570 YYYYYYYYYYYYYWYYWYYSYYYSY 571 @SLXA-B3_649_FC8437_R1_1_1_362_549 572 GGAAACAAAGTTTTTCTCAACATAG 573 +SLXA-B3_649_FC8437_R1_1_1_362_549 574 YYYYYYYYYYYYYYYYYYWWWWYWY 575 @SLXA-B3_649_FC8437_R1_1_1_183_714 576 GTATTATTTAATGGCATACACTCAA 577 +SLXA-B3_649_FC8437_R1_1_1_183_714 578 YYYYYYYYYYWYYYYWYWWUWWWQQ 579 580 Using this module directly you might run: 581 582 >>> handle = open("Quality/solexa_example.fastq", "rU") 583 >>> for record in FastqSolexaIterator(handle) : 584 ... print record.id, record.seq 585 SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA 586 SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA 587 SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG 588 SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG 589 SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA 590 >>> handle.close() 591 592 Typically however, you would call this via Bio.SeqIO instead with "fastq" as 593 the format: 594 595 >>> from Bio import SeqIO 596 >>> handle = open("Quality/solexa_example.fastq", "rU") 597 >>> for record in SeqIO.parse(handle, "fastq-solexa") : 598 ... print record.id, record.seq 599 SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA 600 SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA 601 SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG 602 SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG 603 SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA 604 >>> handle.close() 605 606 If you want to look at the qualities, they are recorded in each record's 607 per-letter-annotation dictionary as a simple list of integers: 608 609 >>> print record.letter_annotations["solexa_quality"] 610 [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] 611 612 These scores aren't very good, but they are high enough that they map 613 almost exactly onto PHRED scores: 614 615 >>> print "%0.2f" % phred_quality_from_solexa(25) 616 25.01 617 618 Let's look at another example read which is even worse, where there are 619 more noticeable differences between the Solexa and PHRED scores:: 620 621 @slxa_0013_1_0001_24 622 ACAAAAATCACAAGCATTCTTATACACC 623 +slxa_0013_1_0001_24 624 ??????????????????:??<?<-6%. 625 626 Again, you would typically use Bio.SeqIO to read this file in (rather than 627 calling the Bio.SeqIO.QualtityIO module directly). Most FASTQ files will 628 contain thousands of reads, so you would normally use Bio.SeqIO.parse() 629 as shown above. This example has only as one entry, so instead we can 630 use the Bio.SeqIO.read() function: 631 632 >>> from Bio import SeqIO 633 >>> handle = open("Quality/solexa.fastq", "rU") 634 >>> record = SeqIO.read(handle, "fastq-solexa") 635 >>> handle.close() 636 >>> print record.id, record.seq 637 slxa_0013_1_0001_24 ACAAAAATCACAAGCATTCTTATACACC 638 >>> print record.letter_annotations["solexa_quality"] 639 [-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] 640 641 These quality scores are so low that when converted from the Solexa scheme 642 into PHRED scores they look quite different: 643 644 >>> print "%0.2f" % phred_quality_from_solexa(-1) 645 2.54 646 647 Note you can use the Bio.SeqIO.write() function or the SeqRecord's format 648 method to output the record(s): 649 650 >>> print record.format("fastq-solexa") 651 @slxa_0013_1_0001_24 652 ACAAAAATCACAAGCATTCTTATACACC 653 + 654 ??????????????????:??<?<-6%. 655 <BLANKLINE> 656 657 Note this output is slightly different from the input file as Biopython 658 has left out the optional repetition of the sequence identifier on the "+" 659 line. If you want the to use PHRED scores, use "fastq" or "qual" as the 660 output format instead, and Biopython will do the conversion for you: 661 662 >>> print record.format("fastq") 663 @slxa_0013_1_0001_24 664 ACAAAAATCACAAGCATTCTTATACACC 665 + 666 $$$$$$$$$$$$$$$$$$"$$"$"!!!! 667 <BLANKLINE> 668 669 >>> print record.format("qual") 670 >slxa_0013_1_0001_24 671 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 672 <BLANKLINE> 673 """ 674 for title_line, seq_string, quality_string in FastqGeneralIterator(handle) : 675 if title2ids : 676 id, name, descr = title_line 677 else : 678 descr = title_line 679 id = descr.split()[0] 680 name = id 681 record = SeqRecord(Seq(seq_string, alphabet), 682 id=id, name=name, description=descr) 683 qualities = [ord(letter)-SOLEXA_SCORE_OFFSET for letter in quality_string] 684 #DO NOT convert these into PHRED qualities automatically! 685 record.letter_annotations["solexa_quality"] = qualities 686 yield record
687
688 -def QualPhredIterator(handle, alphabet = single_letter_alphabet, title2ids = None) :
689 """For QUAL files which include PHRED quality scores, but no sequence. 690 691 For example, consider this short QUAL file:: 692 693 >EAS54_6_R1_2_1_413_324 694 26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 695 26 26 26 23 23 696 >EAS54_6_R1_2_1_540_792 697 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26 698 26 18 26 23 18 699 >EAS54_6_R1_2_1_443_348 700 26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18 701 24 18 18 18 18 702 703 Using this module directly you might run: 704 705 >>> handle = open("Quality/example.qual", "rU") 706 >>> for record in QualPhredIterator(handle) : 707 ... print record.id, record.seq 708 EAS54_6_R1_2_1_413_324 ????????????????????????? 709 EAS54_6_R1_2_1_540_792 ????????????????????????? 710 EAS54_6_R1_2_1_443_348 ????????????????????????? 711 >>> handle.close() 712 713 Typically however, you would call this via Bio.SeqIO instead with "qual" 714 as the format: 715 716 >>> from Bio import SeqIO 717 >>> handle = open("Quality/example.qual", "rU") 718 >>> for record in SeqIO.parse(handle, "qual") : 719 ... print record.id, record.seq 720 EAS54_6_R1_2_1_413_324 ????????????????????????? 721 EAS54_6_R1_2_1_540_792 ????????????????????????? 722 EAS54_6_R1_2_1_443_348 ????????????????????????? 723 >>> handle.close() 724 725 Becase QUAL files don't contain the sequence string itself, the seq 726 property is set to an UnknownSeq object. As no alphabet was given, this 727 has defaulted to a generic single letter alphabet and the character "?" 728 used. 729 730 By specifying a nucleotide alphabet, "N" is used instead: 731 732 >>> from Bio import SeqIO 733 >>> from Bio.Alphabet import generic_dna 734 >>> handle = open("Quality/example.qual", "rU") 735 >>> for record in SeqIO.parse(handle, "qual", alphabet=generic_dna) : 736 ... print record.id, record.seq 737 EAS54_6_R1_2_1_413_324 NNNNNNNNNNNNNNNNNNNNNNNNN 738 EAS54_6_R1_2_1_540_792 NNNNNNNNNNNNNNNNNNNNNNNNN 739 EAS54_6_R1_2_1_443_348 NNNNNNNNNNNNNNNNNNNNNNNNN 740 >>> handle.close() 741 742 However, the quality scores themselves are available as a list of integers 743 in each record's per-letter-annotation: 744 745 >>> print record.letter_annotations["phred_quality"] 746 [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] 747 748 You can still slice one of these SeqRecord objects with an UnknownSeq: 749 750 >>> sub_record = record[5:10] 751 >>> print sub_record.id, sub_record.letter_annotations["phred_quality"] 752 EAS54_6_R1_2_1_443_348 [26, 26, 26, 26, 26] 753 """ 754 #Skip any text before the first record (e.g. blank lines, comments) 755 while True : 756 line = handle.readline() 757 if line == "" : return #Premature end of file, or just empty? 758 if line[0] == ">" : 759 break 760 761 while True : 762 if line[0]!=">" : 763 raise ValueError("Records in Fasta files should start with '>' character") 764 if title2ids : 765 id, name, descr = title2ids(line[1:].rstrip()) 766 else : 767 descr = line[1:].rstrip() 768 id = descr.split()[0] 769 name = id 770 771 qualities = [] 772 line = handle.readline() 773 while True: 774 if not line : break 775 if line[0] == ">": break 776 qualities.extend([int(word) for word in line.split()]) 777 line = handle.readline() 778 779 if qualities : 780 if min(qualities) < 0 or max(qualities) > 90 : 781 raise ValueError(("Quality score range for %s is %i to %i, outside the " \ 782 +"expected 0 to 90. Perhaps these are Solexa/Illumina " \ 783 +"scores, and not PHRED scores?") \ 784 % (id, min(qualities), max(qualities))) 785 786 #Return the record and then continue... 787 record = SeqRecord(UnknownSeq(len(qualities), alphabet), 788 id = id, name = name, description = descr) 789 record.letter_annotations["phred_quality"] = qualities 790 yield record 791 792 if not line : return #StopIteration 793 assert False, "Should not reach this line"
794
795 -class FastqPhredWriter(SequentialSequenceWriter):
796 """Class to write FASTQ format files (using PHRED quality scores). 797 798 Although you can use this class directly, you are strongly encouraged 799 to use the Bio.SeqIO.write() function instead. For example, this code 800 reads in a FASTQ (PHRED) file and re-saves it as another FASTQ (PHRED) 801 file: 802 803 >>> from Bio import SeqIO 804 >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq") 805 >>> out_handle = open("Quality/temp.fastq", "w") 806 >>> SeqIO.write(record_iterator, out_handle, "fastq") 807 3 808 >>> out_handle.close() 809 810 You might want to do this if the original file included extra line breaks, 811 which while valid may not be supported by all tools. The output file from 812 Biopython will have each sequence on a single line, and each quality 813 string on a single line (which is considered desirable for maximum 814 compatibility). 815 816 In this next example, a Solexa FASTQ file is converted into a standard 817 Sanger style FASTQ file using PHRED qualities: 818 819 >>> from Bio import SeqIO 820 >>> record_iterator = SeqIO.parse(open("Quality/solexa.fastq"), "fastq-solexa") 821 >>> out_handle = open("Quality/temp.fastq", "w") 822 >>> SeqIO.write(record_iterator, out_handle, "fastq") 823 1 824 >>> out_handle.close() 825 826 This code is also called if you use the .format("fastq") method of a 827 SeqRecord. 828 829 P.S. To avoid cluttering up your working directory, you can delete this 830 temporary file now: 831 832 >>> import os 833 >>> os.remove("Quality/temp.fastq") 834 835 """
836 - def write_record(self, record):
837 """Write a single FASTQ record to the file.""" 838 assert self._header_written 839 assert not self._footer_written 840 self._record_written = True 841 842 #TODO - Is an empty sequence allowed in FASTQ format? 843 assert SANGER_SCORE_OFFSET == ord("!") 844 #This rounds to the nearest integer: 845 qualities = "".join([chr(int(round(q+SANGER_SCORE_OFFSET,0))) for q \ 846 in _get_phred_quality(record)]) 847 if record.seq is None: 848 raise ValueError("No sequence for record %s" % record.id) 849 if len(qualities) != len(record) : 850 raise ValueError("Record %s has sequence length %i but %i quality scores" \ 851 % (record.id, len(record), len(qualities))) 852 853 title = self.clean(record.id) #TODO - add the description too? cf Fasta output 854 self.handle.write("@%s\n%s\n+\n%s\n" % (title, record.seq, qualities))
855
856 -class QualPhredWriter(SequentialSequenceWriter):
857 """Class to write QUAL format files (using PHRED quality scores). 858 859 Although you can use this class directly, you are strongly encouraged 860 to use the Bio.SeqIO.write() function instead. For example, this code 861 reads in a FASTQ file and saves the quality scores into a QUAL file: 862 863 >>> from Bio import SeqIO 864 >>> record_iterator = SeqIO.parse(open("Quality/example.fastq"), "fastq") 865 >>> out_handle = open("Quality/temp.qual", "w") 866 >>> SeqIO.write(record_iterator, out_handle, "qual") 867 3 868 >>> out_handle.close() 869 870 This code is also called if you use the .format("qual") method of a 871 SeqRecord. 872 873 P.S. Don't forget to clean up the temp file if you don't need it anymore: 874 875 >>> import os 876 >>> os.remove("Quality/temp.qual") 877 """
878 - def __init__(self, handle, wrap=60, record2title=None):
879 """Create a QUAL writer. 880 881 Arguments: 882 - handle - Handle to an output file, e.g. as returned 883 by open(filename, "w") 884 - wrap - Optional line length used to wrap sequence lines. 885 Defaults to wrapping the sequence at 60 characters 886 Use zero (or None) for no wrapping, giving a single 887 long line for the sequence. 888 - record2title - Optional function to return the text to be 889 used for the title line of each record. By default 890 a combination of the record.id and record.description 891 is used. If the record.description starts with the 892 record.id, then just the record.description is used. 893 894 The record2title argument is present for consistency with the 895 Bio.SeqIO.FastaIO writer class. 896 """ 897 SequentialSequenceWriter.__init__(self, handle) 898 #self.handle = handle 899 self.wrap = None 900 if wrap : 901 if wrap < 1 : 902 raise ValueError 903 self.wrap = wrap 904 self.record2title = record2title
905
906 - def write_record(self, record):
907 """Write a single QUAL record to the file.""" 908 assert self._header_written 909 assert not self._footer_written 910 self._record_written = True 911 912 if self.record2title : 913 title=self.clean(self.record2title(record)) 914 else : 915 id = self.clean(record.id) 916 description = self.clean(record.description) 917 918 #if description[:len(id)]==id : 919 if description and description.split(None,1)[0]==id : 920 #The description includes the id at the start 921 title = description 922 else : 923 title = "%s %s" % (id, description) 924 925 assert "\n" not in title 926 assert "\r" not in title 927 self.handle.write(">%s\n" % title) 928 929 #This rounds to the nearest integer. 930 #TODO - can we put a float in a qual file? 931 qualities = [("%i" % round(q,0)) for q in _get_phred_quality(record)] 932 933 if self.wrap : 934 while qualities : 935 line=qualities.pop(0) 936 while qualities \ 937 and len(line) + 1 + len(qualities[0]) < self.wrap : 938 line += " " + qualities.pop(0) 939 self.handle.write(line + "\n") 940 else : 941 data = " ".join(qualities) 942 self.handle.write(data + "\n")
943
944 -class FastqSolexaWriter(SequentialSequenceWriter):
945 """Class to write FASTQ format files (using Solexa quality scores). 946 947 Although you can use this class directly, you are strongly encouraged 948 to use the Bio.SeqIO.write() function instead. For example, this code 949 reads in a FASTQ file and re-saves it as another FASTQ file: 950 951 >>> from Bio import SeqIO 952 >>> record_iterator = SeqIO.parse(open("Quality/solexa.fastq"), "fastq-solexa") 953 >>> out_handle = open("Quality/temp.fastq", "w") 954 >>> SeqIO.write(record_iterator, out_handle, "fastq-solexa") 955 1 956 >>> out_handle.close() 957 958 You might want to do this if the original file included extra line 959 breaks, which (while valid) may not be supported by all tools. The 960 output file from Biopython will have each sequence on a single line, and 961 each quality string on a single line (which is considered desirable for 962 maximum compatibility). 963 964 This code is also called if you use the .format("fastq-solexa") method of 965 a SeqRecord. 966 967 P.S. Don't forget to delete the temp file if you don't need it anymore: 968 969 >>> import os 970 >>> os.remove("Quality/temp.fastq") 971 """
972 - def write_record(self, record):
973 """Write a single FASTQ record to the file.""" 974 assert self._header_written 975 assert not self._footer_written 976 self._record_written = True 977 978 #TODO - Is an empty sequence allowed in FASTQ format? 979 qualities = "".join([chr(int(round(q+SOLEXA_SCORE_OFFSET,0))) for q \ 980 in _get_solexa_quality(record)]) 981 if record.seq is None: 982 raise ValueError("No sequence for record %s" % record.id) 983 if len(qualities) != len(record) : 984 raise ValueError("Record %s has sequence length %i but %i quality scores" \ 985 % (record.id, len(record), len(qualities))) 986 987 title = self.clean(record.id) #TODO - add the description too? cf Fasta output 988 self.handle.write("@%s\n%s\n+\n%s\n" % (title, record.seq, qualities))
989
990 -def PairedFastaQualIterator(fasta_handle, qual_handle, alphabet = single_letter_alphabet, title2ids = None) :
991 """Iterate over matched FASTA and QUAL files as SeqRecord objects. 992 993 For example, consider this short QUAL file:: 994 995 >EAS54_6_R1_2_1_413_324 996 26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 997 26 26 26 23 23 998 >EAS54_6_R1_2_1_540_792 999 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26 1000 26 18 26 23 18 1001 >EAS54_6_R1_2_1_443_348 1002 26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18 1003 24 18 18 18 18 1004 1005 And a matching FASTA file:: 1006 1007 >EAS54_6_R1_2_1_413_324 1008 CCCTTCTTGTCTTCAGCGTTTCTCC 1009 >EAS54_6_R1_2_1_540_792 1010 TTGGCAGGCCAAGGCCGATGGATCA 1011 >EAS54_6_R1_2_1_443_348 1012 GTTGCTTCTGGCGTGGGTGGGGGGG 1013 1014 You can parse these separately using Bio.SeqIO with the "qual" and 1015 "fasta" formats, but then you'll get a group of SeqRecord objects with 1016 no sequence, and a matching group with the sequence but not the 1017 qualities. Because it only deals with one input file handle, Bio.SeqIO 1018 can't be used to read the two files together - but this function can! 1019 For example, 1020 1021 >>> rec_iter = PairedFastaQualIterator(open("Quality/example.fasta", "rU"), 1022 ... open("Quality/example.qual", "rU")) 1023 >>> for record in rec_iter : 1024 ... print record.id, record.seq 1025 EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC 1026 EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA 1027 EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG 1028 1029 As with the FASTQ or QUAL parsers, if you want to look at the qualities, 1030 they are in each record's per-letter-annotation dictionary as a simple 1031 list of integers: 1032 1033 >>> print record.letter_annotations["phred_quality"] 1034 [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] 1035 1036 If you have access to data as a FASTQ format file, using that directly 1037 would be simpler and more straight forward. Note that you can easily use 1038 this function to convert paired FASTA and QUAL files into FASTQ files: 1039 1040 >>> from Bio import SeqIO 1041 >>> rec_iter = PairedFastaQualIterator(open("Quality/example.fasta", "rU"), 1042 ... open("Quality/example.qual", "rU")) 1043 >>> out_handle = open("Quality/temp.fastq", "w") 1044 >>> SeqIO.write(rec_iter, out_handle, "fastq") 1045 3 1046 >>> out_handle.close() 1047 1048 And don't forget to clean up the temp file if you don't need it anymore: 1049 1050 >>> import os 1051 >>> os.remove("Quality/temp.fastq") 1052 """ 1053 from Bio.SeqIO.FastaIO import FastaIterator 1054 fasta_iter = FastaIterator(fasta_handle, alphabet=alphabet, \ 1055 title2ids=title2ids) 1056 qual_iter = QualPhredIterator(qual_handle, alphabet=alphabet, \ 1057 title2ids=title2ids) 1058 1059 #Using zip(...) would create a list loading everything into memory! 1060 #It would also not catch any extra records found in only one file. 1061 while True : 1062 try : 1063 f_rec = fasta_iter.next() 1064 except StopIteration : 1065 f_rec = None 1066 try : 1067 q_rec = qual_iter.next() 1068 except StopIteration : 1069 q_rec = None 1070 if f_rec is None and q_rec is None : 1071 #End of both files 1072 break 1073 if f_rec is None : 1074 raise ValueError("FASTA file has more entries than the QUAL file.") 1075 if q_rec is None : 1076 raise ValueError("QUAL file has more entries than the FASTA file.") 1077 if f_rec.id != q_rec.id : 1078 raise ValueError("FASTA and QUAL entries do not match (%s vs %s)." \ 1079 % (f_rec.id, q_rec.id)) 1080 if len(f_rec) != len(q_rec.letter_annotations["phred_quality"]) : 1081 raise ValueError("Sequence length and number of quality scores disagree for %s" \ 1082 % f_rec.id) 1083 #Merge the data.... 1084 f_rec.letter_annotations["phred_quality"] = q_rec.letter_annotations["phred_quality"] 1085 yield f_rec
1086 #Done 1087 1088
1089 -def _test():
1090 """Run the Bio.SeqIO module's doctests. 1091 1092 This will try and locate the unit tests directory, and run the doctests 1093 from there in order that the relative paths used in the examples work. 1094 """ 1095 import doctest 1096 import os 1097 if os.path.isdir(os.path.join("..","..","Tests")) : 1098 print "Runing doctests..." 1099 cur_dir = os.path.abspath(os.curdir) 1100 os.chdir(os.path.join("..","..","Tests")) 1101 assert os.path.isfile("Quality/example.fastq") 1102 assert os.path.isfile("Quality/example.fasta") 1103 assert os.path.isfile("Quality/example.qual") 1104 assert os.path.isfile("Quality/tricky.fastq") 1105 assert os.path.isfile("Quality/solexa.fastq") 1106 doctest.testmod() 1107 os.chdir(cur_dir) 1108 del cur_dir 1109 print "Done"
1110 1111 if __name__ == "__main__" : 1112 _test() 1113