This file is indexed.

/usr/share/pyshared/cogent/parse/genbank.py is in python-cogent 1.5.3-2.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
#!/usr/bin/env python
from cogent.parse.record import FieldWrapper
from cogent.parse.record_finder import DelimitedRecordFinder, \
    LabeledRecordFinder
from cogent.core.genetic_code import GeneticCodes
from string import maketrans, strip, rstrip
from cogent.core.moltype import PROTEIN, DNA, ASCII
from cogent.core.annotation import Feature
from cogent.core.info import Info

__author__ = "Rob Knight"
__copyright__ = "Copyright 2007-2012, The Cogent Project"
__credits__ = ["Rob Knight", "Peter Maxwell", "Matthew Wakefield",
                    "Gavin Huttley"]
__license__ = "GPL"
__version__ = "1.5.3"
__maintainer__ = "Rob Knight"
__email__ = "rob@spot.colorado.edu"
__status__ = "Production"

all_chars = maketrans('','')
dna_lc =     'utacgrywsmkbdhvn'
dna_lc_cmp = 'aatgcyrwskmvhdbn'
dna_trans = maketrans(dna_lc+dna_lc.upper(),dna_lc_cmp+dna_lc_cmp.upper())
rna_lc =     'utacgrywsmkbdhvn'
rna_lc_cmp = 'aaugcyrwskmvhdbn'
rna_trans = maketrans(rna_lc+rna_lc.upper(),rna_lc_cmp+rna_lc_cmp.upper())

locus_fields = [None, 'locus','length', None, 'mol_type','topology','db','date']
_locus_parser = FieldWrapper(locus_fields)

#need to turn off line stripping, because whitespace is significant
GbFinder = DelimitedRecordFinder('//', constructor=rstrip)

class PartialRecordError(Exception):
    pass

def parse_locus(line):
    """Parses a locus line, including conversion of Length to an int.
    
    WARNING: Gives incorrect results on legacy records that omit the topology. 
    All records spot-checked on 8/30/05 had been updated to include the topology
    even when prior versions omitted it.
    """
    result = _locus_parser(line)
    try:
        result['length'] = int(result['length'])
    except KeyError, e:
        raise PartialRecordError, e

    if None in result:
        del result[None]
    return result

def parse_single_line(line):
    """Generic parser: splits off the label, and return the rest."""
    label, data = line.split(None, 1)
    return data.rstrip()

def indent_splitter(lines):
    """Yields the lines whenever it hits a line with same indent level as first.
    """
    first_line = True
    curr = []
    for line in lines:
        #skip blank lines
        line = line.rstrip()
        if not line:
            continue
        #need to figure out indent if first line
        if first_line:
            indent = len(line) - len(line.lstrip())
            curr.append(line)
            first_line = False
        elif len(line) > indent and line[indent].isspace():
            curr.append(line)
        else:   #got a line that doesn't match the indent
            yield curr
            curr = [line]
    if curr:
        yield curr

def parse_sequence(lines, constructor=''.join):
    """Parses a GenBank sequence block. Doesn't care about ORIGIN line."""
    result = []
    for i in lines:
        if i.startswith('ORIGIN'):
            continue
        result.append(i.translate(all_chars, '0123456789 \t\n\r/'))
    return constructor(result)

def block_consolidator(lines):
    """Takes block with label and multiline data, and returns (label, [data]).
    
    [data] will be list of lines of data, including first line w/o label.
    """
    data = []
    first = True
    label = None
    for line in lines:
        if first:   #find label
            line = line.split(None, 1)
            if len(line) == 2:
                label, curr = line
            else:
                label = line[0]
                curr = ""
            data.append(curr)
            first = False
        else:
            data.append(line)
    return label, data

def parse_organism(lines):
    """Takes ORGANISM block. Returns organism, [taxonomy].
    
    NOTE: Adds species to end of taxonomy if identifiable.
    """
    label, data = block_consolidator(lines)
    #get 'species'
    species = data[0].strip()
    #get rest of taxonomy
    taxonomy = ' '.join(data[1:])
    #normalize whitespace, including deleting newlines
    taxonomy = ' '.join(taxonomy.split())
    #separate by semicolons
    taxa = map(strip, taxonomy.split(';'))  #get rid of leading/trailing spaces
    #delete trailing period if present
    last = taxa[-1]
    if last.endswith('.'):
        taxa[-1] = last[:-1]
    return species, taxa

