Package Bio :: Package GenBank :: Module Scanner
[hide private]
[frames] | no frames]

Source Code for Module Bio.GenBank.Scanner

   1  # Copyright 2007-2016 by Peter Cock.  All rights reserved. 
   2  # Revisions copyright 2010 by Uri Laserson.  All rights reserved. 
   3  # This code is part of the Biopython distribution and governed by its 
   4  # license.  Please see the LICENSE file that should have been included 
   5  # as part of this package. 
   6  """Internal code for parsing GenBank and EMBL files (PRIVATE). 
   7   
   8  This code is NOT intended for direct use.  It provides a basic scanner 
   9  (for use with a event consumer such as Bio.GenBank._FeatureConsumer) 
  10  to parse a GenBank or EMBL file (with their shared INSDC feature table). 
  11   
  12  It is used by Bio.GenBank to parse GenBank files 
  13  It is also used by Bio.SeqIO to parse GenBank and EMBL files 
  14   
  15  Feature Table Documentation: 
  16  http://www.insdc.org/files/feature_table.html 
  17  http://www.ncbi.nlm.nih.gov/projects/collab/FT/index.html 
  18  ftp://ftp.ncbi.nih.gov/genbank/docs/ 
  19  """ 
  20  # 17-MAR-2009: added wgs, wgs_scafld for GenBank whole genome shotgun master records. 
  21  # These are GenBank files that summarize the content of a project, and provide lists of 
  22  # scaffold and contig files in the project. These will be in annotations['wgs'] and 
  23  # annotations['wgs_scafld']. These GenBank files do not have sequences. See 
  24  # http://groups.google.com/group/bionet.molbio.genbank/browse_thread/thread/51fb88bf39e7dc36 
  25  # http://is.gd/nNgk 
  26  # for more details of this format, and an example. 
  27  # Added by Ying Huang & Iddo Friedberg 
  28   
  29  from __future__ import print_function 
  30   
  31  import warnings 
  32  import re 
  33  from collections import OrderedDict 
  34  from Bio.Seq import Seq 
  35  from Bio.SeqRecord import SeqRecord 
  36  from Bio.Alphabet import generic_protein 
  37  from Bio import BiopythonParserWarning 
