source: orange-bioinformatics/Orange/bioinformatics/obiKEGG.py @ 1625:cefeb35cbfc9

Revision 1625:cefeb35cbfc9, 41.8 KB checked in by mitar, 2 years ago (diff)

Moving files around.

Line 
1"""
2=======
3obiKEGG
4=======
5obiKEGG module is an interface to `Kyoto Encyclopedia of Genes and Genomes
6<http://www.genome.jp/kegg/>`_ that allows easy access to KEGG pathway
7and genes data.
8
9Example ::
10   
11    >>> genes = KEGGOrganism("human")
12    >>> genes
13"""
14try:
15    import Image, ImageDraw, ImageMath
16except:
17    pass
18
19import cStringIO
20import math
21import time
22import os, sys, tarfile
23import re
24
25import obiData
26import obiProb
27import orngServerFiles
28
29from cPickle import load, loads, dump
30from collections import defaultdict
31
32
33import xml.dom.minidom as minidom
34
35from functools import partial, wraps
36import obiTaxonomy
37import orngEnviron
38import urllib2
39import cPickle
40
41DEFAULT_DATABASE_PATH = orngServerFiles.localpath("KEGG")
42KEGG_FTP_PATH = "ftp://ftp.genome.jp/pub/kegg/"
43
44_join = lambda *args: os.path.join(DEFAULT_DATABASE_PATH, *args)
45
46def special_tags(domain, filename):
47    return dict([tuple(tag.split(":", 1)) for tag in orngServerFiles.info(domain, filename)["tags"] \
48                 if tag.startswith("#") and ":" in tag])
49def loads(domain, filename, version=None):
50    """ A function decorator for a function that loads data from server files.
51    Makes sure the file is downloaded and of the right version
52    """
53    def decorator(func):
54        @wraps(func)
55        def wrapper(*args, **kwargs):
56            try:
57                ## First try downloading the required server file 
58                file = filename if isinstance(filename, basestring) else filename(*args, **kwargs)
59                orngServerFiles.localpath_download(domain, file)
60                if version is not None and special_tags(domain, file).get("#version", str(version)) < str(version):
61                    orngServerFiles.update(domain, file) 
62            except Exception, ex:
63                print ex, args, kwargs
64            return func(*args, **kwargs)
65        return wrapper
66    return decorator
67           
68def deprecated(first_arg, text=None):
69    import warnings
70    def wrapper(*args, **kwargs):
71        warnings.warn_explicit(("Call to deprecated function %s. " % first_arg.__name__) + (text or ""), 
72                               category=DeprecationWarning,
73                               filename=first_arg.func_code.co_filename,
74                               lineno=first_arg.func_code.co_firstlineno + 1
75                      )
76#        warnings.warn(("Call to deprecated function %s." % first_arg.__name__) + (text or ""),
77#                               category=DeprecationWarning,;
78#                               stacklevel =  -1 #-2 if text is None else -3
79#                        )
80        return first_arg(*args, **kwargs)
81    if isinstance(first_arg, basestring): #deprecated called with a string as first argument
82        return partial(deprecated, text=first_arg)
83    else:
84        return wraps(first_arg)(wrapper)
85       
86def cached(func):
87    @wraps(func)
88    def wrapper(self, *args, **kwargs):
89        sig = args + tuple(sorted(kwargs.items()))
90        if sig in wrapper:
91            pass
92    return wrapper
93           
94def cached_property(func):
95    @property
96    @wraps(func)
97    def wrapper_func(self):
98        cache_name = "_" + func.__name__
99        if not hasattr(self, cache_name):
100            setattr(self, cache_name, func(self))
101        return getattr(self, cache_name)
102    return wrapper_func
103
104def cached_method(func, cache_name="_cached_method_cache", store=None):
105    def wrapper(self, *args, **kwargs):
106        sig = (func.__name__,) + args + tuple(sorted(kwargs.items()))
107        if not hasattr(self, cache_name):
108            setattr(self, cache_name, store() if store is not None else {})
109        if sig not in getattr(self, cache_name):
110            getattr(self, cache_name)[sig] = func(self, *args, **kwargs)
111        return getattr(self, cache_name)[sig]
112    return wrapper
113
114class pickle_cache(dict):
115    def __init__(self, filename, sync=True):
116        dict.__init__(self)
117        self.filename = filename
118        if os.path.exists(self.filename):
119            try:
120                self.update(cPickle.load(open(self.filename, "rb")))
121            except Exception, ex:
122                print self.filename, "loading failed", ex
123               
124    def __setitem__(self, key, value):
125        super(pickle_cache, self).__setitem__(key, value)
126        self.sync()
127       
128    def __delitem__(self, key):
129        super(pickle_cache, self).__delitem__(key)
130        self.sync()
131       
132    def sync(self):
133#        print "sync", self, self.filename
134        cPickle.dump(dict(self.items()), open(self.filename, "wb"), 2)
135       
136def persistent_cached_method(func):
137    return cached_method(func, cache_name="_persistent_cached_method_cache")
138
139def persistent_cached_class(cls, store=pickle_cache):
140    def cache_loader(self):
141        if not hasattr(self, "_cache_filename"):
142            setattr(self, "_cache_filename", getattr(self, "filename") + ".cache")
143#        print "creating cache", self._cache_filename
144        try:
145#            print os.stat(self.filename).st_mtime, os.stat(self._cache_filename).st_mtime
146            if os.stat(self.filename).st_mtime > os.stat(self._cache_filename).st_mtime:
147#                print "cache obsolete"
148                os.remove(self._cache_filename)
149        except Exception, ex:
150            pass
151        return store(self._cache_filename)
152    cls._persistent_cached_method_cache = cached_property(cache_loader)
153    return cls
154
155class cache_manager(object):
156    """ takes an fuction that returns a dict like object (cache)
157    """
158    def __init__(self, opener):
159        self.opener = opener
160        self.name = opener.__name__
161       
162    def __call__(self, *args, **kwargs):
163        return self.opener(self, *args, **kwargs)
164           
165    def __get__(self, instance, type_=None):
166        """ Return an cache store from instance """
167        ## TODO: check if instance is a type
168        if getattr(instance, self.name):
169            return getattr(instance, self.name)
170        return self.opener.__get__(instance, type)
171   
172    def cached_method(manager, wrapped):
173        def wrapper(self, *args, **kwargs):
174            sig = (wrapped.__name__,) + args + tuple(sorted(kwargs.items()))
175            cache = getattr(self, manager.__name__)
176            if sig not in cache:
177                cache[sig] = func(self, *args, **kwargs)
178            return cache[sig]
179
180def proxy_dict_decorator(cls, proxy_name):
181    for name in ["__len__", "__contains__", "__getitem__", "__iter__", "get", "has_key", "items", 
182                 "iteritems", "iterkeys", "itervalues", "keys", "values", "update"]:
183        def proxy(func):
184            def f(self, *args, **kwargs):
185                return func(getattr(self, proxy_name), *args, **kwargs)
186            return f
187        setattr(cls, name, proxy(getattr(dict, name)))
188    return cls
189
190def _create_path_for_file(target): 
191    try:
192        os.makedirs(os.path.dirname(target))
193    except OSError:
194        pass
195
196def downloader(func):
197    def wrapper(*args, **kwargs):
198        import shutil
199        def download(src, dst):
200            if isinstance(dst, basestring):
201                _create_path_for_file(dst)
202                dst = open(dst, "wb")
203            shutil.copyfileobj(src, dst)
204        jobs = func(*args, **kwargs)
205        for src, dst in [jobs] if type(jobs) == tuple else jobs:
206            download(src, dst)
207    return wrapper
208
209from orngMisc import with_gc_disabled
210   
211           
212SLINE, MLINE, LINELIST = range(1, 4)
213
214class _field_mutators(object):
215    @classmethod
216    def deindent(cls, text, field):
217        pad = len(text) - len(text[len(field):].lstrip())
218        return "\n".join([line[pad:] for line in text.splitlines()])
219   
220    @classmethod
221    def to_list(cls, text, field):
222        return [line for line in text.splitlines() if line.strip()]
223   
224    @classmethod
225    def to_id_list(cls, text, field):
226        return [tuple(line.split(": ", 1)) for line in cls.to_list(cls.deindent(text, field), field)]
227   
228    @classmethod
229    def to_dict(cls, text, field):
230        text = cls.deindent(text, field).replace("\n ", " ")
231        lines = [t.split(": ",1) for t in text.split("\n")]
232        return dict([(key, [t.strip() for t in value.split(" ") if t]) for key, value in lines])
233   
234    @classmethod
235    def to_ids(cls, text, field):
236        lines = [line.split(" ", 2)[:2] for line in cls.deindent(text, field).split("\n") if ":" in line.split(" ", 1)[0]]
237        return [db.lower() + id for db, id in lines] 
238#        ids = reduce(lambda col, line:1, text.split("\n"), [])
239        return ids
240    default = deindent
241#    entry = classmethod(lambda cls, text, field: cls.deindent(text, field).split()[0])
242    orthology = to_id_list
243    pathway = classmethod(lambda cls, text, field: [t[1].split()[0] for t in cls.to_id_list(text, field)])
244    codon_usage = classmethod(lambda cls, text, field: text.replace(field, " " * len(field)))
245    structure = to_id_list
246    motif = to_dict
247    pathway = to_ids
248    orthology = to_ids
249   
250    @classmethod
251    def genes(cls, text, field):
252        genes = [line.split(": ", 1) for line in cls.deindent(text, field).replace("\n    ", "").splitlines()]
253        return dict([(org.lower(), [gene.split("(")[0] for gene in genes.split(" ")]) for org, genes in genes])
254   
255    @classmethod
256    def dblinks(cls, text, field):
257        links = [line.split(": ", 1) for line in cls.deindent(text, field).replace("\n   ", "").splitlines()]
258        return dict([(db.lower(), [name for name in links.split(" ")]) for db, links in links])
259   
260    dblinks = to_dict
261
262def entry_decorate(cls):
263    reserved_map = {"class": "class_"}
264    for field in cls.FIELDS:
265        mutator = getattr(_field_mutators, field.lower(), _field_mutators.default)
266        def construct(field, mutator):
267            def get(self):
268                text = self.get(field)
269                return mutator(text, field) if text else None
270            return get
271        if not hasattr(cls, field.lower()):
272            setattr(cls, reserved_map.get(field.lower(), field.lower()), property(construct(field, mutator)))
273        else:
274            import warnings
275            warnings.warn(str(cls)+ "already contains " + field.lower())
276    return cls
277
278def borg_class(cls):
279    borg__init__ = cls.__init__
280    cls._borg_instance = None
281    def __init__(self, *args, **kwargs):
282        if not cls._borg_instance:
283            borg__init__(self, *args, **kwargs)
284            cls._borg_instance = args + tuple(sorted(kwargs.items())), self
285        elif cls._borg_instance[0] == args + tuple(sorted(kwargs.items())):
286            self.__dict__ = cls._borg_instance[1].__dict__
287        else:
288            print "Warning. Atempting to create an individual instance of '%s'.\nWe are borg. You will be assimilated. Resistance is futile." % cls.__name__
289            self.__dict__ = cls._borg_instance[1].__dict__
290    cls.__init__ = __init__
291    return cls
292
293def defaultpath(func):
294#    @wraps
295    def wrapper(*args, **kwargs):
296        if "path" not in kwargs:
297            kwargs["path"] = DEFAULT_DATABASE_PATH
298        return func(*args, **kwargs)
299    return wrapper
300             
301       
302class KEGGDBEntry(object):
303    FIELDS = []
304    MULTIPLE_FIELDS = []
305    def __init__(self, entrytext, index=None):
306        self.entrytext = entrytext
307        if not self.FIELDS:
308            self.FIELDS = set([line.split()[0] for line in entrytext.split("\n") if line.strip() and not line.startswith(" ")])
309        self._index(index)
310       
311    def _indices(self, field):
312        text = self.entrytext
313        fieldlen = len(field)
314        return reduce(lambda indices,s: indices.append(text.find(field, indices[-1] + fieldlen)) or indices, range(text.count(field)), [-fieldlen])[1:]
315         
316    def _index(self, index=None):
317        if index is not None:
318            index, self.partitions = index
319            self.index = dict(zip(self.FIELDS, index))
320            return 
321        index = dict([(field, self._indices(field)) for field in self.FIELDS])
322        sorted_indices = sorted(reduce(set.union, index.values(), set())) + [-1]
323        self.partitions = sorted_indices #[(sorted_indices[i], sorted_indices[i+1]) for i in range(len(sorted_indices) - 1)]
324        self.index = dict([(key, [sorted_indices.index(ind) for ind in indices]) for key, indices in index.items()])
325       
326    def section_iter(self, field):
327        for ind in self.index[field]:
328            yield self.entrytext[self.partitions[ind]: self.partitions[ind+1]]
329       
330    def get(self, key, default=None):
331        if key not in self.FIELDS:
332            return default
333        sections = [section for section in self.section_iter(key)]
334        if key in self.MULTIPLE_FIELDS:
335            return sections
336        elif sections:
337            return sections[0]
338        else: 
339            return default
340       
341    def entry_key(self):
342        return self.entry.split()[0]
343       
344class KEGGDataBase(object):
345    Entry = KEGGDBEntry
346    VERSION = "v1.0"
347    def __init__(self, file=None):
348        self.load(file)
349       
350    @with_gc_disabled
351    def load(self, file=None):
352        if file is None:
353            file = self.FILENAME % {"path":DEFAULT_DATABASE_PATH}
354        self.filename = file
355        data = open(file, "rb").read().split("///\n")
356        file_timestamp = os.stat(self.filename).st_mtime
357#        print file_timestamp
358        build_index = False
359#        if os.path.exists(self.filename + ".index"):
360        try:
361            version, timestamp, index = cPickle.load(open(self.filename + ".index", "rb"))
362            assert(version == self.VERSION and timestamp == file_timestamp)
363        except Exception, ex:
364#            print "error", ex
365            build_index = True
366            index = [None] * len(data)
367        self.entrys = [self.Entry(entry, ind) for entry, ind in zip(data, index) if entry.strip()]
368        if build_index:
369#            self.entrys = [self.Entry(entry) for entry in data if entry.strip()]
370            index = [([e.index[field] for field in e.FIELDS], e.partitions) for e in self.entrys]
371            cPickle.dump((self.VERSION, file_timestamp, index), open(self.filename + ".index", "wb"), cPickle.HIGHEST_PROTOCOL)
372       
373        self.entry_dict = dict([(entry.entry_key(), entry) for entry in self.entrys])
374       
375proxy_dict_decorator(KEGGDataBase, "entry_dict")
376       
377class KEGGGeneEntry(KEGGDBEntry):
378    FIELDS = ["ENTRY", "NAME", "DEFINITION","ORTHOLOGY", "PATHWAY", "CLASS", "POSITION",
379              "MOTIF", "DBLINKS", "STRUCTURE", "CODON_USAGE", "AASEQ", "NTSEQ"]
380   
381    def aliases(self):
382        return [self.entry_key()] + (self.name.split(",") if self.name else []) + [link[1][0] for link in self.dblinks.items() if self.dblinks]
383   
384    @property
385    def alt_names(self):
386        return self.aliases()
387       
388entry_decorate(KEGGGeneEntry)
389   
390class KEGGGenes(KEGGDataBase):
391    Entry = KEGGGeneEntry
392    FILENAME = "genes/organisms/%(org_code)s/%(org_name)s.ent"
393   
394    @loads("KEGG", lambda self, org: "kegg_genes_%s.tar.gz" % org)
395    def load(self, org):
396        super(KEGGGenes, self).load(_join(self.filename(org)))
397        self.entry_dict = dict([(org + ":" + key, value) for key, value in self.entry_dict.items()])
398        self.org_code = org
399       
400    def gene_aliases(self):
401        aliases = {}
402        for entry in self.entrys:
403            aliases.update(dict.fromkeys(entry.aliases(), self.org_code + ":" + entry.entry_key()))
404        return aliases
405           
406    @classmethod
407    def filename(cls, org):
408        return cls.FILENAME % {"org_code":org, "org_name":KEGGGenome()[org].name.split(",")[1].strip()}
409     
410    @classmethod
411    @downloader
412    def download(cls, org):
413#        local_dir = orngServerFiles.localpath("KEGG") if local_dir is None else local_dir
414        filename = cls.filename(org)
415        return (urllib2.urlopen(KEGG_FTP_PATH + filename), 
416                _join(filename))
417   
418class KEGGGenomeEntry(KEGGDBEntry):
419    FIELDS = ["ENTRY", "NAME", "DEFINITION", "ANNOTATION", "TAXONOMY", "DATA_SOURCE",
420              "ORIGINAL_DB", "CHROMOSOME", "STATISTICS", "REFERENCE"]
421    MULTIPLE_FIELDS = ["REFERENCE"]
422   
423    def org_code(self):
424        if self.name is not None:
425            return self.name.split(",")[0]
426        else:
427            return self.entry.split()[0]
428       
429    def entry_key(self):
430        return self.org_code()
431   
432entry_decorate(KEGGGenomeEntry)
433       
434class KEGGGenome(KEGGDataBase):
435    Entry = KEGGGenomeEntry
436    VERSION = "v1.1"
437    TAXID_MAP = {"4932" : "559292"}
438    def name(self):
439        return super(KEGGGeneome, self).name()
440   
441    @loads("KEGG", "kegg_genome.tar.gz", version=VERSION)
442    def load(self, file=None):
443        filename = _join("genes", "genome")
444        super(KEGGGenome, self).load(filename)
445   
446    @classmethod
447    def common_organisms(cls):
448#        genome = KEGGGenome()
449#        id_map = {"562":"511145", "2104":"272634", "5833":"36329", "4896":"284812", "11103":None, "4754":None, "4577":None}
450#        return [genome.get(id_map.get(taxid, taxid)).entry_key() for taxid in obiTaxonomy.common_taxids() if id_map.get(taxid, taxid) is not None]
451        return ['ath', 'bta', 'cel', 'cre', 'dre', 'ddi', 'dme', 'eco', 'hsa', 'mmu', 'mpn', 'osa',
452                'pfa', 'rno', 'sce', 'spo', 'zma', 'xla']
453       
454    @classmethod
455    def essential_organisms(cls):
456        genome = KEGGGenome()
457       
458        return [genome.get(cls.TAXID_MAP.get(taxid, taxid)).entry_key() for taxid in obiTaxonomy.essential_taxids()]
459           
460    def get(self, key, *args):
461        if key in self.TAXID_MAP:
462            key = self.TAXID_MAP[key]
463        if key not in self:
464            keys = self.search(key)
465            keys = [k for k in keys if key in self[k].name or key in self[k].taxonomy]
466            key = keys.pop(0) if keys else key
467        return super(KEGGGenome, self).get(key, *args)
468   
469    def search(self, string, relevance=False):
470        if string in self.TAXID_MAP:
471            string = self.TAXID_MAP[string]
472           
473        def match(entry, string):
474            rel = 0
475            if string in entry.entrytext:
476                weight_f = [(entry.definition or "", 2), (entry.name or "", 4), (entry.taxonomy or "", 8), (entry.entry_key(), 16)]
477                rel += 1 + sum([w for text, w in weight_f if string in text]) + (16 if entry.entry_key() in self.common_organisms() else 0)
478            return rel
479        matched = sorted(zip([match(entry, string) for entry in self.entrys], self.entrys), reverse=True)
480        return [(entry.entry_key(), rel) if relevance else entry.entry_key() for rel, entry in matched if rel]
481           
482    @classmethod
483    @downloader
484    def download(cls, file=None):
485        return (urllib2.urlopen(KEGG_FTP_PATH + "/genes/genome"), 
486                _join("genes", "genome") if file is None else file)
487       
488borg_class(KEGGGenome)
489       
490class KEGGCompoundEntry(KEGGDBEntry):
491    FIELDS = ["ENTRY", "NAME", "FORMULA", "MASS", "REMARK", "REACTION", "PATHWAY",
492              "ENZYME", "DBLINKS", "ATOM", "BOND"]
493    def entry_key(self):
494        return "cpd:" + super(KEGGCompoundEntry, self).entry_key()
495   
496entry_decorate(KEGGCompoundEntry)
497   
498class KEGGCompounds(KEGGDataBase):
499    Entry = KEGGCompoundEntry
500    FILENAME = "%(path)s/ligand/compound/compound"
501   
502    @loads("KEGG", "kegg_ligand.tar.gz")
503    def load(self, file=None):
504        super(KEGGCompounds, self).load(file)
505       
506    @classmethod
507    @downloader
508    def download(cls, file=None):
509        return (urllib2.urlopen(KEGG_FTP_PATH + "/ligand/compound/compound"),
510                _join("ligand", "compound", "compound") if file is None else file)
511   
512borg_class(KEGGCompounds)
513   
514class KEGGEnzymeEntry(KEGGDBEntry):
515    FIELDS = ["ENTRY", "NAME", "CLASS", "SYSNAME", "REACTION", "ALL_REAC", "SUBSTRATE",
516              "PRODUCT", "COFACTOR", "COMMENT", "REFERENCE", "PATHWAY", "ORTHOLOGY", 
517              "GENES", "STRUCTURES", "DBLINKS"]
518    MULTIPLE_FIELDS = ["REFERENCE"]
519    def entry_key(self):
520        return "EC:" + self.entry.split()[1]
521   
522entry_decorate(KEGGEnzymeEntry)
523
524class KEGGEnzymes(KEGGDataBase):
525    Entry = KEGGEnzymeEntry
526    FILENAME = "%(path)s/ligand/enzyme/enzyme"
527
528    @loads("KEGG", "kegg_ligand.tar.gz")
529    def load(self, file=None):
530        super(KEGGEnzymes, self).load(file)
531   
532    @classmethod
533    @downloader
534    def download(cls, file=None):
535        return (urllib2.urlopen(KEGG_FTP_PATH + "ligand/enzyme/enzyme"),
536                _join("ligand", "enzyme", "enzyme"))
537   
538borg_class(KEGGEnzymes)
539
540class KEGGReactionEntry(KEGGDBEntry):
541    FIELDS = ["ENTRY", "NAME", "DEFINITION", "EQUATION", "COMMENT", "RPAIR", "PATHWAY", 
542              "ENZYME", "ORTHOLOGY"]
543    def entry_key(self):
544        return "rn:" + super(KEGGReactionEntry, self).entry_key()
545   
546entry_decorate(KEGGReactionEntry)
547
548class KEGGReactions(KEGGDataBase):
549    Entry = KEGGReactionEntry
550    FILENAME = "%(path)s/ligand/reaction/reaction"
551   
552    @loads("KEGG", "kegg_ligand.tar.gz")
553    def load(self, file=None):
554        super(KEGGReactions, self).load(file)
555       
556    @classmethod
557    @downloader
558    def download(cls, file=None):
559        return (urllib2.urlopen(KEGG_FTP_PATH + "ligand/reaction/reaction"),
560                _join("ligand", "reaction", "reaction"))
561       
562borg_class(KEGGReactions)
563
564class KEGGKOEntry(KEGGDBEntry):
565    FIELDS = ["ENTRY", "NAME", "DEFINITION", "CLASS", "DBLINKS", "GENES"]
566   
567    def entry_key(self):
568        return "ko:" + super(KEGGKOEntry, self).entry_key()
569   
570entry_decorate(KEGGKOEntry)
571   
572class KEGGKO(KEGGDataBase):
573    Entry = KEGGKOEntry
574    FILENAME = "%(path)s/genes/ko"
575   
576    @loads("KEGG", "kegg_orthology.tar.gz")
577    def load(self, file=None):
578        super(KEGGKO, self).load(file)
579       
580    @classmethod
581    @downloader
582    def download(cls, file=None):
583        return (urllib2.urlopen(KEGG_FTP_PATH + "genes/ko"),
584                _join("genes", "ko"))
585       
586borg_class(KEGGKO)
587
588class KEGGBriteEntry(object):
589    _search_re = {"ids": re.compile('(?P<ids>\[.*:.*\])'),
590                  "title": re.compile(r'(<[Bb]>)?(?P<title>\b[a-zA-Z0-9_/\s,;:.+=\-\[\]{}\(\)]+?)(?(1)</[Bb]>)$'),
591                  "links": re.compile('(?P<links><a href=".+?">.*?</a>)')}
592    def __init__(self, line, entrys=None):
593        self.entrys = [] #entrys if entrys is not None else []
594        self.line = line[1:].strip()
595        for name, re in self._search_re.items():
596            search = re.search(self.line)
597            setattr(self, name, search.group(name) if search else None)
598#        self.identifiers = self.groups.get("ids")
599#        self.title = self.groups. # TODO: parse line (html, kegg identifiers ...)
600
601    def __iter__(self):
602        return iter(self.entrys)
603
604class KEGGBrite(KEGGBriteEntry):
605    VERSION = "v1.0"
606    def __init__(self, brite_id):
607        super(KEGGBrite, self).__init__("")
608        self.load(brite_id)
609       
610    @loads("KEGG", "kegg_brite.tar.gz")
611    def load(self, brite_id):
612        file = self.filename(brite_id)
613        lines = open(file, "rb").read().split("\n!\n")[1].splitlines()
614        def collect(lines, depth, collection):
615            while lines:
616                line = lines[0]
617                if line.startswith("#"):
618                    lines.pop(0)
619                elif line.startswith(depth) and len(line.strip()) > 1:
620                    collection.append(KEGGBriteEntry(lines.pop(0))) 
621                elif line[0] > depth:
622                    collect(lines, line[0], collection[-1].entrys)
623                elif line[0] < depth:
624                    return
625                else:
626                    lines.pop(0)
627                       
628        collect([line for line in lines if not line.startswith("#") and len(line) > 1], "A", self.entrys)
629   
630    @classmethod
631    @defaultpath
632    def filename(cls, brite_id, path=DEFAULT_DATABASE_PATH):
633        if path is None:
634            path = "%(path)s"
635        if brite_id.startswith("br"):
636            return ("%(path)s/brite/br/" + brite_id + ".keg") % dict(path=path)
637        elif brite_id.startswith("ko"):
638            return ("%(path)s/brite/ko/" + brite_id + ".keg") % dict(path=path)
639        else:
640            org = brite[:-5]
641            return ("%(path)s/brite/organisms/" + org + "/" + brite_id + ".keg") % dict(path=path)
642       
643    @classmethod
644    @downloader
645    def download(cls, brite_id):
646        filename = cls.filename(brite_id, path="")
647        return (urllib2.urlopen(cls.filename(brite_id, path=KEGG_FTP_PATH.rstrip("/"))),
648                cls.filename(brite_id))
649
650class KEGGOrganism(object):
651    VERSION = "v2.0"
652    DOMAIN = "KEGG"
653    def __init__(self, org, genematcher=None):
654        self.load(org)
655        self.genematcher = genematcher
656       
657#    @loads("KEGG", lambda self, org: "kegg_genes_%s.tar.gz" % self.organism_name_search(org))
658    def load(self, org):
659        org = self.organism_name_search(org)
660        self.org_code = org
661   
662    @property
663#    @deprecated("Use org_code instead")
664    def org(self):
665        return self.org_code
666   
667    @cached_property
668    def genes(self):
669        return KEGGGenes(self.org_code)
670       
671    @cached_property
672    def gene_aliases(self):
673        return self.genes.gene_aliases()
674   
675    def pathways(self, with_ids=None):
676        return [pathway for pathway, values in KEGGPathway.list(self.org_code).items() if all([id in values for id in (with_ids or [])])]
677   
678#    @deprecated("Use KEGGOrganism.pathways instead")
679    def list_pathways(self):
680        return self.pathways()
681   
682#    @deprecated
683    def get_linked_pathways(self, pathway_id):
684        return KEGGPathway(pathway_id).genes()
685   
686    def enzymes(self, genes=None):
687        enzymes = KEGGEnzymes()
688        return [enzyme.entry_key() for enzyme in enzymes.itervalues() if enzyme.genes and self.org_code in enzyme.genes]
689   
690    def get_enriched_pathways(self, genes, reference=None, prob=obiProb.Binomial(), callback=None):
691        """Return a dictionary with enriched pathways ids as keys and (list_of_genes, p_value, num_of_reference_genes) tuples as items."""
692        allPathways = defaultdict(lambda :[[], 1.0, []])
693        import orngMisc
694        milestones = orngMisc.progressBarMilestones(len(genes), 100)
695        for i, gene in enumerate(genes):
696            pathways = self.pathways([gene])
697            for pathway in pathways:
698                allPathways[pathway][0].append(gene)
699            if callback and i in milestones:
700                callback(i*100.0/len(genes))
701        reference = set(reference if reference is not None else self.genes.keys())
702        for p_id, entry in allPathways.items():
703            entry[2].extend(reference.intersection(KEGGPathway(p_id).genes()))
704            entry[1] = prob.p_value(len(entry[0]), len(reference), len(entry[2]), len(genes))
705        return dict([(pid, (genes, p, len(ref))) for pid, (genes, p, ref) in allPathways.items()])
706   
707    def get_genes_by_enzyme(self, enzyme):
708        enzyme = KEGGEnzymes()[enzyme]
709        return enzyme.genes.get(self.org_code, []) if enzyme.genes else []
710   
711    def get_genes_by_pathway(self, pathway_id):
712        return KEGGPathway(pathway_id).genes()
713   
714    def get_enzymes_by_pathway(self, pathway_id):
715        return KEGGPathway(pathway_id).enzymes()
716   
717    def get_compounds_by_pathway(self, pathway_id):
718        return KEGGPathway(pathway_id).compounds()
719   
720    def get_pathways_by_genes(self, gene_ids):
721        gene_ids = set(gene_ids)
722        pathways = [self.genes[id].pathway for id in gene_ids if self.genes[id].pathway]
723        pathways = reduce(set.union, pathways, set())
724        return [id for id in pathways if gene_ids.issubset(KEGGPathway(id).genes())] 
725   
726    def get_pathways_by_enzymes(self, enzyme_ids):
727        enzyme_ids = set(enzyme_ids)
728        pathways = [KEGGEnzymes()[id].pathway for id in enzyme_ids]
729        pathwats = reduce(set.union, pathways, set())
730        return [id for id in pathways if enzyme_ids.issubset(KEGGPathway(id).enzymes())]
731   
732    def get_pathways_by_compounds(self, compound_ids):
733        compound_ids = set(compound_ids)
734        pathways = [KEGGCompounds()[id].pathway for id in compound_ids]
735        pathwats = reduce(set.union, pathways, set())
736        return [id for id in pathways if compound_ids.issubset(KEGGPathway(id).compounds())]
737   
738    def get_enzymes_by_compound(self, compound_id):
739        return KEGGCompound()[compound_id].enzyme
740   
741    def get_enzymes_by_gene(self, gene_id):
742        return self.genes[gene_id].enzymes
743   
744    def get_compounds_by_enzyme(self, enzyme_id):
745        return self._enzymes_to_compounds.get(enzyme_id)
746   
747    def get_unique_gene_ids(self, genes, caseSensitive=True):
748        """Return a tuple with three elements. The first is a dictionary mapping from unique gene
749        ids to gene names in genes, the second is a list of conflicting gene names and the third is a list
750        of unknown genes.
751        """
752        unique, conflicting, unknown = {}, [], []
753        for gene in genes:
754            names = self.genematcher.match(gene)
755            if len(names) == 1:
756                unique[names[0]] = gene
757            elif len(names) == 0:
758                unknown.append(gene)
759            else:
760                conflicting.append(gene)
761        return unique, conflicting, unknown
762   
763    @property
764    @cached_method
765    def _enzymes_to_compounds(self):
766        dd = {}
767        for val in KEGGCompounds().values():
768            dd.update(dict.fromkeys(val.enzymes, val.entry_key()))
769        return dd
770   
771    def _set_genematcher(self, genematcher):
772        setattr(self, "_genematcher", genematcher)
773       
774    def _get_genematcher(self):
775        if getattr(self, "_genematcher", None) == None:
776            import obiGene
777            if self.org_code == "ddi":
778                self._genematcher = obiGene.matcher([obiGene.GMKEGG(self.org_code), obiGene.GMDicty(),
779                                                     [obiGene.GMKEGG(self.org_code), obiGene.GMDicty()]])
780            else:
781                self._genematcher = obiGene.matcher([obiGene.GMKEGG(self.org_code)])
782            self._genematcher.set_targets(self.genes.keys())
783        return self._genematcher
784    genematcher = property(_get_genematcher, _set_genematcher)
785   
786    def get_genes(self):
787        return self.genes
788   
789    @classmethod
790    def organism_name_search(cls, name):
791        genome = KEGGGenome()
792        if name not in genome:
793            ids = genome.search(name)
794            if not ids:
795                import obiTaxonomy
796                ids = obiTaxonomy.search(name)
797                ids = [id for id in ids if genome.search(id)]
798            name = ids.pop(0) if ids else name
799           
800        try:
801            return KEGGGenome()[name].entry_key()
802        except KeyError:
803            raise ValueError("Organims with name='%s' not found in KEGG." % name)
804       
805    @classmethod
806    def organism_version(cls, name):
807        name = cls.organism_name_search(name)
808        orngServerFiles.localpath_download("KEGG", "kegg_genes_%s.tar.gz" % name)
809        return orngServerFiles.info("KEGG", "kegg_genes_%s.tar.gz" % name)["datetime"]
810             
811class KEGGPathway(object):
812    PNG_FILENAME = "%(path)s/pathway/organisms/%(org)s/%(org)s%(map_id)s.png"
813    KGML_FILENAME = "%(path)s/xml/kgml/metabolic/organisms/%(org)s/%(org)s%(map_id)s.xml"
814    VERSION = "v2.0"
815   
816    class entry(object):
817        def __init__(self, dom_element):
818            self.__dict__.update(dom_element.attributes.items())
819            self.graphics = ()
820            self.components = []
821            self.graphics = dict(dom_element.getElementsByTagName("graphics")[0].attributes.items())
822            self.components = [node.getAttribute("id") for node in dom_element.getElementsByTagName("component")]
823    class reaction(object):
824        def __init__(self, dom_element):
825            self.__dict__.update(dom_element.attributes.items())
826            self.substrates = [node.getAttribute("name") for node in dom_element.getElementsByTagName("substrate")]
827            self.products = [node.getAttribute("name") for node in dom_element.getElementsByTagName("product")]
828           
829    class relation(object):
830        def __init__(self, dom_element):
831            self.__dict__.update(dom_element.attributes.items())
832            self.subtypes = [node.attributes.items() for node in dom_element.getElementsByTagName("subtype")]
833       
834    def __init__(self, file):
835        if not os.path.exists(file):
836            path, org, map_id = self.split_pathway_id(file)
837#            file = self.KGML_FILENAME % dict(path=DEFAULT_DATABASE_PATH, org=org, map_id=id)
838            file = self.filename_kgml(org, map_id)
839        if not os.path.exists(file) and os.path.exists(file.replace("metabolic", "non-metabolic")):
840            file = file.replace("metabolic", "non-metabolic")
841        self.filename = file
842#        self.load(file)
843       
844    @persistent_cached_method
845    def pathway_attributes(self):
846        return dict(self.pathway_dom().attributes.items())
847   
848    @property
849    def name(self):
850        return self.pathway_attributes().get("name")
851   
852    @property
853    def org(self):
854        return self.pathway_attributes().get("org")
855   
856    @property
857    def number(self):
858        return self.pathway_attributes().get("number")
859   
860    @property   
861    def title(self):
862        return self.pathway_attributes().get("title")
863   
864    @property
865    def image(self):
866        return self.pathway_attributes().get("image")
867   
868    @property
869    def link(self):
870        return self.pathway_attributes().get("link")
871   
872    @cached_method
873    @loads("KEGG", lambda self: "kegg_pathways_%s.tar.gz" % self.filename.split("/")[-1][:-9])
874    def pathway_dom(self):
875        return minidom.parse(self.filename).getElementsByTagName("pathway")[0]
876   
877    @cached_method
878    def entrys(self):
879        return [self.entry(e) for e in self.pathway_dom().getElementsByTagName("entry")]
880   
881    @cached_method
882    def reactions(self):
883        return [self.reaction(e) for e in self.pathway_dom().getElementsByTagName("reaction")]
884   
885    @cached_method
886    def relations(self):
887        return [self.relation(e) for e in self.pathway_dom().getElementsByTagName("relation")]
888   
889    @classmethod
890    def split_pathway_id(cls, id):
891        path, id = id.split(":") if ":" in id else ("path", id)
892        org, id = id[:-5], id[-5:]
893        return path, org, id 
894       
895    def __iter__(self):
896        """ Iterate over all elements in the pathway
897        """
898        return iter(self.all_elements())
899   
900    def __contains__(self, element):
901        """ Retrurn true if element in the pathway
902        """
903        return element in self.all_elements()
904   
905    def __getitem__(self, key):
906        return
907   
908    @cached_method
909    def all_elements(self):
910        """ Return all elements
911        """
912        return reduce(list.__add__, [self.genes(), self.compounds(), self.enzmes(), self.reactions()], [])
913   
914    def _get_entrys_by_type(self, type):
915        return reduce(set.union, [entry.name.split() for entry in self.entrys() if entry.type == type], set())
916   
917    @persistent_cached_method
918    def genes(self):
919        """ Return all genes on the pathway
920        """
921        return self._get_entrys_by_type("gene")
922   
923    @persistent_cached_method
924    def compounds(self):
925        """ Return all compounds on the pathway
926        """
927        return self._get_entrys_by_type("compound")
928   
929    @persistent_cached_method
930    def enzymes(self):
931        """ Return all enzymes on the pathway
932        """
933        return self._get_entrys_by_type("enzyme")
934   
935    @persistent_cached_method
936    def orthologs(self):
937        """ Return all orthologs on the pathway
938        """
939        return self._get_entrys_by_type("ortholog")
940   
941    @persistent_cached_method
942    def maps(self):
943        """ Return all linked maps on the pathway
944        """
945        return self._get_entrys_by_type("map")
946   
947    @persistent_cached_method
948    def groups(self):
949        """ Return all groups on the pathway
950        """
951        return self._get_entrys_by_type("ortholog")
952
953    def get_image(self):
954        """ Return an image of the pathway
955        """
956        return self.filename_png(self.org, self.number) % dict(path=DEFAULT_DATABASE_PATH)
957   
958    @persistent_cached_method
959    def get_bounding_box_dict(self):
960        return dict([(element.id, element.graphics) for element in self.entrys() if element.graphics])
961   
962    @persistent_cached_method
963    def graphics(self, item):
964        return [entry.graphics for entry in self.entrys() if item in entry.name and entry.graphics]
965       
966    @classmethod
967    @partial(cached_method, cache_name="_cls_cached_method_cache_")
968    @loads("KEGG", lambda cls, org: "kegg_pathways_%s.tar.gz" % org)
969    def list(cls, org):
970        file = cls.directory_png(org) + org + ".list"
971        data = [line.split() for line in open(file, "rb").read().splitlines()]
972        return reduce(lambda dict, line: dict[line[0]].update(line[1:]) or dict, data, defaultdict(set))
973   
974    @classmethod
975    @defaultpath
976    def filename_kgml(cls, org, map_id, path=DEFAULT_DATABASE_PATH):
977        path = "%(path)s" if path is None else path
978        if org in ["map", "ec"]:
979            return "%(path)s/xml/kgml/metabolic/ec/ec%(map_id)s.xml" % dict(org=org, map_id=map_id, path=path)
980        elif org == "ko":
981            return "%(path)s/xml/kgml/metabolic/ko/ko%(map_id)s.xml" % dict(org=org, map_id=map_id, path=path)
982        else:
983            return "%(path)s/xml/kgml/metabolic/organisms/%(org)s/%(org)s%(map_id)s.xml" % dict(org=org, map_id=map_id, path=path)
984       
985    @classmethod
986    @defaultpath
987    def filename_png(cls, org, map_id, path=DEFAULT_DATABASE_PATH):
988        path = "%(path)s" if path is None else path
989        if org in ["map", "ec", "ko", "rn"]:
990            return "%(path)s/pathway/%(org)s/%(org)s%(map_id)s.png" % dict(org=org, map_id=map_id, path=path)
991        else:
992            return "%(path)s/pathway/organisms/%(org)s/%(org)s%(map_id)s.png" % dict(org=org, map_id=map_id, path=path)
993       
994    @classmethod
995    @defaultpath
996    def directory_png(cls, org, path=DEFAULT_DATABASE_PATH):
997        return cls.filename_png(org, "", path=path).rsplit("/", 1)[0] + "/" 
998   
999    @classmethod
1000    @defaultpath
1001    def directory_kgml(cls, org, path=DEFAULT_DATABASE_PATH):
1002        return cls.filename_kgml(org, "", path=path).rsplit("/", 1)[0] + "/"
1003       
1004    @classmethod
1005    def pathways(cls, org):
1006        file = cls.directory_png(org) + org + ".list"
1007        pathways = [line.split()[0] for line in open(file, "rb").read().splitlines() if line.strip()]
1008        return sorted(set(pathways))
1009       
1010    @classmethod
1011    def download(cls, pathway, target=None):
1012        """ Download the pathway (xml and png files) to target directory
1013        """
1014        org = pathway[:3]
1015        import urllib2
1016        xml_data = urllib2.urlopen("ftp://ftp.genome.jp/pub/kegg/xml/kgml/metabolic/organisms/%s/%s.xml" % (org, pathway)).read()
1017        open(os.path.join(target, pathway + ".xml"), "wb").write(xml_data)
1018        png_data = urllib2.urlopen("ftp://ftp.genome.jp/pub/kegg/pathway/organisms/%s/%s.png" % (org, pathway)).read()
1019        open(os.path.join(target, pathway + ".png"), "wb").write(png_data)
1020       
1021    @classmethod
1022    @downloader
1023    def download_pathways(cls, org):
1024        png_path = cls.directory_png(org, path=None).replace("%(path)s/", "")
1025        xml_path = cls.directory_kgml(org, path=None).replace("%(path)s/", "")
1026        data = urllib2.urlopen("ftp://ftp.genome.jp/pub/kegg/" + png_path + org + ".list").read()
1027        pathways = sorted(set([line.split()[0][5:] for line in data.splitlines() if line]))
1028        from obiData import FtpDownloader
1029        ftp = FtpDownloader("ftp.genome.jp", _join(), "pub/kegg/", numOfThreads=15)
1030        ftp.massRetrieve([png_path + pathway + ".png" for pathway in pathways])
1031           
1032        if org != "map" :
1033            pathways = ftp.listdir(xml_path)
1034            pathways = [p for p in pathways if p.endswith(".xml")]
1035            ftp.massRetrieve([xml_path + pathway for pathway in pathways])
1036            if org not in ["ec", "rn"]: # Add non-metabolic pathways to the download list
1037                xml_path = xml_path.replace("metabolic", "non-metabolic")
1038                pathways = ftp.listdir(xml_path)
1039                pathways = [p for p in pathways if p.endswith(".xml")] 
1040                ftp.massRetrieve([xml_path + pathway for pathway in pathways])
1041        ftp.massRetrieve([png_path + org + ".list"])
1042        return []
1043   
1044persistent_cached_class(KEGGPathway)
1045       
1046def organism_name_search(name):
1047    return KEGGOrganism.organism_name_search(name)
1048
1049def pathways(org):
1050    return KEGGPathway.list(org)
1051
1052def organisms():
1053    return KEGGOrganism.organisms()
1054
1055def to_taxid(name):
1056    genome = KEGGGenome()
1057    names = genome.search(name)
1058    if genome[names[0]].taxonomy:
1059        return re.findall(r"TAX:\s*(\d+)", genome[names[0]].taxonomy)[0]
1060    else:
1061        return None #Should this raise an error?
1062
1063def from_taxid(taxid):
1064    return KEGGGenome().search(taxid)[0]
1065   
1066def test():
1067    p = KEGGPathway("sce00010.xml")
1068    print p.genes
1069    print p.reactions
1070    print p.compounds
1071    print p.image
1072    g = KEGGGenome()
1073    org = KEGGOrganism("Homo sapiens")
1074    print list(org.genes)[:10]
1075    org.gene_aliases
1076    print org.pathways(with_ids=org.genes.keys()[:5])
1077    print org.enzymes()
1078    print org.enriched_pathways(org.genes.keys()[:10])
1079    print org.genematcher
1080
1081if __name__=="__main__":
1082    test()
1083#    org1 = KEGGOrganism("ddi")
1084#    org2 = KEGGOrganism("ddi")
1085#    org2.api = KEGGInterface()
1086#    tests = [("get_genes", ()),
1087#             ("get_genes_by_enzyme", ("ec:1.1.1.1",)),
1088#             ("get_genes_by_pathway", ("path:ddi00010",)),
1089#             ("get_pathways_by_genes", (["ddi:DDB_0191256"],)),
1090#             ("get_pathways_by_enzymes", (["ec:1.1.1.1"],)),
1091#             ("get_pathways_by_compounds", (["cpd:C00001"],)),
1092#             ("get_linked_pathways", ("path:ddi00010",)),
1093#             ("list_pathways", ()),
1094#             ("get_compounds_by_enzyme", ("ec:1.1.1.1",)),
1095#             ("get_compounds_by_pathway", ("path:ddi00010",)),
1096#             ("get_enzymes_by_compound", ("cpd:C00001",)),
1097#             ("get_enzymes_by_pathway", ("path:ddi00010",)),
1098#             ("get_enzymes_by_gene", ("ddi:DDB_0191256",))]
1099#    for name, args in tests:
1100#        s1 = set(getattr(org1, name)(*args))
1101#        s2 = set(getattr(org2, name)(*args))
1102#        if s1 and s2:
1103#            print name
1104#            print s1-s2
1105#            print s2-s1
1106#        else:
1107#            print name
1108#            print "both empty"
Note: See TracBrowser for help on using the repository browser.