def is_feature_component_start(line):
    """Checks if a line starts with '/', ignoring whitespace."""
    return line.lstrip().startswith('/')

feature_component_iterator = LabeledRecordFinder(is_feature_component_start)

_join_with_empty = dict.fromkeys(['translation'])
_leave_as_lines = {}

def parse_feature(lines):
    """Parses a feature. Doesn't handle subfeatures.
    
    Returns dict containing:
    'type': source, gene, CDS, etc.
    'location': unparsed location string
    ...then, key-value pairs for each annotation, 
        e.g. '/gene="MNBH"' -> {'gene':['MNBH']} (i.e. quotes stripped)
    All relations are assumed 'to many', and order will be preserved.
    """
    result = {}
    type_, data = block_consolidator(lines)
    result['type'] = type_
    location = []
    found_feature = False
    for curr_line_idx, line in enumerate(data):
        if line.lstrip().startswith('/'):
            found_feature = True
            break
        else:
            location.append(line)
    result['raw_location'] = location
    try:
        result['location'] = \
            parse_location_line(location_line_tokenizer(location))
    except (TypeError, ValueError):
        result['location'] = None
    if not found_feature:
        return result
    fci = feature_component_iterator
    for feature_component in fci(data[curr_line_idx:]):
        first = feature_component[0].lstrip()[1:]   #remove leading space, '/'
        try:
            label, first_line = first.split('=', 1)
        except ValueError:   #sometimes not delimited by =
            label, first_line = first, ''
        #chop off leading quote if appropriate
        if first_line.startswith('"'):
            first_line = first_line[1:]
        remainder = [first_line] + feature_component[1:]
        #chop off trailing quote, if appropriate
        last_line = remainder[-1].rstrip()
        if last_line.endswith('"'):
            remainder[-1] = last_line[:-1]
        if label in _join_with_empty:
            curr_data = ''.join(map(strip, remainder))
        elif label in _leave_as_lines:
            curr_data = remainder
        else:
            curr_data = ' '.join(map(strip, remainder))
        if label.lower() == 'type':
            # some source features have /type=...
            label = 'type_field'
        if label not in result:
            result[label.lower()] = []

        result[label.lower()].append(curr_data)
    return result

def location_line_tokenizer(lines):
    """Tokenizes location lines into spans, joins and complements."""
    curr = []
    text = ' '.join(map(strip, lines))
    for char in text:
        if char == '(':
            yield ''.join(curr).strip() + char
            curr = []
        elif char == ')':
            if curr:
                yield ''.join(curr).strip()
            yield char
            curr = []
        elif char == ',':
            if curr:
                yield ''.join(curr).strip()
            yield ','
            curr = []
        else:
            curr.append(char)
    if curr:
        yield ''.join(curr).strip()

def parse_simple_location_segment(segment):
    """Parses location segment of form a..b or a, incl. '<' and '>'."""
    first_ambiguity, second_ambiguity = None, None
    if '..' in segment:
        first, second = segment.split('..')
        if not first[0].isdigit():
            first_ambiguity = first[0]
            first = long(first[1:])
        else:
            first = long(first)
        if not second[0].isdigit():
            second_ambiguity = second[0]
            second = long(second[1:])
        else:
            second = long(second)
        
        return Location([Location(first, Ambiguity=first_ambiguity), \
            Location(second, Ambiguity=second_ambiguity)])
    else:
        if not segment[0].isdigit():
            first_ambiguity = segment[0]
            segment = segment[1:]
        return Location(long(segment), Ambiguity=first_ambiguity)