38 39 40 -class InsdcScanner(object):
41 """Basic functions for breaking up a GenBank/EMBL file into sub sections. 42 43 The International Nucleotide Sequence Database Collaboration (INSDC) 44 between the DDBJ, EMBL, and GenBank. These organisations all use the 45 same "Feature Table" layout in their plain text flat file formats. 46 47 However, the header and sequence sections of an EMBL file are very 48 different in layout to those produced by GenBank/DDBJ.""" 49 50 # These constants get redefined with sensible values in the sub classes: 51 RECORD_START = "XXX" # "LOCUS " or "ID " 52 HEADER_WIDTH = 3 # 12 or 5 53 FEATURE_START_MARKERS = ["XXX***FEATURES***XXX"] 54 FEATURE_END_MARKERS = ["XXX***END FEATURES***XXX"] 55 FEATURE_QUALIFIER_INDENT = 0 56 FEATURE_QUALIFIER_SPACER = "" 57 SEQUENCE_HEADERS = ["XXX"] # with right hand side spaces removed 58
59 - def __init__(self, debug=0):
60 assert len(self.RECORD_START) == self.HEADER_WIDTH 61 for marker in self.SEQUENCE_HEADERS: 62 assert marker == marker.rstrip() 63 assert len(self.FEATURE_QUALIFIER_SPACER) == self.FEATURE_QUALIFIER_INDENT 64 self.debug = debug 65 self.handle = None 66 self.line = None
67
68 - def set_handle(self, handle):
69 self.handle = handle 70 self.line = ""
71
72 - def find_start(self):
73 """Read in lines until find the ID/LOCUS line, which is returned. 74 75 Any preamble (such as the header used by the NCBI on ``*.seq.gz`` archives) 76 will we ignored.""" 77 while True: 78 if self.line: 79 line = self.line 80 self.line = "" 81 else: 82 line = self.handle.readline() 83 if not line: 84 if self.debug: 85 print("End of file") 86 return None 87 if line[:self.HEADER_WIDTH] == self.RECORD_START: 88 if self.debug > 1: 89 print("Found the start of a record:\n" + line) 90 break 91 line = line.rstrip() 92 if line == "//": 93 if self.debug > 1: 94 print("Skipping // marking end of last record") 95 elif line == "": 96 if self.debug > 1: 97 print("Skipping blank line before record") 98 else: 99 # Ignore any header before the first ID/LOCUS line. 100 if self.debug > 1: 101 print("Skipping header line before record:\n" + line) 102 self.line = line 103 return line
104
105 - def parse_header(self):
106 """Return list of strings making up the header 107 108 New line characters are removed. 109 110 Assumes you have just read in the ID/LOCUS line. 111 """ 112 assert self.line[:self.HEADER_WIDTH] == self.RECORD_START, \ 113 "Not at start of record" 114 115 header_lines = [] 116 while True: 117 line = self.handle.readline() 118 if not line: 119 raise ValueError("Premature end of line during sequence data") 120 line = line.rstrip() 121 if line in self.FEATURE_START_MARKERS: 122 if self.debug: 123 print("Found feature table") 124 break 125 # if line[:self.HEADER_WIDTH]==self.FEATURE_START_MARKER[:self.HEADER_WIDTH]: 126 # if self.debug : print("Found header table (?)") 127 # break 128 if line[:self.HEADER_WIDTH].rstrip() in self.SEQUENCE_HEADERS: 129 if self.debug: 130 print("Found start of sequence") 131 break 132 if line == "//": 133 raise ValueError("Premature end of sequence data marker '//' found") 134 header_lines.append(line) 135 self.line = line 136 return header_lines
137
138 - def parse_features(self, skip=False):
139 """Return list of tuples for the features (if present) 140 141 Each feature is returned as a tuple (key, location, qualifiers) 142 where key and location are strings (e.g. "CDS" and 143 "complement(join(490883..490885,1..879))") while qualifiers 144 is a list of two string tuples (feature qualifier keys and values). 145 146 Assumes you have already read to the start of the features table. 147 """ 148 if self.line.rstrip() not in self.FEATURE_START_MARKERS: 149 if self.debug: 150 print("Didn't find any feature table") 151 return [] 152 153 while self.line.rstrip() in self.FEATURE_START_MARKERS: 154 self.line = self.handle.readline() 155 156 features = [] 157 line = self.line 158 while True: 159 if not line: 160 raise ValueError("Premature end of line during features table") 161 if line[:self.HEADER_WIDTH].rstrip() in self.SEQUENCE_HEADERS: 162 if self.debug: 163 print("Found start of sequence") 164 break 165 line = line.rstrip() 166 if line == "//": 167 raise ValueError("Premature end of features table, marker '//' found") 168 if line in self.FEATURE_END_MARKERS: 169 if self.debug: 170 print("Found end of features") 171 line = self.handle.readline() 172 break 173 if line[2:self.FEATURE_QUALIFIER_INDENT].strip() == "": 174 # This is an empty feature line between qualifiers. Empty 175 # feature lines within qualifiers are handled below (ignored). 176 line = self.handle.readline() 177 continue 178 if len(line) < self.FEATURE_QUALIFIER_INDENT: 179 warnings.warn("line too short to contain a feature: %r" % line, 180 BiopythonParserWarning) 181 line = self.handle.readline() 182 continue 183 184 if skip: 185 line = self.handle.readline() 186 while line[:self.FEATURE_QUALIFIER_INDENT] == self.FEATURE_QUALIFIER_SPACER: 187 line = self.handle.readline() 188 else: 189 # Build up a list of the lines making up this feature: 190 if line[self.FEATURE_QUALIFIER_INDENT] != " " \ 191 and " " in line[self.FEATURE_QUALIFIER_INDENT:]: 192 # The feature table design enforces a length limit on the feature keys. 193 # Some third party files (e.g. IGMT's EMBL like files) solve this by 194 # over indenting the location and qualifiers. 195 feature_key, line = line[2:].strip().split(None, 1) 196 feature_lines = [line] 197 warnings.warn("Overindented %s feature?" % feature_key, 198 BiopythonParserWarning) 199 else: 200 feature_key = line[2:self.FEATURE_QUALIFIER_INDENT].strip() 201 feature_lines = [line[self.FEATURE_QUALIFIER_INDENT:]] 202 line = self.handle.readline() 203 while line[:self.FEATURE_QUALIFIER_INDENT] == self.FEATURE_QUALIFIER_SPACER \ 204 or (line != '' and line.rstrip() == ""): # cope with blank lines in the midst of a feature 205 # Use strip to remove any harmless trailing white space AND and leading 206 # white space (e.g. out of spec files with too much indentation) 207 feature_lines.append(line[self.FEATURE_QUALIFIER_INDENT:].strip()) 208 line = self.handle.readline() 209 features.append(self.parse_feature(feature_key, feature_lines)) 210 self.line = line 211 return features
212
213 - def parse_feature(self, feature_key, lines):
214 r"""Expects a feature as a list of strings, returns a tuple (key, location, qualifiers) 215 216 For example given this GenBank feature:: 217 218 CDS complement(join(490883..490885,1..879)) 219 /locus_tag="NEQ001" 220 /note="conserved hypothetical [Methanococcus jannaschii]; 221 COG1583:Uncharacterized ACR; IPR001472:Bipartite nuclear 222 localization signal; IPR002743: Protein of unknown 223 function DUF57" 224 /codon_start=1 225 /transl_table=11 226 /product="hypothetical protein" 227 /protein_id="NP_963295.1" 228 /db_xref="GI:41614797" 229 /db_xref="GeneID:2732620" 230 /translation="MRLLLELKALNSIDKKQLSNYLIQGFIYNILKNTEYSWLHNWKK 231 EKYFNFTLIPKKDIIENKRYYLIISSPDKRFIEVLHNKIKDLDIITIGLAQFQLRKTK 232 KFDPKLRFPWVTITPIVLREGKIVILKGDKYYKVFVKRLEELKKYNLIKKKEPILEEP 233 IEISLNQIKDGWKIIDVKDRYYDFRNKSFSAFSNWLRDLKEQSLRKYNNFCGKNFYFE 234 EAIFEGFTFYKTVSIRIRINRGEAVYIGTLWKELNVYRKLDKEEREFYKFLYDCGLGS 235 LNSMGFGFVNTKKNSAR" 236 237 Then should give input key="CDS" and the rest of the data as a list of strings 238 lines=["complement(join(490883..490885,1..879))", ..., "LNSMGFGFVNTKKNSAR"] 239 where the leading spaces and trailing newlines have been removed. 240 241 Returns tuple containing: (key as string, location string, qualifiers as list) 242 as follows for this example: 243 244 key = "CDS", string 245 location = "complement(join(490883..490885,1..879))", string 246 qualifiers = list of string tuples: 247 248 [('locus_tag', '"NEQ001"'), 249 ('note', '"conserved hypothetical [Methanococcus jannaschii];\nCOG1583:..."'), 250 ('codon_start', '1'), 251 ('transl_table', '11'), 252 ('product', '"hypothetical protein"'), 253 ('protein_id', '"NP_963295.1"'), 254 ('db_xref', '"GI:41614797"'), 255 ('db_xref', '"GeneID:2732620"'), 256 ('translation', '"MRLLLELKALNSIDKKQLSNYLIQGFIYNILKNTEYSWLHNWKK\nEKYFNFT..."')] 257 258 In the above example, the "note" and "translation" were edited for compactness, 259 and they would contain multiple new line characters (displayed above as \n) 260 261 If a qualifier is quoted (in this case, everything except codon_start and 262 transl_table) then the quotes are NOT removed. 263 264 Note that no whitespace is removed. 265 """ 266 # Skip any blank lines 267 iterator = (x for x in lines if x) 268 try: 269 line = next(iterator) 270 271 feature_location = line.strip() 272 while feature_location[-1:] == ",": 273 # Multiline location, still more to come! 274 line = next(iterator) 275 feature_location += line.strip() 276 if feature_location.count("(") > feature_location.count(")"): 277 # Including the prev line in warning would be more explicit, 278 # but this way get one-and-only-one warning shown by default: 279 warnings.warn("Non-standard feature line wrapping (didn't break on comma)?", 280 BiopythonParserWarning) 281 while feature_location[-1:] == "," or feature_location.count("(") > feature_location.count(")"): 282 line = next(iterator) 283 feature_location += line.strip() 284 285 qualifiers = [] 286 287 for line_number, line in enumerate(iterator): 288 # check for extra wrapping of the location closing parentheses 289 if line_number == 0 and line.startswith(")"): 290 feature_location += line.strip() 291 elif line[0] == "/": 292 # New qualifier 293 i = line.find("=") 294 key = line[1:i] # does not work if i==-1 295 value = line[i + 1:] # we ignore 'value' if i==-1 296 if i == -1: 297 # Qualifier with no key, e.g. /pseudo 298 key = line[1:] 299 qualifiers.append((key, None)) 300 elif not value: 301 # ApE can output /note= 302 qualifiers.append((key, "")) 303 elif value == '"': 304 # One single quote 305 if self.debug: 306 print("Single quote %s:%s" % (key, value)) 307 # DO NOT remove the quote... 308 qualifiers.append((key, value)) 309 elif value[0] == '"': 310 # Quoted... 311 value_list = [value] 312 while value_list[-1][-1] != '"': 313 value_list.append(next(iterator)) 314 value = '\n'.join(value_list) 315 # DO NOT remove the quotes... 316 qualifiers.append((key, value)) 317 else: 318 # Unquoted 319 # if debug : print("Unquoted line %s:%s" % (key,value)) 320 qualifiers.append((key, value)) 321 else: 322 # Unquoted continuation 323 assert len(qualifiers) > 0 324 assert key == qualifiers[-1][0] 325 # if debug : print("Unquoted Cont %s:%s" % (key, line)) 326 if qualifiers[-1][1] is None: 327 raise StopIteration 328 qualifiers[-1] = (key, qualifiers[-1][1] + "\n" + line) 329 return feature_key, feature_location, qualifiers 330 except StopIteration: 331 # Bummer 332 raise ValueError("Problem with '%s' feature:\n%s" 333 % (feature_key, "\n".join(lines)))
334 357
358 - def _feed_first_line(self, consumer, line):
359 """Handle the LOCUS/ID line, passing data to the comsumer 360 361 This should be implemented by the EMBL / GenBank specific subclass 362 363 Used by the parse_records() and parse() methods. 364 """ 365 pass
366
367 - def _feed_header_lines(self, consumer, lines):
368 """Handle the header lines (list of strings), passing data to the comsumer 369 370 This should be implemented by the EMBL / GenBank specific subclass 371 372 Used by the parse_records() and parse() methods. 373 """ 374 pass
375 376 @staticmethod
377 - def _feed_feature_table(consumer, feature_tuples):
378 """Handle the feature table (list of tuples), passing data to the comsumer 379 380 Used by the parse_records() and parse() methods. 381 """ 382 consumer.start_feature_table() 383 for feature_key, location_string, qualifiers in feature_tuples: 384 consumer.feature_key(feature_key) 385 consumer.location(location_string) 386 for q_key, q_value in qualifiers: 387 if q_value is None: 388 consumer.feature_qualifier(q_key, q_value) 389 else: 390 consumer.feature_qualifier(q_key, q_value.replace("\n", " "))
391
392 - def _feed_misc_lines(self, consumer, lines):
393 """Handle any lines between features and sequence (list of strings), passing data to the consumer 394 395 This should be implemented by the EMBL / GenBank specific subclass 396 397 Used by the parse_records() and parse() methods. 398 """ 399 pass
400
401 - def feed(self, handle, consumer, do_features=True):
402 """Feed a set of data into the consumer. 403 404 This method is intended for use with the "old" code in Bio.GenBank 405 406 Arguments: 407 408 - handle - A handle with the information to parse. 409 - consumer - The consumer that should be informed of events. 410 - do_features - Boolean, should the features be parsed? 411 Skipping the features can be much faster. 412 413 Return values: 414 415 - true - Passed a record 416 - false - Did not find a record 417 """ 418 # Should work with both EMBL and GenBank files provided the 419 # equivalent Bio.GenBank._FeatureConsumer methods are called... 420 self.set_handle(handle) 421 if not self.find_start(): 422 # Could not find (another) record 423 consumer.data = None 424 return False 425 426 # We use the above class methods to parse the file into a simplified format. 427 # The first line, header lines and any misc lines after the features will be 428 # dealt with by GenBank / EMBL specific derived classes. 429 430 # First line and header: 431 self._feed_first_line(consumer, self.line) 432 self._feed_header_lines(consumer, self.parse_header()) 433 434 # Features (common to both EMBL and GenBank): 435 if do_features: 436 self._feed_feature_table(consumer, self.parse_features(skip=False)) 437 else: 438 self.parse_features(skip=True) # ignore the data 439 440 # Footer and sequence 441 misc_lines, sequence_string = self.parse_footer() 442 self._feed_misc_lines(consumer, misc_lines) 443 444 consumer.sequence(sequence_string) 445 # Calls to consumer.base_number() do nothing anyway 446 consumer.record_end("//") 447 448 assert self.line == "//" 449 450 # And we are done 451 return True
452
453 - def parse(self, handle, do_features=True):
454 """Returns a SeqRecord (with SeqFeatures if do_features=True) 455 456 See also the method parse_records() for use on multi-record files. 457 """ 458 from Bio.GenBank import _FeatureConsumer 459 from Bio.GenBank.utils import FeatureValueCleaner 460 461 consumer = _FeatureConsumer(use_fuzziness=1, 462 feature_cleaner=FeatureValueCleaner()) 463 464 if self.feed(handle, consumer, do_features): 465 return consumer.data 466 else: 467 return None
468
469 - def parse_records(self, handle, do_features=True):
470 """Returns a SeqRecord object iterator 471 472 Each record (from the ID/LOCUS line to the // line) becomes a SeqRecord 473 474 The SeqRecord objects include SeqFeatures if do_features=True 475 476 This method is intended for use in Bio.SeqIO 477 """ 478 # This is a generator function 479 while True: 480 record = self.parse(handle, do_features) 481 if record is None: 482 break 483 if record.id is None: 484 raise ValueError("Failed to parse the record's ID. Invalid ID line?") 485 if record.name == "<unknown name>": 486 raise ValueError("Failed to parse the record's name. Invalid ID line?") 487 if record.description == "<unknown description>": 488 raise ValueError("Failed to parse the record's description") 489 yield record
490
491 - def parse_cds_features(self, handle, 492 alphabet=generic_protein, 493 tags2id=('protein_id', 'locus_tag', 'product')):
494 """Returns SeqRecord object iterator 495 496 Each CDS feature becomes a SeqRecord. 497 498 - alphabet - Used for any sequence found in a translation field. 499 - tags2id - Tupple of three strings, the feature keys to use 500 for the record id, name and description, 501 502 This method is intended for use in Bio.SeqIO 503 """ 504 self.set_handle(handle) 505 while self.find_start(): 506 # Got an EMBL or GenBank record... 507 self.parse_header() # ignore header lines! 508 feature_tuples = self.parse_features() 509 # self.parse_footer() # ignore footer lines! 510 while True: 511 line = self.handle.readline() 512 if not line: 513 break 514 if line[:2] == "//": 515 break 516 self.line = line.rstrip() 517 518 # Now go though those features... 519 for key, location_string, qualifiers in feature_tuples: 520 if key == "CDS": 521 # Create SeqRecord 522 # ================ 523 # SeqRecord objects cannot be created with annotations, they 524 # must be added afterwards. So create an empty record and 525 # then populate it: 526 record = SeqRecord(seq=None) 527 annotations = record.annotations 528 529 # Should we add a location object to the annotations? 530 # I *think* that only makes sense for SeqFeatures with their 531 # sub features... 532 annotations['raw_location'] = location_string.replace(' ', '') 533 534 for (qualifier_name, qualifier_data) in qualifiers: 535 if qualifier_data is not None \ 536 and qualifier_data[0] == '"' and qualifier_data[-1] == '"': 537 # Remove quotes 538 qualifier_data = qualifier_data[1:-1] 539 # Append the data to the annotation qualifier... 540 if qualifier_name == "translation": 541 assert record.seq is None, "Multiple translations!" 542 record.seq = Seq(qualifier_data.replace("\n", ""), alphabet) 543 elif qualifier_name == "db_xref": 544 # its a list, possibly empty. Its safe to extend 545 record.dbxrefs.append(qualifier_data) 546 else: 547 if qualifier_data is not None: 548 qualifier_data = qualifier_data.replace("\n", " ").replace(" ", " ") 549 try: 550 annotations[qualifier_name] += " " + qualifier_data 551 except KeyError: 552 # Not an addition to existing data, its the first bit 553 annotations[qualifier_name] = qualifier_data 554 555 # Fill in the ID, Name, Description 556 # ================================= 557 try: 558 record.id = annotations[tags2id[0]] 559 except KeyError: 560 pass 561 try: 562 record.name = annotations[tags2id[1]] 563 except KeyError: 564 pass 565 try: 566 record.description = annotations[tags2id[2]] 567 except KeyError: 568 pass 569 570 yield record
571
572 573 -class EmblScanner(InsdcScanner):
574 """For extracting chunks of information in EMBL files""" 575 576 RECORD_START = "ID " 577 HEADER_WIDTH = 5 578 FEATURE_START_MARKERS = ["FH Key Location/Qualifiers", "FH"] 579 FEATURE_END_MARKERS = ["XX"] # XX can also mark the end of many things! 580 FEATURE_QUALIFIER_INDENT = 21 581 FEATURE_QUALIFIER_SPACER = "FT" + " " * (FEATURE_QUALIFIER_INDENT - 2) 582 SEQUENCE_HEADERS = ["SQ", "CO"] # Remove trailing spaces 583 584 EMBL_INDENT = HEADER_WIDTH 585 EMBL_SPACER = " " * EMBL_INDENT 586 623
624 - def _feed_first_line(self, consumer, line):
625 assert line[:self.HEADER_WIDTH].rstrip() == "ID" 626 if line[self.HEADER_WIDTH:].count(";") == 6: 627 # Looks like the semi colon separated style introduced in 2006 628 self._feed_first_line_new(consumer, line) 629 elif line[self.HEADER_WIDTH:].count(";") == 3: 630 if line.rstrip().endswith(" SQ"): 631 # EMBL-bank patent data 632 self._feed_first_line_patents(consumer, line) 633 else: 634 # Looks like the pre 2006 style 635 self._feed_first_line_old(consumer, line) 636 elif line[self.HEADER_WIDTH:].count(";") == 2: 637 # Looks like KIKO patent data 638 self._feed_first_line_patents_kipo(consumer, line) 639 else: 640 raise ValueError('Did not recognise the ID line layout:\n' + line)
641
642 - def _feed_first_line_patents(self, consumer, line):
643 # Old style EMBL patent records where ID line ended SQ 644 # Not 100% sure that PRT here is really molecule type and 645 # not the data file division... 646 # 647 # Either Non-Redundant Level 1 database records, 648 # ID <accession>; <molecule type>; <non-redundant level 1>; <cluster size L1> 649 # e.g. ID NRP_AX000635; PRT; NR1; 15 SQ 650 # 651 # Or, Non-Redundant Level 2 database records: 652 # ID <L2-accession>; <molecule type>; <non-redundant level 2>; <cluster size L2> 653 # e.g. ID NRP0000016E; PRT; NR2; 5 SQ 654 # e.g. ID NRP_AX000635; PRT; NR1; 15 SQ 655 fields = [data.strip() for data in line[self.HEADER_WIDTH:].strip()[:-3].split(";")] 656 assert len(fields) == 4 657 consumer.locus(fields[0]) 658 consumer.residue_type(fields[1]) # semi-redundant 659 consumer.data_file_division(fields[2])
660 # TODO - Record cluster size? 661
662 - def _feed_first_line_patents_kipo(self, consumer, line):
663 # EMBL format patent sequence from KIPO, e.g. 664 # ftp://ftp.ebi.ac.uk/pub/databases/patentdata/kipo_prt.dat.gz 665 # 666 # e.g. ID DI500001 STANDARD; PRT; 111 AA. 667 # 668 # This follows the style of _feed_first_line_old 669 assert line[:self.HEADER_WIDTH].rstrip() == "ID" 670 fields = [line[self.HEADER_WIDTH:].split(None, 1)[0]] 671 fields.extend(line[self.HEADER_WIDTH:].split(None, 1)[1].split(";")) 672 fields = [entry.strip() for entry in fields] 673 """ 674 The tokens represent: 675 676 0. Primary accession number 677 (space sep) 678 1. ??? (e.g. standard) 679 (semi-colon) 680 2. Molecule type (protein)? Division? Always 'PRT' 681 3. Sequence length (e.g. '111 AA.') 682 """ 683 consumer.locus(fields[0]) # Should we also call the accession consumer? 684 # consumer.molecule_type(fields[2]) 685 self._feed_seq_length(consumer, fields[3])
686
687 - def _feed_first_line_old(self, consumer, line):
688 # Expects an ID line in the style before 2006, e.g. 689 # ID SC10H5 standard; DNA; PRO; 4870 BP. 690 # ID BSUB9999 standard; circular DNA; PRO; 4214630 BP. 691 assert line[:self.HEADER_WIDTH].rstrip() == "ID" 692 fields = [line[self.HEADER_WIDTH:].split(None, 1)[0]] 693 fields.extend(line[self.HEADER_WIDTH:].split(None, 1)[1].split(";")) 694 fields = [entry.strip() for entry in fields] 695 """ 696 The tokens represent: 697 698 0. Primary accession number 699 (space sep) 700 1. ??? (e.g. standard) 701 (semi-colon) 702 2. Topology and/or Molecule type (e.g. 'circular DNA' or 'DNA') 703 3. Taxonomic division (e.g. 'PRO') 704 4. Sequence length (e.g. '4639675 BP.') 705 """ 706 consumer.locus(fields[0]) # Should we also call the accession consumer? 707 consumer.residue_type(fields[2]) 708 if "circular" in fields[2]: 709 consumer.topology("circular") 710 consumer.molecule_type(fields[2].replace("circular", "").strip()) 711 elif "linear" in fields[2]: 712 consumer.topology("linear") 713 consumer.molecule_type(fields[2].replace("linear", "").strip()) 714 else: 715 consumer.molecule_type(fields[2].strip()) 716 consumer.data_file_division(fields[3]) 717 self._feed_seq_length(consumer, fields[4])
718
719 - def _feed_first_line_new(self, consumer, line):
720 # Expects an ID line in the style introduced in 2006, e.g. 721 # ID X56734; SV 1; linear; mRNA; STD; PLN; 1859 BP. 722 # ID CD789012; SV 4; linear; genomic DNA; HTG; MAM; 500 BP. 723 assert line[:self.HEADER_WIDTH].rstrip() == "ID" 724 fields = [data.strip() for data in line[self.HEADER_WIDTH:].strip().split(";")] 725 assert len(fields) == 7 726 """ 727 The tokens represent: 728 729 0. Primary accession number 730 1. Sequence version number 731 2. Topology: 'circular' or 'linear' 732 3. Molecule type (e.g. 'genomic DNA') 733 4. Data class (e.g. 'STD') 734 5. Taxonomic division (e.g. 'PRO') 735 6. Sequence length (e.g. '4639675 BP.') 736 """ 737 738 consumer.locus(fields[0]) 739 740 # Call the accession consumer now, to make sure we record 741 # something as the record.id, in case there is no AC line 742 consumer.accession(fields[0]) 743 744 # TODO - How to deal with the version field? At the moment the consumer 745 # will try and use this for the ID which isn't ideal for EMBL files. 746 version_parts = fields[1].split() 747 if len(version_parts) == 2 \ 748 and version_parts[0] == "SV" \ 749 and version_parts[1].isdigit(): 750 consumer.version_suffix(version_parts[1]) 751 752 # Based on how the old GenBank parser worked, merge these two: 753 consumer.residue_type(" ".join(fields[2:4])) # Semi-obsolete 754 755 consumer.topology(fields[2]) 756 consumer.molecule_type(fields[3]) 757 758 # consumer.xxx(fields[4]) # TODO - What should we do with the data class? 759 760 consumer.data_file_division(fields[5]) 761 762 self._feed_seq_length(consumer, fields[6])
763 764 @staticmethod
765 - def _feed_seq_length(consumer, text):
766 length_parts = text.split() 767 assert len(length_parts) == 2, "Invalid sequence length string %r" % text 768 assert length_parts[1].upper() in ["BP", "BP.", "AA", "AA."] 769 consumer.size(length_parts[0])
770
771 - def _feed_header_lines(self, consumer, lines):
772 consumer_dict = { 773 'AC': 'accession', 774 'SV': 'version', # SV line removed in June 2006, now part of ID line 775 'DE': 'definition', 776 # 'RN' : 'reference_num', 777 # 'RC' : reference comment... TODO 778 # 'RP' : 'reference_bases', 779 # 'RX' : reference cross reference... DOI or Pubmed 780 'RG': 'consrtm', # optional consortium 781 # 'RA' : 'authors', 782 # 'RT' : 'title', 783 'RL': 'journal', 784 'OS': 'organism', 785 'OC': 'taxonomy', 786 # 'DR' : data reference 787 'CC': 'comment', 788 # 'XX' : splitter 789 } 790 # We have to handle the following specially: 791 # RX (depending on reference type...) 792 for line in lines: 793 line_type = line[:self.EMBL_INDENT].strip() 794 data = line[self.EMBL_INDENT:].strip() 795 if line_type == 'XX': 796 pass 797 elif line_type == 'RN': 798 # Reformat reference numbers for the GenBank based consumer 799 # e.g. '[1]' becomes '1' 800 if data[0] == "[" and data[-1] == "]": 801 data = data[1:-1] 802 consumer.reference_num(data) 803 elif line_type == 'RP': 804 if data.strip() == "[-]": 805 # Patent EMBL files from KIPO just use: RN [-] 806 pass 807 else: 808 # Reformat reference numbers for the GenBank based consumer 809 # e.g. '1-4639675' becomes '(bases 1 to 4639675)' 810 # and '160-550, 904-1055' becomes '(bases 160 to 550; 904 to 1055)' 811 # Note could be multi-line, and end with a comma 812 parts = [bases.replace("-", " to ").strip() for bases in data.split(",") if bases.strip()] 813 consumer.reference_bases("(bases %s)" % "; ".join(parts)) 814 elif line_type == 'RT': 815 # Remove the enclosing quotes and trailing semi colon. 816 # Note the title can be split over multiple lines. 817 if data.startswith('"'): 818 data = data[1:] 819 if data.endswith('";'): 820 data = data[:-2] 821 consumer.title(data) 822 elif line_type == 'RX': 823 # EMBL support three reference types at the moment: 824 # - PUBMED PUBMED bibliographic database (NLM) 825 # - DOI Digital Object Identifier (International DOI Foundation) 826 # - AGRICOLA US National Agriculture Library (NAL) of the US Department 827 # of Agriculture (USDA) 828 # 829 # Format: 830 # RX resource_identifier; identifier. 831 # 832 # e.g. 833 # RX DOI; 10.1016/0024-3205(83)90010-3. 834 # RX PUBMED; 264242. 835 # 836 # Currently our reference object only supports PUBMED and MEDLINE 837 # (as these were in GenBank files?). 838 key, value = data.split(";", 1) 839 if value.endswith("."): 840 value = value[:-1] 841 value = value.strip() 842 if key == "PUBMED": 843 consumer.pubmed_id(value) 844 # TODO - Handle other reference types (here and in BioSQL bindings) 845 elif line_type == 'CC': 846 # Have to pass a list of strings for this one (not just a string) 847 consumer.comment([data]) 848 elif line_type == 'DR': 849 # Database Cross-reference, format: 850 # DR database_identifier; primary_identifier; secondary_identifier. 851 # 852 # e.g. 853 # DR MGI; 98599; Tcrb-V4. 854 # 855 # TODO - How should we store any secondary identifier? 856 parts = data.rstrip(".").split(";") 857 # Turn it into "database_identifier:primary_identifier" to 858 # mimic the GenBank parser. e.g. "MGI:98599" 859 consumer.dblink("%s:%s" % (parts[0].strip(), 860 parts[1].strip())) 861 elif line_type == 'RA': 862 # Remove trailing ; at end of authors list 863 consumer.authors(data.rstrip(";")) 864 elif line_type == 'PR': 865 # In the EMBL patent files, this is a PR (PRiority) line which 866 # provides the earliest active priority within the family. 867 # The priority number comes first, followed by the priority date. 868 # 869 # e.g. 870 # PR JP19990377484 16-DEC-1999 871 # 872 # However, in most EMBL files this is a PR (PRoject) line which 873 # gives the BioProject reference number. 874 # 875 # e.g. 876 # PR Project:PRJNA60715; 877 # 878 # In GenBank files this corresponds to the old PROJECT line 879 # which was later replaced with the DBLINK line. 880 if data.startswith("Project:"): 881 # Remove trailing ; at end of the project reference 882 consumer.project(data.rstrip(";")) 883 elif line_type == 'KW': 884 consumer.keywords(data.rstrip(";")) 885 elif line_type in consumer_dict: 886 # Its a semi-automatic entry! 887 getattr(consumer, consumer_dict[line_type])(data) 888 else: 889 if self.debug: 890 print("Ignoring EMBL header line:\n%s" % line)
891
892 - def _feed_misc_lines(self, consumer, lines):
893 # TODO - Should we do something with the information on the SQ line(s)? 894 lines.append("") 895 line_iter = iter(lines) 896 try: 897 for line in line_iter: 898 if line.startswith("CO "): 899 line = line[5:].strip() 900 contig_location = line 901 while True: 902 line = next(line_iter) 903 if not line: 904 break 905 elif line.startswith("CO "): 906 # Don't need to preseve the whitespace here. 907 contig_location += line[5:].strip() 908 else: 909 raise ValueError('Expected CO (contig) continuation line, got:\n' + line) 910 consumer.contig_location(contig_location) 911 if line.startswith("SQ Sequence "): 912 # e.g. 913 # SQ Sequence 219 BP; 82 A; 48 C; 33 G; 45 T; 11 other; 914 # 915 # Or, EMBL-bank patent, e.g. 916 # SQ Sequence 465 AA; 3963407aa91d3a0d622fec679a4524e0; MD5; 917 self._feed_seq_length(consumer, line[14:].rstrip().rstrip(";").split(";", 1)[0]) 918 # TODO - Record the checksum etc? 919 return 920 except StopIteration: 921 raise ValueError("Problem in misc lines before sequence")
922
923 924 -class _ImgtScanner(EmblScanner):
925 """For extracting chunks of information in IMGT (EMBL like) files (PRIVATE). 926 927 IMGT files are like EMBL files but in order to allow longer feature types 928 the features should be indented by 25 characters not 21 characters. In 929 practice the IMGT flat files tend to use either 21 or 25 characters, so we 930 must cope with both. 931 932 This is private to encourage use of Bio.SeqIO rather than Bio.GenBank. 933 """ 934 935 FEATURE_START_MARKERS = ["FH Key Location/Qualifiers", 936 "FH Key Location/Qualifiers (from EMBL)", 937 "FH Key Location/Qualifiers", 938 "FH"] 939
940 - def _feed_first_line(self, consumer, line):
941 assert line[:self.HEADER_WIDTH].rstrip() == "ID" 942 if line[self.HEADER_WIDTH:].count(";") != 5: 943 # Assume its an older EMBL-like line, 944 return EmblScanner._feed_first_line(self, consumer, line) 945 # Otherwise assume its the new (circa 2016) IMGT style 946 # as used in the IPD-IMGT/HLA Database 947 # 948 # https://github.com/ANHIG/IMGTHLA/ 949 # 950 # The key changes post 3.16 are the addition of an SV value 951 # to the ID line, these additions should make the format more 952 # similar to the ENA style. 953 # 954 # ID HLA00001 standard; DNA; HUM; 3503 BP. 955 # 956 # becomes 957 # 958 # ID HLA00001; SV 1; standard; DNA; HUM; 3503 BP. 959 fields = [data.strip() for data in line[self.HEADER_WIDTH:].strip().split(";")] 960 assert len(fields) == 6 961 """ 962 The tokens represent: 963 964 0. Primary accession number (eg 'HLA00001') 965 1. Sequence version number (eg 'SV 1') 966 2. ??? eg 'standard' 967 3. Molecule type (e.g. 'DNA') 968 4. Taxonomic division (e.g. 'HUM') 969 5. Sequence length (e.g. '3503 BP.') 970 """ 971 consumer.locus(fields[0]) 972 973 # See TODO on the EMBL _feed_first_line_new about version field 974 version_parts = fields[1].split() 975 if len(version_parts) == 2 \ 976 and version_parts[0] == "SV" \ 977 and version_parts[1].isdigit(): 978 consumer.version_suffix(version_parts[1]) 979 980 consumer.residue_type(fields[3]) 981 if "circular" in fields[3]: 982 consumer.topology("circular") 983 consumer.molecule_type(fields[3].replace("circular", "").strip()) 984 elif "linear" in fields[3]: 985 consumer.topology("linear") 986 consumer.molecule_type(fields[3].replace("linear", "").strip()) 987 else: 988 consumer.molecule_type(fields[3].strip()) 989 consumer.data_file_division(fields[4]) 990 self._feed_seq_length(consumer, fields[5])
991
992 - def parse_features(self, skip=False):
993 """Return list of tuples for the features (if present) 994 995 Each feature is returned as a tuple (key, location, qualifiers) 996 where key and location are strings (e.g. "CDS" and 997 "complement(join(490883..490885,1..879))") while qualifiers 998 is a list of two string tuples (feature qualifier keys and values). 999 1000 Assumes you have already read to the start of the features table. 1001 """ 1002 if self.line.rstrip() not in self.FEATURE_START_MARKERS: 1003 if self.debug: 1004 print("Didn't find any feature table") 1005 return [] 1006 1007 while self.line.rstrip() in self.FEATURE_START_MARKERS: 1008 self.line = self.handle.readline() 1009 1010 bad_position_re = re.compile(r'([0-9]+)>') 1011 1012 features = [] 1013 line = self.line 1014 while True: 1015 if not line: 1016 raise ValueError("Premature end of line during features table") 1017 if line[:self.HEADER_WIDTH].rstrip() in self.SEQUENCE_HEADERS: 1018 if self.debug: 1019 print("Found start of sequence") 1020 break 1021 line = line.rstrip() 1022 if line == "//": 1023 raise ValueError("Premature end of features table, marker '//' found") 1024 if line in self.FEATURE_END_MARKERS: 1025 if self.debug: 1026 print("Found end of features") 1027 line = self.handle.readline() 1028 break 1029 if line[2:self.FEATURE_QUALIFIER_INDENT].strip() == "": 1030 # This is an empty feature line between qualifiers. Empty 1031 # feature lines within qualifiers are handled below (ignored). 1032 line = self.handle.readline() 1033 continue 1034 1035 if skip: 1036 line = self.handle.readline() 1037 while line[:self.FEATURE_QUALIFIER_INDENT] == self.FEATURE_QUALIFIER_SPACER: 1038 line = self.handle.readline() 1039 else: 1040 assert line[:2] == "FT" 1041 try: 1042 feature_key, location_start = line[2:].strip().split() 1043 except ValueError: 1044 # e.g. "FT TRANSMEMBRANE-REGION2163..2240\n" 1045 # Assume indent of 25 as per IMGT spec, with the location 1046 # start in column 26 (one-based). 1047 feature_key = line[2:25].strip() 1048 location_start = line[25:].strip() 1049 feature_lines = [location_start] 1050 line = self.handle.readline() 1051 while line[:self.FEATURE_QUALIFIER_INDENT] == self.FEATURE_QUALIFIER_SPACER \ 1052 or line.rstrip() == "": # cope with blank lines in the midst of a feature 1053 # Use strip to remove any harmless trailing white space AND and leading 1054 # white space (copes with 21 or 26 indents and orther variants) 1055 assert line[:2] == "FT" 1056 feature_lines.append(line[self.FEATURE_QUALIFIER_INDENT:].strip()) 1057 line = self.handle.readline() 1058 feature_key, location, qualifiers = \ 1059 self.parse_feature(feature_key, feature_lines) 1060 # Try to handle known problems with IMGT locations here: 1061 if ">" in location: 1062 # Nasty hack for common IMGT bug, should be >123 not 123> 1063 # in a location string. At least here the meaning is clear, 1064 # and since it is so common I don't want to issue a warning 1065 # warnings.warn("Feature location %s is invalid, " 1066 # "moving greater than sign before position" 1067 # % location, BiopythonParserWarning) 1068 location = bad_position_re.sub(r'>\1', location) 1069 features.append((feature_key, location, qualifiers)) 1070 self.line = line 1071 return features
1072
1073 1074 -class GenBankScanner(InsdcScanner):
1075 """For extracting chunks of information in GenBank files""" 1076 RECORD_START = "LOCUS " 1077 HEADER_WIDTH = 12 1078 FEATURE_START_MARKERS = ["FEATURES Location/Qualifiers", "FEATURES"] 1079 FEATURE_END_MARKERS = [] 1080 FEATURE_QUALIFIER_INDENT = 21 1081 FEATURE_QUALIFIER_SPACER = " " * FEATURE_QUALIFIER_INDENT 1082 SEQUENCE_HEADERS = ["CONTIG", "ORIGIN", "BASE COUNT", "WGS"] # trailing spaces removed 1083 1084 GENBANK_INDENT = HEADER_WIDTH 1085 GENBANK_SPACER = " " * GENBANK_INDENT 1086 1087 STRUCTURED_COMMENT_START = "-START##" 1088 STRUCTURED_COMMENT_END = "-END##" 1089 STRUCTURED_COMMENT_DELIM = " :: " 1090 1143
1144 - def _feed_first_line(self, consumer, line):
1145 """Scan over and parse GenBank LOCUS line (PRIVATE). 1146 1147 This must cope with several variants, primarily the old and new column 1148 based standards from GenBank. Additionally EnsEMBL produces GenBank 1149 files where the LOCUS line is space separated rather that following 1150 the column based layout. 1151 1152 We also try to cope with GenBank like files with partial LOCUS lines. 1153 """ 1154 ##################################### 1155 # LOCUS line # 1156 ##################################### 1157 assert line[0:self.GENBANK_INDENT] == 'LOCUS ', \ 1158 'LOCUS line does not start correctly:\n' + line 1159 1160 # Have to break up the locus line, and handle the different bits of it. 1161 # There are at least two different versions of the locus line... 1162 if line[29:33] in [' bp ', ' aa ', ' rc '] and line[55:62] == ' ': 1163 # Old... note we insist on the 55:62 being empty to avoid trying 1164 # to parse space separated LOCUS lines from Ensembl etc, see below. 1165 # 1166 # Positions Contents 1167 # --------- -------- 1168 # 00:06 LOCUS 1169 # 06:12 spaces 1170 # 12:?? Locus name 1171 # ??:?? space 1172 # ??:29 Length of sequence, right-justified 1173 # 29:33 space, bp, space 1174 # 33:41 strand type / molecule type, e.g. DNA 1175 # 41:42 space 1176 # 42:51 Blank (implies linear), linear or circular 1177 # 51:52 space 1178 # 52:55 The division code (e.g. BCT, VRL, INV) 1179 # 55:62 space 1180 # 62:73 Date, in the form dd-MMM-yyyy (e.g., 15-MAR-1991) 1181 # 1182 # assert line[29:33] in [' bp ', ' aa ',' rc '] , \ 1183 # 'LOCUS line does not contain size units at expected position:\n' + line 1184 assert line[41:42] == ' ', \ 1185 'LOCUS line does not contain space at position 42:\n' + line 1186 assert line[42:51].strip() in ['', 'linear', 'circular'], \ 1187 'LOCUS line does not contain valid entry (linear, circular, ...):\n' + line 1188 assert line[51:52] == ' ', \ 1189 'LOCUS line does not contain space at position 52:\n' + line 1190 # assert line[55:62] == ' ', \ 1191 # 'LOCUS line does not contain spaces from position 56 to 62:\n' + line 1192 if line[62:73].strip(): 1193 assert line[64:65] == '-', \ 1194 'LOCUS line does not contain - at position 65 in date:\n' + line 1195 assert line[68:69] == '-', \ 1196 'LOCUS line does not contain - at position 69 in date:\n' + line 1197 1198 name_and_length_str = line[self.GENBANK_INDENT:29] 1199 while ' ' in name_and_length_str: 1200 name_and_length_str = name_and_length_str.replace(' ', ' ') 1201 name_and_length = name_and_length_str.split(' ') 1202 assert len(name_and_length) <= 2, \ 1203 'Cannot parse the name and length in the LOCUS line:\n' + line 1204 assert len(name_and_length) != 1, \ 1205 'Name and length collide in the LOCUS line:\n' + line 1206 # Should be possible to split them based on position, if 1207 # a clear definition of the standard exists THAT AGREES with 1208 # existing files. 1209 name, length = name_and_length 1210 if len(name) > 16: 1211 # As long as the sequence is short, can steal its leading spaces 1212 # to extend the name over the current 16 character limit. 1213 # However, that deserves a warning as it is out of spec. 1214 warnings.warn("GenBank LOCUS line identifier over 16 characters", 1215 BiopythonParserWarning) 1216 consumer.locus(name) 1217 consumer.size(length) 1218 # consumer.residue_type(line[33:41].strip()) 1219 1220 if line[33:51].strip() == "" and line[29:33] == ' aa ': 1221 # Amino acids -> protein (even if there is no residue type given) 1222 # We want to use a protein alphabet in this case, rather than a 1223 # generic one. Not sure if this is the best way to achieve this, 1224 # but it works because the scanner checks for this: 1225 consumer.residue_type("PROTEIN") 1226 else: 1227 consumer.residue_type(line[33:51].strip()) 1228 1229 consumer.molecule_type(line[33:41].strip()) 1230 consumer.topology(line[42:51].strip()) 1231 consumer.data_file_division(line[52:55]) 1232 if line[62:73].strip(): 1233 consumer.date(line[62:73]) 1234 elif line[40:44] in [' bp ', ' aa ', ' rc '] \ 1235 and line[54:64].strip() in ['', 'linear', 'circular']: 1236 # New... linear/circular/big blank test should avoid EnsEMBL style 1237 # LOCUS line being treated like a proper column based LOCUS line. 1238 # 1239 # Positions Contents 1240 # --------- -------- 1241 # 00:06 LOCUS 1242 # 06:12 spaces 1243 # 12:?? Locus name 1244 # ??:?? space 1245 # ??:40 Length of sequence, right-justified 1246 # 40:44 space, bp, space 1247 # 44:47 Blank, ss-, ds-, ms- 1248 # 47:54 Blank, DNA, RNA, tRNA, mRNA, uRNA, snRNA, cDNA 1249 # 54:55 space 1250 # 55:63 Blank (implies linear), linear or circular 1251 # 63:64 space 1252 # 64:67 The division code (e.g. BCT, VRL, INV) 1253 # 67:68 space 1254 # 68:79 Date, in the form dd-MMM-yyyy (e.g., 15-MAR-1991) 1255 # 1256 assert line[40:44] in [' bp ', ' aa ', ' rc '], \ 1257 'LOCUS line does not contain size units at expected position:\n' + line 1258 assert line[44:47] in [' ', 'ss-', 'ds-', 'ms-'], \ 1259 'LOCUS line does not have valid strand type (Single stranded, ...):\n' + line 1260 assert line[47:54].strip() == "" \ 1261 or 'DNA' in line[47:54].strip().upper() \ 1262 or 'RNA' in line[47:54].strip().upper(), \ 1263 'LOCUS line does not contain valid sequence type (DNA, RNA, ...):\n' + line 1264 assert line[54:55] == ' ', \ 1265 'LOCUS line does not contain space at position 55:\n' + line 1266 assert line[55:63].strip() in ['', 'linear', 'circular'], \ 1267 'LOCUS line does not contain valid entry (linear, circular, ...):\n' + line 1268 assert line[63:64] == ' ', \ 1269 'LOCUS line does not contain space at position 64:\n' + line 1270 assert line[67:68] == ' ', \ 1271 'LOCUS line does not contain space at position 68:\n' + line 1272 if line[68:79].strip(): 1273 assert line[70:71] == '-', \ 1274 'LOCUS line does not contain - at position 71 in date:\n' + line 1275 assert line[74:75] == '-', \ 1276 'LOCUS line does not contain - at position 75 in date:\n' + line 1277 1278 name_and_length_str = line[self.GENBANK_INDENT:40] 1279 while ' ' in name_and_length_str: 1280 name_and_length_str = name_and_length_str.replace(' ', ' ') 1281 name_and_length = name_and_length_str.split(' ') 1282 assert len(name_and_length) <= 2, \ 1283 'Cannot parse the name and length in the LOCUS line:\n' + line 1284 assert len(name_and_length) != 1, \ 1285 'Name and length collide in the LOCUS line:\n' + line 1286 # Should be possible to split them based on position, if 1287 # a clear definition of the stand exists THAT AGREES with 1288 # existing files. 1289 consumer.locus(name_and_length[0]) 1290 consumer.size(name_and_length[1]) 1291 1292 if line[44:54].strip() == "" and line[40:44] == ' aa ': 1293 # Amino acids -> protein (even if there is no residue type given) 1294 # We want to use a protein alphabet in this case, rather than a 1295 # generic one. Not sure if this is the best way to achieve this, 1296 # but it works because the scanner checks for this: 1297 consumer.residue_type(("PROTEIN " + line[54:63]).strip()) 1298 else: 1299 consumer.residue_type(line[44:63].strip()) 1300 1301 consumer.molecule_type(line[44:54].strip()) 1302 consumer.topology(line[55:63].strip()) 1303 consumer.data_file_division(line[64:67]) 1304 if line[68:79].strip(): 1305 consumer.date(line[68:79]) 1306 elif line[self.GENBANK_INDENT:].strip().count(" ") == 0: 1307 # Truncated LOCUS line, as produced by some EMBOSS tools - see bug 1762 1308 # 1309 # e.g. 1310 # 1311 # "LOCUS U00096" 1312 # 1313 # rather than: 1314 # 1315 # "LOCUS U00096 4639675 bp DNA circular BCT" 1316 # 1317 # Positions Contents 1318 # --------- -------- 1319 # 00:06 LOCUS 1320 # 06:12 spaces 1321 # 12:?? Locus name 1322 if line[self.GENBANK_INDENT:].strip() != "": 1323 consumer.locus(line[self.GENBANK_INDENT:].strip()) 1324 else: 1325 # Must just have just "LOCUS ", is this even legitimate? 1326 # We should be able to continue parsing... we need real world testcases! 1327 warnings.warn("Minimal LOCUS line found - is this " 1328 "correct?\n:%r" % line, BiopythonParserWarning) 1329 elif len(line.split()) == 8 and line.split()[3] in ("aa", "bp") \ 1330 and line.split()[5] in ('linear', 'circular'): 1331 # Cope with invalidly spaced GenBank LOCUS lines like 1332 # LOCUS AB070938 6497 bp DNA linear BCT 11-OCT-2001 1333 splitline = line.split() 1334 consumer.locus(splitline[1]) 1335 consumer.size(splitline[2]) 1336 consumer.residue_type(splitline[4]) 1337 consumer.topology(splitline[5]) 1338 consumer.data_file_division(splitline[6]) 1339 consumer.date(splitline[7]) 1340 warnings.warn("Attempting to parse malformed locus line:\n%r\n" 1341 "Found locus %r size %r residue_type %r\n" 1342 "Some fields may be wrong." 1343 % (line, splitline[1], splitline[2], splitline[4]), 1344 BiopythonParserWarning) 1345 elif len(line.split()) == 7 and line.split()[3] in ["aa", "bp"]: 1346 # Cope with EnsEMBL genbank files which use space separation rather 1347 # than the expected column based layout. e.g. 1348 # LOCUS HG531_PATCH 1000000 bp DNA HTG 18-JUN-2011 1349 # LOCUS HG531_PATCH 759984 bp DNA HTG 18-JUN-2011 1350 # LOCUS HG506_HG1000_1_PATCH 814959 bp DNA HTG 18-JUN-2011 1351 # LOCUS HG506_HG1000_1_PATCH 1219964 bp DNA HTG 18-JUN-2011 1352 # Notice that the 'bp' can occur in the position expected by either 1353 # the old or the new fixed column standards (parsed above). 1354 splitline = line.split() 1355 consumer.locus(splitline[1]) 1356 consumer.size(splitline[2]) 1357 consumer.residue_type(splitline[4]) 1358 consumer.data_file_division(splitline[5]) 1359 consumer.date(splitline[6]) 1360 elif len(line.split()) >= 4 and line.split()[3] in ["aa", "bp"]: 1361 # Cope with EMBOSS seqret output where it seems the locus id can cause 1362 # the other fields to overflow. We just IGNORE the other fields! 1363 warnings.warn("Malformed LOCUS line found - is this " 1364 "correct?\n:%r" % line, BiopythonParserWarning) 1365 consumer.locus(line.split()[1]) 1366 consumer.size(line.split()[2]) 1367 elif len(line.split()) >= 4 and line.split()[-1] in ["aa", "bp"]: 1368 # Cope with pseudo-GenBank files like this: 1369 # "LOCUS RNA5 complete 1718 bp" 1370 # Treat everything between LOCUS and the size as the identifier. 1371 warnings.warn("Malformed LOCUS line found - is this " 1372 "correct?\n:%r" % line, BiopythonParserWarning) 1373 consumer.locus(line[5:].rsplit(None, 2)[0].strip()) 1374 consumer.size(line.split()[-2]) 1375 else: 1376 raise ValueError('Did not recognise the LOCUS line layout:\n' + line)
1377
1378 - def _feed_header_lines(self, consumer, lines):
1379 # Following dictionary maps GenBank lines to the associated 1380 # consumer methods - the special cases like LOCUS where one 1381 # genbank line triggers several consumer calls have to be 1382 # handled individually. 1383 consumer_dict = { 1384 'DEFINITION': 'definition', 1385 'ACCESSION': 'accession', 1386 'NID': 'nid', 1387 'PID': 'pid', 1388 'DBSOURCE': 'db_source', 1389 'KEYWORDS': 'keywords', 1390 'SEGMENT': 'segment', 1391 'SOURCE': 'source', 1392 'AUTHORS': 'authors', 1393 'CONSRTM': 'consrtm', 1394 'PROJECT': 'project', 1395 'TITLE': 'title', 1396 'JOURNAL': 'journal', 1397 'MEDLINE': 'medline_id', 1398 'PUBMED': 'pubmed_id', 1399 'REMARK': 'remark'} 1400 # We have to handle the following specially: 1401 # ORIGIN (locus, size, residue_type, data_file_division and date) 1402 # COMMENT (comment) 1403 # VERSION (version and gi) 1404 # DBLINK (database links like projects, newlines important) 1405 # REFERENCE (eference_num and reference_bases) 1406 # ORGANISM (organism and taxonomy) 1407 lines = [_f for _f in lines if _f] 1408 lines.append("") # helps avoid getting StopIteration all the time 1409 line_iter = iter(lines) 1410 try: 1411 line = next(line_iter) 1412 while True: 1413 if not line: 1414 break 1415 line_type = line[:self.GENBANK_INDENT].strip() 1416 data = line[self.GENBANK_INDENT:].strip() 1417 1418 if line_type == 'VERSION': 1419 # Need to call consumer.version(), and maybe also consumer.gi() as well. 1420 # e.g. 1421 # VERSION AC007323.5 GI:6587720 1422 while ' ' in data: 1423 data = data.replace(' ', ' ') 1424 if ' GI:' not in data: 1425 consumer.version(data) 1426 else: 1427 if self.debug: 1428 print("Version [" + data.split(' GI:')[0] + "], gi [" + data.split(' GI:')[1] + "]") 1429 consumer.version(data.split(' GI:')[0]) 1430 consumer.gi(data.split(' GI:')[1]) 1431 # Read in the next line! 1432 line = next(line_iter) 1433 elif line_type == 'DBLINK': 1434 # Need to call consumer.dblink() for each line, e.g. 1435 # DBLINK Project: 57779 1436 # BioProject: PRJNA57779 1437 consumer.dblink(data.strip()) 1438 # Read in the next line, and see if its more of the DBLINK section: 1439 while True: 1440 line = next(line_iter) 1441 if line[:self.GENBANK_INDENT] == self.GENBANK_SPACER: 1442 # Add this continuation to the data string 1443 consumer.dblink(line[self.GENBANK_INDENT:].strip()) 1444 else: 1445 # End of the DBLINK, leave this text in the variable "line" 1446 break 1447 elif line_type == 'REFERENCE': 1448 if self.debug > 1: 1449 print("Found reference [" + data + "]") 1450 # Need to call consumer.reference_num() and consumer.reference_bases() 1451 # e.g. 1452 # REFERENCE 1 (bases 1 to 86436) 1453 # 1454 # Note that this can be multiline, see Bug 1968, e.g. 1455 # 1456 # REFERENCE 42 (bases 1517 to 1696; 3932 to 4112; 17880 to 17975; 21142 to 1457 # 28259) 1458 # 1459 # For such cases we will call the consumer once only. 1460 data = data.strip() 1461 1462 # Read in the next line, and see if its more of the reference: 1463 while True: 1464 line = next(line_iter) 1465 if line[:self.GENBANK_INDENT] == self.GENBANK_SPACER: 1466 # Add this continuation to the data string 1467 data += " " + line[self.GENBANK_INDENT:] 1468 if self.debug > 1: 1469 print("Extended reference text [" + data + "]") 1470 else: 1471 # End of the reference, leave this text in the variable "line" 1472 break 1473 1474 # We now have all the reference line(s) stored in a string, data, 1475 # which we pass to the consumer 1476 while ' ' in data: 1477 data = data.replace(' ', ' ') 1478 if ' ' not in data: 1479 if self.debug > 2: 1480 print('Reference number \"' + data + '\"') 1481 consumer.reference_num(data) 1482 else: 1483 if self.debug > 2: 1484 print('Reference number \"' + data[:data.find(' ')] + '\", \"' + data[data.find(' ') + 1:] + '\"') 1485 consumer.reference_num(data[:data.find(' ')]) 1486 consumer.reference_bases(data[data.find(' ') + 1:]) 1487 elif line_type == 'ORGANISM': 1488 # Typically the first line is the organism, and subsequent lines 1489 # are the taxonomy lineage. However, given longer and longer 1490 # species names (as more and more strains and sub strains get 1491 # sequenced) the oragnism name can now get wrapped onto multiple 1492 # lines. The NCBI say we have to recognise the lineage line by 1493 # the presence of semi-colon delimited entries. In the long term, 1494 # they are considering adding a new keyword (e.g. LINEAGE). 1495 # See Bug 2591 for details. 1496 organism_data = data 1497 lineage_data = "" 1498 while True: 1499 line = next(line_iter) 1500 if line[0:self.GENBANK_INDENT] == self.GENBANK_SPACER: 1501 if lineage_data or ";" in line: 1502 lineage_data += " " + line[self.GENBANK_INDENT:] 1503 elif line[self.GENBANK_INDENT:].strip() == ".": 1504 # No lineage data, just . place holder 1505 pass 1506 else: 1507 organism_data += " " + line[self.GENBANK_INDENT:].strip() 1508 else: 1509 # End of organism and taxonomy 1510 break 1511 consumer.organism(organism_data) 1512 if lineage_data.strip() == "" and self.debug > 1: 1513 print("Taxonomy line(s) missing or blank") 1514 consumer.taxonomy(lineage_data.strip()) 1515 del organism_data, lineage_data 1516 elif line_type == 'COMMENT': 1517 # A COMMENT can either be plain text or tabular (Structured Comment), 1518 # or contain both. Multi-line comments are common. The code calls 1519 # consumer.comment() once with a list where each entry 1520 # is a line. If there's a structured comment consumer.structured_comment() 1521 # is called with a dict of dicts where the secondary key/value pairs are 1522 # the same as those in the structured comment table. The primary key is 1523 # the title or header of the table (e.g. Assembly-Data, FluData). See 1524 # http://www.ncbi.nlm.nih.gov/genbank/structuredcomment 1525 # for more information on Structured Comments. 1526 data = line[self.GENBANK_INDENT:] 1527 if self.debug > 1: 1528 print("Found comment") 1529 comment_list = [] 1530 structured_comment_dict = OrderedDict() 1531 regex = r"([^#]+){0}$".format(self.STRUCTURED_COMMENT_START) 1532 structured_comment_key = re.search(regex, data) 1533 if structured_comment_key is not None: 1534 structured_comment_key = structured_comment_key.group(1) 1535 if self.debug > 1: 1536 print("Found Structured Comment") 1537 else: 1538 comment_list.append(data) 1539 1540 while True: 1541 line = next(line_iter) 1542 data = line[self.GENBANK_INDENT:] 1543 if line[0:self.GENBANK_INDENT] == self.GENBANK_SPACER: 1544 if self.STRUCTURED_COMMENT_START in data: 1545 regex = r"([^#]+){0}$".format(self.STRUCTURED_COMMENT_START) 1546 structured_comment_key = re.search(regex, data) 1547 if structured_comment_key is not None: 1548 structured_comment_key = structured_comment_key.group(1) 1549 else: 1550 comment_list.append(data) 1551 elif structured_comment_key is not None and self.STRUCTURED_COMMENT_DELIM in data: 1552 match = re.search(r"(.+?)\s*{0}\s*(.+)".format(self.STRUCTURED_COMMENT_DELIM), data) 1553 structured_comment_dict.setdefault(structured_comment_key, OrderedDict()) 1554 structured_comment_dict[structured_comment_key][match.group(1)] = match.group(2) 1555 if self.debug > 2: 1556 print("Structured Comment continuation [" + data + "]") 1557 elif self.STRUCTURED_COMMENT_END not in data: 1558 comment_list.append(data) 1559 if self.debug > 2: 1560 print("Comment continuation [" + data + "]") 1561 else: 1562 # End of the comment 1563 break 1564 if comment_list: 1565 consumer.comment(comment_list) 1566 if structured_comment_dict: 1567 consumer.structured_comment(structured_comment_dict) 1568 del comment_list, structured_comment_key, structured_comment_dict 1569 elif line_type in consumer_dict: 1570 # It's a semi-automatic entry! 1571 # Now, this may be a multi line entry... 1572 while True: 1573 line = next(line_iter) 1574 if line[0:self.GENBANK_INDENT] == self.GENBANK_SPACER: 1575 data += ' ' + line[self.GENBANK_INDENT:] 1576 else: 1577 # We now have all the data for this entry: 1578 1579 # The DEFINITION field must ends with a period 1580 # # see ftp://ftp.ncbi.nih.gov/genbank/gbrel.txt [3.4.5] 1581 # and discussion https://github.com/biopython/biopython/pull/616 1582 # We consider this period belong to the syntax, not to the data 1583 # So remove it if it exist 1584 if line_type == 'DEFINITION' and data.endswith('.'): 1585 data = data[:-1] 1586 getattr(consumer, consumer_dict[line_type])(data) 1587 # End of continuation - return to top of loop! 1588 break 1589 else: 1590 if self.debug: 1591 print("Ignoring GenBank header line:\n" % line) 1592 # Read in next line 1593 line = next(line_iter) 1594 except StopIteration: 1595 raise ValueError("Problem in header")
1596
1597 - def _feed_misc_lines(self, consumer, lines):
1598 # Deals with a few misc lines between the features and the sequence 1599 lines.append("") 1600 line_iter = iter(lines) 1601 try: 1602 for line in line_iter: 1603 if line.startswith('BASE COUNT'): 1604 line = line[10:].strip() 1605 if line: 1606 if self.debug: 1607 print("base_count = " + line) 1608 consumer.base_count(line) 1609 if line.startswith('ORIGIN'): 1610 line = line[6:].strip() 1611 if line: 1612 if self.debug: 1613 print("origin_name = " + line) 1614 consumer.origin_name(line) 1615 if line.startswith('WGS '): 1616 line = line[3:].strip() 1617 consumer.wgs(line) 1618 if line.startswith('WGS_SCAFLD'): 1619 line = line[10:].strip() 1620 consumer.add_wgs_scafld(line) 1621 if line.startswith('CONTIG'): 1622 line = line[6:].strip() 1623 contig_location = line 1624 while True: 1625 line = next(line_iter) 1626 if not line: 1627 break 1628 elif line[:self.GENBANK_INDENT] == self.GENBANK_SPACER: 1629 # Don't need to preseve the whitespace here. 1630 contig_location += line[self.GENBANK_INDENT:].rstrip() 1631 elif line.startswith('ORIGIN'): 1632 # Strange, seen this in GenPept files via Entrez gbwithparts 1633 line = line[6:].strip() 1634 if line: 1635 consumer.origin_name(line) 1636 break 1637 else: 1638 raise ValueError('Expected CONTIG continuation line, got:\n' + line) 1639 consumer.contig_location(contig_location) 1640 return 1641 except StopIteration: 1642 raise ValueError("Problem in misc lines before sequence")
1643 1644 if __name__ == "__main__": 1645 from Bio._py3k import StringIO 1646 1647 gbk_example = \ 1648 """LOCUS SCU49845 5028 bp DNA PLN 21-JUN-1999 1649 DEFINITION Saccharomyces cerevisiae TCP1-beta gene, partial cds, and Axl2p 1650 (AXL2) and Rev7p (REV7) genes, complete cds. 1651 ACCESSION U49845 1652 VERSION U49845.1 GI:1293613 1653 KEYWORDS . 1654 SOURCE Saccharomyces cerevisiae (baker's yeast) 1655 ORGANISM Saccharomyces cerevisiae 1656 Eukaryota; Fungi; Ascomycota; Saccharomycotina; Saccharomycetes; 1657 Saccharomycetales; Saccharomycetaceae; Saccharomyces. 1658 REFERENCE 1 (bases 1 to 5028) 1659 AUTHORS Torpey,L.E., Gibbs,P.E., Nelson,J. and Lawrence,C.W. 1660 TITLE Cloning and sequence of REV7, a gene whose function is required for 1661 DNA damage-induced mutagenesis in Saccharomyces cerevisiae 1662 JOURNAL Yeast 10 (11), 1503-1509 (1994) 1663 PUBMED 7871890 1664 REFERENCE 2 (bases 1 to 5028) 1665 AUTHORS Roemer,T., Madden,K., Chang,J. and Snyder,M. 1666 TITLE Selection of axial growth sites in yeast requires Axl2p, a novel 1667 plasma membrane glycoprotein 1668 JOURNAL Genes Dev. 10 (7), 777-793 (1996) 1669 PUBMED 8846915 1670 REFERENCE 3 (bases 1 to 5028) 1671 AUTHORS Roemer,T. 1672 TITLE Direct Submission 1673 JOURNAL Submitted (22-FEB-1996) Terry Roemer, Biology, Yale University, New 1674 Haven, CT, USA 1675 FEATURES Location/Qualifiers 1676 source 1..5028 1677 /organism="Saccharomyces cerevisiae" 1678 /db_xref="taxon:4932" 1679 /chromosome="IX" 1680 /map="9" 1681 CDS <1..206 1682 /codon_start=3 1683 /product="TCP1-beta" 1684 /protein_id="AAA98665.1" 1685 /db_xref="GI:1293614" 1686 /translation="SSIYNGISTSGLDLNNGTIADMRQLGIVESYKLKRAVVSSASEA 1687 AEVLLRVDNIIRARPRTANRQHM" 1688 gene 687..3158 1689 /gene="AXL2" 1690 CDS 687..3158 1691 /gene="AXL2" 1692 /note="plasma membrane glycoprotein" 1693 /codon_start=1 1694 /function="required for axial budding pattern of S. 1695 cerevisiae" 1696 /product="Axl2p" 1697 /protein_id="AAA98666.1" 1698 /db_xref="GI:1293615" 1699 /translation="MTQLQISLLLTATISLLHLVVATPYEAYPIGKQYPPVARVNESF 1700 TFQISNDTYKSSVDKTAQITYNCFDLPSWLSFDSSSRTFSGEPSSDLLSDANTTLYFN 1701 VILEGTDSADSTSLNNTYQFVVTNRPSISLSSDFNLLALLKNYGYTNGKNALKLDPNE 1702 VFNVTFDRSMFTNEESIVSYYGRSQLYNAPLPNWLFFDSGELKFTGTAPVINSAIAPE 1703 TSYSFVIIATDIEGFSAVEVEFELVIGAHQLTTSIQNSLIINVTDTGNVSYDLPLNYV 1704 YLDDDPISSDKLGSINLLDAPDWVALDNATISGSVPDELLGKNSNPANFSVSIYDTYG 1705 DVIYFNFEVVSTTDLFAISSLPNINATRGEWFSYYFLPSQFTDYVNTNVSLEFTNSSQ 1706 DHDWVKFQSSNLTLAGEVPKNFDKLSLGLKANQGSQSQELYFNIIGMDSKITHSNHSA 1707 NATSTRSSHHSTSTSSYTSSTYTAKISSTSAAATSSAPAALPAANKTSSHNKKAVAIA 1708 CGVAIPLGVILVALICFLIFWRRRRENPDDENLPHAISGPDLNNPANKPNQENATPLN 1709 NPFDDDASSYDDTSIARRLAALNTLKLDNHSATESDISSVDEKRDSLSGMNTYNDQFQ 1710 SQSKEELLAKPPVQPPESPFFDPQNRSSSVYMDSEPAVNKSWRYTGNLSPVSDIVRDS 1711 YGSQKTVDTEKLFDLEAPEKEKRTSRDVTMSSLDPWNSNISPSPVRKSVTPSPYNVTK 1712 HRNRHLQNIQDSQSGKNGITPTTMSTSSSDDFVPVKDGENFCWVHSMEPDRRPSKKRL 1713 VDFSNKSNVNVGQVKDIHGRIPEML" 1714 gene complement(3300..4037) 1715 /gene="REV7" 1716 CDS complement(3300..4037) 1717 /gene="REV7" 1718 /codon_start=1 1719 /product="Rev7p" 1720 /protein_id="AAA98667.1" 1721 /db_xref="GI:1293616" 1722 /translation="MNRWVEKWLRVYLKCYINLILFYRNVYPPQSFDYTTYQSFNLPQ 1723 FVPINRHPALIDYIEELILDVLSKLTHVYRFSICIINKKNDLCIEKYVLDFSELQHVD 1724 KDDQIITETEVFDEFRSSLNSLIMHLEKLPKVNDDTITFEAVINAIELELGHKLDRNR 1725 RVDSLEEKAEIERDSNWVKCQEDENLPDNNGFQPPKIKLTSLVGSDVGPLIIHQFSEK 1726 LISGDDKILNGVYSQYEEGESIFGSLF" 1727 ORIGIN 1728 1 gatcctccat atacaacggt atctccacct caggtttaga tctcaacaac ggaaccattg 1729 61 ccgacatgag acagttaggt atcgtcgaga gttacaagct aaaacgagca gtagtcagct 1730 121 ctgcatctga agccgctgaa gttctactaa gggtggataa catcatccgt gcaagaccaa 1731 181 gaaccgccaa tagacaacat atgtaacata tttaggatat acctcgaaaa taataaaccg 1732 241 ccacactgtc attattataa ttagaaacag aacgcaaaaa ttatccacta tataattcaa 1733 301 agacgcgaaa aaaaaagaac aacgcgtcat agaacttttg gcaattcgcg tcacaaataa 1734 361 attttggcaa cttatgtttc ctcttcgagc agtactcgag ccctgtctca agaatgtaat 1735 421 aatacccatc gtaggtatgg ttaaagatag catctccaca acctcaaagc tccttgccga 1736 481 gagtcgccct cctttgtcga gtaattttca cttttcatat gagaacttat tttcttattc 1737 541 tttactctca catcctgtag tgattgacac tgcaacagcc accatcacta gaagaacaga 1738 601 acaattactt aatagaaaaa ttatatcttc ctcgaaacga tttcctgctt ccaacatcta 1739 661 cgtatatcaa gaagcattca cttaccatga cacagcttca gatttcatta ttgctgacag 1740 721 ctactatatc actactccat ctagtagtgg ccacgcccta tgaggcatat cctatcggaa 1741 781 aacaataccc cccagtggca agagtcaatg aatcgtttac atttcaaatt tccaatgata 1742 841 cctataaatc gtctgtagac aagacagctc aaataacata caattgcttc gacttaccga 1743 901 gctggctttc gtttgactct agttctagaa cgttctcagg tgaaccttct tctgacttac 1744 961 tatctgatgc gaacaccacg ttgtatttca atgtaatact cgagggtacg gactctgccg 1745 1021 acagcacgtc tttgaacaat acataccaat ttgttgttac aaaccgtcca tccatctcgc 1746 1081 tatcgtcaga tttcaatcta ttggcgttgt taaaaaacta tggttatact aacggcaaaa 1747 1141 acgctctgaa actagatcct aatgaagtct tcaacgtgac ttttgaccgt tcaatgttca 1748 1201 ctaacgaaga atccattgtg tcgtattacg gacgttctca gttgtataat gcgccgttac 1749 1261 ccaattggct gttcttcgat tctggcgagt tgaagtttac tgggacggca ccggtgataa 1750 1321 actcggcgat tgctccagaa acaagctaca gttttgtcat catcgctaca gacattgaag 1751 1381 gattttctgc cgttgaggta gaattcgaat tagtcatcgg ggctcaccag ttaactacct 1752 1441 ctattcaaaa tagtttgata atcaacgtta ctgacacagg taacgtttca tatgacttac 1753 1501 ctctaaacta tgtttatctc gatgacgatc ctatttcttc tgataaattg ggttctataa 1754 1561 acttattgga tgctccagac tgggtggcat tagataatgc taccatttcc gggtctgtcc 1755 1621 cagatgaatt actcggtaag aactccaatc ctgccaattt ttctgtgtcc atttatgata 1756 1681 cttatggtga tgtgatttat ttcaacttcg aagttgtctc cacaacggat ttgtttgcca 1757 1741 ttagttctct tcccaatatt aacgctacaa ggggtgaatg gttctcctac tattttttgc 1758 1801 cttctcagtt tacagactac gtgaatacaa acgtttcatt agagtttact aattcaagcc 1759 1861 aagaccatga ctgggtgaaa ttccaatcat ctaatttaac attagctgga gaagtgccca 1760 1921 agaatttcga caagctttca ttaggtttga aagcgaacca aggttcacaa tctcaagagc 1761 1981 tatattttaa catcattggc atggattcaa agataactca ctcaaaccac agtgcgaatg 1762 2041 caacgtccac aagaagttct caccactcca cctcaacaag ttcttacaca tcttctactt 1763 2101 acactgcaaa aatttcttct acctccgctg ctgctacttc ttctgctcca gcagcgctgc 1764 2161 cagcagccaa taaaacttca tctcacaata aaaaagcagt agcaattgcg tgcggtgttg 1765 2221 ctatcccatt aggcgttatc ctagtagctc tcatttgctt cctaatattc tggagacgca 1766 2281 gaagggaaaa tccagacgat gaaaacttac cgcatgctat tagtggacct gatttgaata 1767 2341 atcctgcaaa taaaccaaat caagaaaacg ctacaccttt gaacaacccc tttgatgatg 1768 2401 atgcttcctc gtacgatgat acttcaatag caagaagatt ggctgctttg aacactttga 1769 2461 aattggataa ccactctgcc actgaatctg atatttccag cgtggatgaa aagagagatt 1770 2521 ctctatcagg tatgaataca tacaatgatc agttccaatc ccaaagtaaa gaagaattat 1771 2581 tagcaaaacc cccagtacag cctccagaga gcccgttctt tgacccacag aataggtctt 1772 2641 cttctgtgta tatggatagt gaaccagcag taaataaatc ctggcgatat actggcaacc 1773 2701 tgtcaccagt ctctgatatt gtcagagaca gttacggatc acaaaaaact gttgatacag 1774 2761 aaaaactttt cgatttagaa gcaccagaga aggaaaaacg tacgtcaagg gatgtcacta 1775 2821 tgtcttcact ggacccttgg aacagcaata ttagcccttc tcccgtaaga aaatcagtaa 1776 2881 caccatcacc atataacgta acgaagcatc gtaaccgcca cttacaaaat attcaagact 1777 2941 ctcaaagcgg taaaaacgga atcactccca caacaatgtc aacttcatct tctgacgatt 1778 3001 ttgttccggt taaagatggt gaaaattttt gctgggtcca tagcatggaa ccagacagaa 1779 3061 gaccaagtaa gaaaaggtta gtagattttt caaataagag taatgtcaat gttggtcaag 1780 3121 ttaaggacat tcacggacgc atcccagaaa tgctgtgatt atacgcaacg atattttgct 1781 3181 taattttatt ttcctgtttt attttttatt agtggtttac agatacccta tattttattt 1782 3241 agtttttata cttagagaca tttaatttta attccattct tcaaatttca tttttgcact 1783 3301 taaaacaaag atccaaaaat gctctcgccc tcttcatatt gagaatacac tccattcaaa 1784 3361 attttgtcgt caccgctgat taatttttca ctaaactgat gaataatcaa aggccccacg 1785 3421 tcagaaccga ctaaagaagt gagttttatt ttaggaggtt gaaaaccatt attgtctggt 1786 3481 aaattttcat cttcttgaca tttaacccag tttgaatccc tttcaatttc tgctttttcc 1787 3541 tccaaactat cgaccctcct gtttctgtcc aacttatgtc ctagttccaa ttcgatcgca 1788 3601 ttaataactg cttcaaatgt tattgtgtca tcgttgactt taggtaattt ctccaaatgc 1789 3661 ataatcaaac tatttaagga agatcggaat tcgtcgaaca cttcagtttc cgtaatgatc 1790 3721 tgatcgtctt tatccacatg ttgtaattca ctaaaatcta aaacgtattt ttcaatgcat 1791 3781 aaatcgttct ttttattaat aatgcagatg gaaaatctgt aaacgtgcgt taatttagaa 1792 3841 agaacatcca gtataagttc ttctatatag tcaattaaag caggatgcct attaatggga 1793 3901 acgaactgcg gcaagttgaa tgactggtaa gtagtgtagt cgaatgactg aggtgggtat 1794 3961 acatttctat aaaataaaat caaattaatg tagcatttta agtataccct cagccacttc 1795 4021 tctacccatc tattcataaa gctgacgcaa cgattactat tttttttttc ttcttggatc 1796 4081 tcagtcgtcg caaaaacgta taccttcttt ttccgacctt ttttttagct ttctggaaaa 1797 4141 gtttatatta gttaaacagg gtctagtctt agtgtgaaag ctagtggttt cgattgactg 1798 4201 atattaagaa agtggaaatt aaattagtag tgtagacgta tatgcatatg tatttctcgc 1799 4261 ctgtttatgt ttctacgtac ttttgattta tagcaagggg aaaagaaata catactattt 1800 4321 tttggtaaag gtgaaagcat aatgtaaaag ctagaataaa atggacgaaa taaagagagg 1801 4381 cttagttcat cttttttcca aaaagcaccc aatgataata actaaaatga aaaggatttg 1802 4441 ccatctgtca gcaacatcag ttgtgtgagc aataataaaa tcatcacctc cgttgccttt 1803 4501 agcgcgtttg tcgtttgtat cttccgtaat tttagtctta tcaatgggaa tcataaattt 1804 4561 tccaatgaat tagcaatttc gtccaattct ttttgagctt cttcatattt gctttggaat 1805 4621 tcttcgcact tcttttccca ttcatctctt tcttcttcca aagcaacgat ccttctaccc 1806 4681 atttgctcag agttcaaatc ggcctctttc agtttatcca ttgcttcctt cagtttggct 1807 4741 tcactgtctt ctagctgttg ttctagatcc tggtttttct tggtgtagtt ctcattatta 1808 4801 gatctcaagt tattggagtc ttcagccaat tgctttgtat cagacaattg actctctaac 1809 4861 ttctccactt cactgtcgag ttgctcgttt ttagcggaca aagatttaat ctcgttttct 1810 4921 ttttcagtgt tagattgctc taattctttg agctgttctc tcagctcctc atatttttct 1811 4981 tgccatgact cagattctaa ttttaagcta ttcaatttct ctttgatc 1812 //""" 1813 1814 # GenBank format protein (aka GenPept) file from: 1815 # http://www.molecularevolution.org/resources/fileformats/ 1816 gbk_example2 = \ 1817 """LOCUS AAD51968 143 aa linear BCT 21-AUG-2001 1818 DEFINITION transcriptional regulator RovA [Yersinia enterocolitica]. 1819 ACCESSION AAD51968 1820 VERSION AAD51968.1 GI:5805369 1821 DBSOURCE locus AF171097 accession AF171097.1 1822 KEYWORDS . 1823 SOURCE Yersinia enterocolitica 1824 ORGANISM Yersinia enterocolitica 1825 Bacteria; Proteobacteria; Gammaproteobacteria; Enterobacteriales; 1826 Enterobacteriaceae; Yersinia. 1827 REFERENCE 1 (residues 1 to 143) 1828 AUTHORS Revell,P.A. and Miller,V.L. 1829 TITLE A chromosomally encoded regulator is required for expression of the 1830 Yersinia enterocolitica inv gene and for virulence 1831 JOURNAL Mol. Microbiol. 35 (3), 677-685 (2000) 1832 MEDLINE 20138369 1833 PUBMED 10672189 1834 REFERENCE 2 (residues 1 to 143) 1835 AUTHORS Revell,P.A. and Miller,V.L. 1836 TITLE Direct Submission 1837 JOURNAL Submitted (22-JUL-1999) Molecular Microbiology, Washington 1838 University School of Medicine, Campus Box 8230, 660 South Euclid, 1839 St. Louis, MO 63110, USA 1840 COMMENT Method: conceptual translation. 1841 FEATURES Location/Qualifiers 1842 source 1..143 1843 /organism="Yersinia enterocolitica" 1844 /mol_type="unassigned DNA" 1845 /strain="JB580v" 1846 /serotype="O:8" 1847 /db_xref="taxon:630" 1848 Protein 1..143 1849 /product="transcriptional regulator RovA" 1850 /name="regulates inv expression" 1851 CDS 1..143 1852 /gene="rovA" 1853 /coded_by="AF171097.1:380..811" 1854 /note="regulator of virulence" 1855 /transl_table=11 1856 ORIGIN 1857 1 mestlgsdla rlvrvwrali dhrlkplelt qthwvtlhni nrlppeqsqi qlakaigieq 1858 61 pslvrtldql eekglitrht candrrakri klteqsspii eqvdgvicst rkeilggisp 1859 121 deiellsgli dklerniiql qsk 1860 // 1861 """ 1862 1863 embl_example = """ID X56734; SV 1; linear; mRNA; STD; PLN; 1859 BP. 1864 XX 1865 AC X56734; S46826; 1866 XX 1867 DT 12-SEP-1991 (Rel. 29, Created) 1868 DT 25-NOV-2005 (Rel. 85, Last updated, Version 11) 1869 XX 1870 DE Trifolium repens mRNA for non-cyanogenic beta-glucosidase 1871 XX 1872 KW beta-glucosidase. 1873 XX 1874 OS Trifolium repens (white clover) 1875 OC Eukaryota; Viridiplantae; Streptophyta; Embryophyta; Tracheophyta; 1876 OC Spermatophyta; Magnoliophyta; eudicotyledons; core eudicotyledons; rosids; 1877 OC eurosids I; Fabales; Fabaceae; Papilionoideae; Trifolieae; Trifolium. 1878 XX 1879 RN [5] 1880 RP 1-1859 1881 RX PUBMED; 1907511. 1882 RA Oxtoby E., Dunn M.A., Pancoro A., Hughes M.A.; 1883 RT "Nucleotide and derived amino acid sequence of the cyanogenic 1884 RT beta-glucosidase (linamarase) from white clover (Trifolium repens L.)"; 1885 RL Plant Mol. Biol. 17(2):209-219(1991). 1886 XX 1887 RN [6] 1888 RP 1-1859 1889 RA Hughes M.A.; 1890 RT ; 1891 RL Submitted (19-NOV-1990) to the EMBL/GenBank/DDBJ databases. 1892 RL Hughes M.A., University of Newcastle Upon Tyne, Medical School, Newcastle 1893 RL Upon Tyne, NE2 4HH, UK 1894 XX 1895 FH Key Location/Qualifiers 1896 FH 1897 FT source 1..1859 1898 FT /organism="Trifolium repens" 1899 FT /mol_type="mRNA" 1900 FT /clone_lib="lambda gt10" 1901 FT /clone="TRE361" 1902 FT /tissue_type="leaves" 1903 FT /db_xref="taxon:3899" 1904 FT CDS 14..1495 1905 FT /product="beta-glucosidase" 1906 FT /EC_number="3.2.1.21" 1907 FT /note="non-cyanogenic" 1908 FT /db_xref="GOA:P26204" 1909 FT /db_xref="InterPro:IPR001360" 1910 FT /db_xref="InterPro:IPR013781" 1911 FT /db_xref="UniProtKB/Swiss-Prot:P26204" 1912 FT /protein_id="CAA40058.1" 1913 FT /translation="MDFIVAIFALFVISSFTITSTNAVEASTLLDIGNLSRSSFPRGFI 1914 FT FGAGSSAYQFEGAVNEGGRGPSIWDTFTHKYPEKIRDGSNADITVDQYHRYKEDVGIMK 1915 FT DQNMDSYRFSISWPRILPKGKLSGGINHEGIKYYNNLINELLANGIQPFVTLFHWDLPQ 1916 FT VLEDEYGGFLNSGVINDFRDYTDLCFKEFGDRVRYWSTLNEPWVFSNSGYALGTNAPGR 1917 FT CSASNVAKPGDSGTGPYIVTHNQILAHAEAVHVYKTKYQAYQKGKIGITLVSNWLMPLD 1918 FT DNSIPDIKAAERSLDFQFGLFMEQLTTGDYSKSMRRIVKNRLPKFSKFESSLVNGSFDF 1919 FT IGINYYSSSYISNAPSHGNAKPSYSTNPMTNISFEKHGIPLGPRAASIWIYVYPYMFIQ 1920 FT EDFEIFCYILKINITILQFSITENGMNEFNDATLPVEEALLNTYRIDYYYRHLYYIRSA 1921 FT IRAGSNVKGFYAWSFLDCNEWFAGFTVRFGLNFVD" 1922 FT mRNA 1..1859 1923 FT /experiment="experimental evidence, no additional details 1924 FT recorded" 1925 XX 1926 SQ Sequence 1859 BP; 609 A; 314 C; 355 G; 581 T; 0 other; 1927 aaacaaacca aatatggatt ttattgtagc catatttgct ctgtttgtta ttagctcatt 60 1928 cacaattact tccacaaatg cagttgaagc ttctactctt cttgacatag gtaacctgag 120 1929 tcggagcagt tttcctcgtg gcttcatctt tggtgctgga tcttcagcat accaatttga 180 1930 aggtgcagta aacgaaggcg gtagaggacc aagtatttgg gataccttca cccataaata 240 1931 tccagaaaaa ataagggatg gaagcaatgc agacatcacg gttgaccaat atcaccgcta 300 1932 caaggaagat gttgggatta tgaaggatca aaatatggat tcgtatagat tctcaatctc 360 1933 ttggccaaga atactcccaa agggaaagtt gagcggaggc ataaatcacg aaggaatcaa 420 1934 atattacaac aaccttatca acgaactatt ggctaacggt atacaaccat ttgtaactct 480 1935 ttttcattgg gatcttcccc aagtcttaga agatgagtat ggtggtttct taaactccgg 540 1936 tgtaataaat gattttcgag actatacgga tctttgcttc aaggaatttg gagatagagt 600 1937 gaggtattgg agtactctaa atgagccatg ggtgtttagc aattctggat atgcactagg 660 1938 aacaaatgca ccaggtcgat gttcggcctc caacgtggcc aagcctggtg attctggaac 720 1939 aggaccttat atagttacac acaatcaaat tcttgctcat gcagaagctg tacatgtgta 780 1940 taagactaaa taccaggcat atcaaaaggg aaagataggc ataacgttgg tatctaactg 840 1941 gttaatgcca cttgatgata atagcatacc agatataaag gctgccgaga gatcacttga 900 1942 cttccaattt ggattgttta tggaacaatt aacaacagga gattattcta agagcatgcg 960 1943 gcgtatagtt aaaaaccgat tacctaagtt ctcaaaattc gaatcaagcc tagtgaatgg 1020 1944 ttcatttgat tttattggta taaactatta ctcttctagt tatattagca atgccccttc 1080 1945 acatggcaat gccaaaccca gttactcaac aaatcctatg accaatattt catttgaaaa 1140 1946 acatgggata cccttaggtc caagggctgc ttcaatttgg atatatgttt atccatatat 1200 1947 gtttatccaa gaggacttcg agatcttttg ttacatatta aaaataaata taacaatcct 1260 1948 gcaattttca atcactgaaa atggtatgaa tgaattcaac gatgcaacac ttccagtaga 1320 1949 agaagctctt ttgaatactt acagaattga ttactattac cgtcacttat actacattcg 1380 1950 ttctgcaatc agggctggct caaatgtgaa gggtttttac gcatggtcat ttttggactg 1440 1951 taatgaatgg tttgcaggct ttactgttcg ttttggatta aactttgtag attagaaaga 1500 1952 tggattaaaa aggtacccta agctttctgc ccaatggtac aagaactttc tcaaaagaaa 1560 1953 ctagctagta ttattaaaag aactttgtag tagattacag tacatcgttt gaagttgagt 1620 1954 tggtgcacct aattaaataa aagaggttac tcttaacata tttttaggcc attcgttgtg 1680 1955 aagttgttag gctgttattt ctattatact atgttgtagt aataagtgca ttgttgtacc 1740 1956 agaagctatg atcataacta taggttgatc cttcatgtat cagtttgatg ttgagaatac 1800 1957 tttgaattaa aagtcttttt ttattttttt aaaaaaaaaa aaaaaaaaaa aaaaaaaaa 1859 1958 // 1959 """ 1960 1961 print("GenBank CDS Iteration") 1962 print("=====================") 1963 1964 g = GenBankScanner() 1965 for record in g.parse_cds_features(StringIO(gbk_example)): 1966 print(record) 1967 1968 g = GenBankScanner() 1969 for record in g.parse_cds_features(StringIO(gbk_example2), 1970 tags2id=('gene', 'locus_tag', 'product')): 1971 print(record) 1972 1973 g = GenBankScanner() 1974 for record in g.parse_cds_features(StringIO(gbk_example + "\n" + gbk_example2), 1975 tags2id=('gene', 'locus_tag', 'product')): 1976 print(record) 1977 1978 print("") 1979 print("GenBank Iteration") 1980 print("=================") 1981 g = GenBankScanner() 1982 for record in g.parse_records(StringIO(gbk_example), do_features=False): 1983 print("%s %s %s" % (record.id, record.name, record.description)) 1984 print(record.seq) 1985 1986 g = GenBankScanner() 1987 for record in g.parse_records(StringIO(gbk_example), do_features=True): 1988 print("%s %s %s" % (record.id, record.name, record.description)) 1989 print(record.seq) 1990 1991 g = GenBankScanner() 1992 for record in g.parse_records(StringIO(gbk_example2), do_features=False): 1993 print("%s %s %s" % (record.id, record.name, record.description)) 1994 print(record.seq) 1995 1996 g = GenBankScanner() 1997 for record in g.parse_records(StringIO(gbk_example2), do_features=True): 1998 print("%s %s %s" % (record.id, record.name, record.description)) 1999 print(record.seq) 2000 2001 print("") 2002 print("EMBL CDS Iteration") 2003 print("==================") 2004 2005 e = EmblScanner() 2006 for record in e.parse_cds_features(StringIO(embl_example)): 2007 print(record) 2008 2009 print("") 2010 print("EMBL Iteration") 2011 print("==============") 2012 e = EmblScanner() 2013 for record in e.parse_records(StringIO(embl_example), do_features=True): 2014 print("%s %s %s" % (record.id, record.name, record.description)) 2015 print(record.seq) 2016