def parse_location_line(tokens, parser=parse_simple_location_segment):
    """Parses location line tokens into location list."""
    stack = []
    curr = stack
    for t in tokens:
        if t .endswith('('):
            new = [curr, t]
            curr.append(new)
            curr = new
        elif t == ',':  #ignore
            continue
        elif t == ')':
            parent, type_ = curr[:2]
            children = curr[2:]
            if type_ == 'complement(':
                children.reverse()
                for c in children:
                    c.Strand *= -1
            curr_index = parent.index(curr)
            del parent[curr_index]
            parent[curr_index:curr_index] = children[:]
            curr = parent
        else:
            curr.append(parser(t))
    return LocationList(stack)
            
class Location(object):
    """GenBank location object. Integer, or low, high, or 2-base bound.
   
    data must either be a long, an object that can be coerced to a long, or a 
        sequence of two BasePosition objects. It can _not_ be two numbers.
    Ambiguity should be None (the default), '>', or '<'.
    IsBetween should be False (the default), or True.
    IsBounds should be False(the default, indicates range), or True.
    Accession should be an accession, or None (default).
    Db should be a database identifier, or None (default).
    Strand should be 1 (forward, default) or -1 (reverse).
    
    WARNING: This Location will allow you to do things that can't happen in
    GenBank, such as having a start and stop that aren't from the same
    accession. No validation is performed to prevent this. All reasonable
    cases should work.
    
    WARNING: Coordinates are based on 1, not 0, as in GenBank format.
    """
    def __init__(self, data, Ambiguity=None, IsBetween=False, IsBounds=False, \
            Accession=None, Db=None, Strand=1):
        """Returns new LocalLocation object."""
        try:
            data = long(data)
        except TypeError:
            pass    #assume was two Location objects.
        self._data = data
        self.Ambiguity = Ambiguity
        self.IsBetween = IsBetween
        self.IsBounds = IsBounds
        self.Accession = Accession
        self.Db = Db
        self.Strand = Strand
    
    def __str__(self):
        """Returns self in string format.
        
        WARNING: More permissive than GenBank's Backus-Naur form allows. If
        you abuse this object, you'll get results that aren't valid GenBank
        locations.
        """
        if self.IsBetween:  #between two bases
            try:
                first, last = self._data
                curr = '%s^%s' % (first, last)
            except TypeError:   #only one base? must be this or the next
                curr = '%s^%s' % (first, first+1)
        else:   #not self.IsBetween
            try:
                data = long(self._data)
                #if the above line succeeds, we've got a single item
                if self.Ambiguity:
                    curr = self.Ambiguity + str(data)
                else:
                    curr = str(data)
            except TypeError:
                #if long conversion failed, should have two LocalLocation objects
                first, last = self._data
                if self.IsBounds:
                    curr = '(%s%s%s)' % (first, '.', last)
                else:
                    curr = '%s%s%s' % (first, '..', last)
        #check if we need to add on the accession and database
        if self.Accession:
            curr = self.Accession + ':' + curr
            #we're only going to add the Db if we got an accession
            if self.Db:
                curr = self.Db + '::' + curr
        #check if it's complemented
        if self.Strand == -1:
            curr = 'complement(%s)' % curr
        return curr
    
    def isAmbiguous(self):
        """Returns True if ambiguous (single-base ambiguity or two locations.)
        """
        if self.Ambiguity:
            return True
        try:
            iter(self._data)
            return True
        except:
            return False
    
    def first(self):
        """Returns first base self could be."""
        try:
            return long(self._data)
        except TypeError:
            return self._data[0].first()
    
    def last(self):
        """Returns last base self could be."""
        try:
            return long(self._data)
        except TypeError:
            return self._data[-1].last()

class LocationList(list):
    """List of Location objects.
    
    WARNING: Coordinates are based on 1, not 0, to match GenBank format.
    """
    BIGNUM = 1e300
    def first(self):
        """Returns first base of self."""
        curr = self.BIGNUM
        for i in self:
            first = i.first()
            if curr > first:
                curr = first
        return curr
    
    def last(self):
        """Returns last base of self."""
        curr = 0
        for i in self:
            last = i.last()
            if last > curr:
                curr = last
        return curr
    
    def strand(self):
        """Returns strand of components: 1=forward, -1=reverse, 0=both
        """
        curr = {}
        for i in self:
            curr[i.Strand] = 1
        if len(curr) >= 2:  #found stuff on both strands
            return 0
        else:
            return curr.keys()[0]
    
    def __str__(self):
        """Returns (normalized) string representation of self."""
        if len(self) == 0:
            return ''
        elif len(self) == 1:
            return str(self[0])
        else:
            return 'join(' + ','.join(map(str, self)) + ')'
    
    def extract(self, sequence, trans_table=dna_trans):
        """Extracts pieces of self from sequence."""
        result = []
        for i in self:
            first, last = i.first(), i.last() + 1   #inclusive, not exclusive
            #translate to 0-based indices and check if it wraps around
            if first < last:
                curr = sequence[first-1:last-1]
            else:
                curr = sequence[first-1:]+sequence[:last-1]
            #reverse-complement if necessary
            if i.Strand == -1:
                curr = curr.translate(trans_table)[::-1]
            result.append(curr)
        return ''.join(result)
    

def parse_feature_table(lines):
    """Simple parser for feature table. Assumes starts with FEATURES line."""
    if not lines:
        return []
    if lines[0].startswith('FEATURES'):
        lines = lines[1:]
    return [parse_feature(f) for f in indent_splitter(lines)]

reference_label_marker = ' ' * 11
reference_field_finder = LabeledRecordFinder(lambda x: \
    not x.startswith(reference_label_marker), constructor=None)

def parse_reference(lines):
    """Simple parser for single reference."""
    result = {}
    for field in reference_field_finder(lines):
        label, data = block_consolidator(field)
        result[label.lower()] = ' '.join(map(strip, data))
    return result

def parse_source(lines):
    """Simple parser for source fields."""
    result = {}
    all_lines = list(lines)
    source_field = reference_field_finder(all_lines).next()
    label, data = block_consolidator(source_field)
    result[label.lower()] = ' '.join(map(strip, data))
    source_length = len(source_field)
    species, taxonomy = parse_organism(lines[source_length:])
    result['species'] = species
    result['taxonomy'] = taxonomy
    return result

#adaptors to update curr with data from each parser
def locus_adaptor(lines, curr):
    curr.update(parse_locus(lines[0]))
def source_adaptor(lines, curr):
    curr.update(parse_source(lines))
def ref_adaptor(lines, curr):
    if 'references' not in curr:
        curr['references'] = []
    curr['references'].append(parse_reference(lines))
def feature_table_adaptor(lines, curr):
    if 'features' not in curr:
        curr['features'] = []
    curr['features'].extend(parse_feature_table(lines))
def sequence_adaptor(lines, curr):
    curr['sequence'] = parse_sequence(lines)
def generic_adaptor(lines, curr):
    label, data = block_consolidator(lines)
    curr[label.lower()] = ' '.join(map(strip, lines))

handlers = {
    'LOCUS': locus_adaptor,
    'SOURCE': source_adaptor,
    'REFERENCE': ref_adaptor,
    'FEATURES': feature_table_adaptor,
    'ORIGIN': sequence_adaptor,
    '//': lambda lines, curr: None,
    '?': lambda lines, curr: None
}

def MinimalGenbankParser(lines, handlers=handlers,\
    default_handler=generic_adaptor):
    for rec in GbFinder(lines):
        curr = {}
        bad_record = False
        for field in indent_splitter(rec):
            first_word = field[0].split(None, 1)[0]
            handler = handlers.get(first_word, default_handler)

            try:
                handler(field, curr)
            except:
                bad_record = True
                break

        if not bad_record:
            yield curr

def parse_location_segment(location_segment):
    """Parses a location segment into its component pieces.
    
    Known possibilities:
    http://www.ebi.ac.uk/embl/Documentation/FT_definitions/feature_table.html
    
    467             single base
    a..b            range from a to b, including a and b
    <a              strictly before a
    >a              strictly after a
    (a.b)           a single base between a and b, inclusive
    a^b             a site between two adjacent bases between a and b
    accession:a     a occurs in accession, not in the current sequence
    db::accession:a a occurrs in accession in db, not in the current sequence
    """
    s = location_segment    #save some typing...
    lsp = parse_location_segment
    #check if it's a range
    if '..' in s:
        first, second = s.split('..')
        return Location([lsp(first), lsp(second)])
    #check if it's between two adjacent bases
    elif '^' in s:
        first, second = s.split('^')
        return Location([lsp(first), lsp(second)],IsBetween=True)
    #check if it's a single base reference -- but don't be fooled by accessions!
    elif '.' in s and s.startswith('(') and s.endswith(')'):
        first, second = s.split('.')
        return Location([lsp(first[1:]), lsp(second[:-1])])

def parse_location_atom(location_atom):
    """Parses a location atom, supposed to be a single-base position."""
    a = location_atom
    if a.startswith('<') or a.startswith('>'):   #fuzzy
        position = long(a[1:])
        return Location(position, Ambiguity = a[0])
    #otherwise, should just be an integer
    return Location(long(a))

wanted_types = dict.fromkeys(['CDS'])

def extract_nt_prot_seqs(rec, wanted=wanted_types):
    """Extracts nucleotide seqs, and, where possible, protein seqs, from recs."""
    rec_seq = rec['sequence']
    for f in rec['features']:
        if f['type'] not in wanted:
            continue
        translation = f['translation'][0]
        raw_seq = f['location'].extract(rec_seq)
        print raw_seq
        seq = raw_seq[long(f['codon_start'][0])-1:]
        print 'dt:', translation
        print 'ct:', GeneticCodes[f.get('transl_table', '1')[0]].translate(seq)
        print 's :', seq

def RichGenbankParser(handle, info_excludes=None, moltype=None,
    skip_contigs=False):
    """Returns annotated sequences from GenBank formatted file.
    
    Arguments:
        - info_excludes: a series of fields to be excluded from the Info object
        - moltype: a MolType instance, such as PROTEIN, DNA. Default is ASCII.
        - skip_contigs: ignores records with no actual sequence data, typically
          a genomic contig."""
    info_excludes = info_excludes or []
    moltype = moltype or ASCII
    for rec in MinimalGenbankParser(handle):
        info = Info()
        # populate the Info object, excluding the sequence
        for label, value in rec.items():
            if label in info_excludes:
                continue
            info[label] = value
        
        if rec['mol_type'] == 'protein':  # which it doesn't for genbank
            moltype = PROTEIN
        elif rec['mol_type'] == 'DNA':
            moltype = DNA
        
        try:
            seq = moltype.makeSequence(rec['sequence'].upper(), Info=info,
                                        Name=rec['locus'])
        except KeyError:
            if not skip_contigs:
                if 'contig' in rec:
                    yield rec['locus'], rec['contig']
                elif 'WGS' in rec:
                    yield rec['locus'], rec['WGS']
                else:
                    yield rec['locus'], None
            continue
        
        for feature in rec['features']:
            spans = []
            reversed = None
            if feature['location'] == None or feature['type'] in ['source', \
                                                                'organism']:
                continue
            for location in feature['location']:
                (lo, hi) = (location.first() - 1, location.last())
                if location.Strand == -1:
                    (lo, hi) = (hi, lo)
                    assert reversed is not False
                    reversed = True
                else:
                    assert reversed is not True
                    reversed = False
                # ensure we don't put in a span that starts beyond the sequence
                if lo > len(seq):
                    continue
                # or that's longer than the sequence
                hi = [hi, len(seq)][hi > len(seq)]
                spans.append((lo, hi))
            if reversed:
                spans.reverse()
            
            for id_field in ['gene', 'note', 'product', 'clone']:
                if id_field in feature:
                    name = feature[id_field]
                    if not isinstance(name, basestring):
                        name = ' '.join(name)
                    break
            else:
                name = None
            seq.addAnnotation(Feature, feature['type'], name, spans)
        yield (rec['locus'], seq)
        
def parse(*args):
    return RichGenbankParser(*args).next()[1]

if __name__ == '__main__':  #demo if called from commandline
    from sys import argv
    rec = parse(open(argv[1], 'U'))
    print len(rec), rec.getName()
    for annot in rec.annotations:
        print annot