source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1770:48af55b8e628

Revision 1770:48af55b8e628, 60.3 KB checked in by Ales Erjavec <ales.erjavec@…>, 12 months ago (diff)

'Ontology.getitem' now raises a KeyError for an unknown term id.

Line 
1"""
2`obiGO` is a Gene Ontology (GO) Handling Library.
3
4"""
5
6from __future__ import absolute_import
7import os
8import sys
9import tarfile
10import gzip
11import re
12import cPickle
13import shutil
14import urllib2
15import warnings
16import copy
17
18from urllib import urlretrieve
19from gzip import GzipFile
20from datetime import datetime
21from collections import defaultdict
22
23from Orange.orng import orngEnviron
24from Orange.orng import orngServerFiles
25from Orange.orng import orngMisc
26
27from . import obiProb, obiGene
28
29default_database_path = os.path.join(orngServerFiles.localpath(), "GO")
30
31
32_CVS_REVISION_RE = re.compile(r"^(rev)?\d+\.\d+$")
33
34evidenceTypes = {
35# Experimental
36    'EXP': 'Inferred from Experiment',
37    'IDA': 'Inferred from Direct Assay',
38    'IPI': 'Inferred from Physical Interaction',  # [with <database:protein_name>]',
39    'IMP': 'Inferred from Mutant Phenotype',
40    'IGI': 'Inferred from Genetic Interaction',  # [with <database:gene_symbol[allele_symbol]>]',
41    'IEP': 'Inferred from Expression Pattern',
42# Computational Analysis Evidence Codes
43    'ISS': 'Inferred from Sequence Similarity',  # [with <database:sequence_id>] ',
44    'ISA': 'Inferred from Sequence Alignment',
45    'ISO': 'Inferred from Sequence Orthology',
46    'ISM': 'Inferred from Sequence Model',
47    'IGC': 'Inferred from Genomic Context',
48    'RCA': 'Inferred from Reviewed Computational Analysis',
49# Author Statement Evidence Codes
50    'TAS': 'Traceable author statement',
51    'NAS': 'Non-traceable author statement',
52# Curatorial Statement Evidence Codes
53    'IC': 'Inferred by curator',
54    'ND': 'No biological data available',
55# Computationally-assigned Evidence Codes
56    'IEA': 'Inferred from electronic annotation',  # [to <database:id>]',
57# Obsolete Evidence Codes
58    'NR': 'Not Recorded(Obsolete)'
59}
60
61
62evidenceDict = defaultdict(int, [(e, 2 ** i) for i, e in
63                                 enumerate(evidenceTypes.keys())])
64
65evidenceTypesOrdered = [
66'EXP',
67'IDA',
68'IPI',
69'IMP',
70'IGI',
71'IEP',
72# Computational Analysis Evidence Codes
73'ISS',
74'ISA',
75'ISO',
76'ISM',
77'IGC',
78'RCA',
79# Author Statement Evidence Codes
80'TAS',
81'NAS',
82# Curatorial Statement Evidence Codes
83'IC',
84'ND',
85# Computationally-assigned Evidence Codes
86'IEA',
87# Obsolete Evidence Codes
88'NR'
89]
90
91multiplicitySet = set(
92    ["alt_id", "is_a", "subset", "synonym", "related_synonym",
93     "exact_synonym", "broad_synonym", "narrow_synonym",
94     "xref_analog", "xref_unknown", "relationship"])
95
96multipleTagSet = multiplicitySet
97
98annotationFields = [
99    "DB", "DB_Object_ID", "DB_Object_Symbol", "Qualifier", "GO_ID",
100    "DB_Reference", "Evidence_Code", "With_From", "Aspect", "DB_Object_Name",
101    "DB_Object_Synonym", "DB_Object_Type", "Taxon", "Date", "Assigned_by"]
102
103annotationFieldsDict = {
104    "DB": 0,
105    "DB_Object_ID": 1,
106    "DB_Object_Symbol": 2,
107    "Qualifier": 3,
108    "GO_ID": 4,
109    "GO ID": 4,
110    "GOID": 4,
111    "DB_Reference": 5,
112    "DB:Reference": 5,
113    "Evidence_Code": 6,
114    "Evidence Code": 6,
115    "Evidence_code": 6,  # compatibility with  older revisions
116    "With_or_From": 7,
117    "With (or) From": 7,
118    "Aspect": 8,
119    "DB_Object_Name": 9,
120    "DB_Object_Synonym": 10,
121    "DB_Object_Type": 11,
122    "taxon": 12,
123    "Date": 13,
124    "Assigned_by": 14
125}
126
127builtinOBOObjects = ["""
128[Typedef]
129id: is_a
130name: is_a
131range: OBO:TERM_OR_TYPE
132domain: OBO:TERM_OR_TYPE
133definition: The basic subclassing relationship [OBO:defs]"""
134,
135"""[Typedef]
136id: disjoint_from
137name: disjoint_from
138range: OBO:TERM
139domain: OBO:TERM
140definition: Indicates that two classes are disjoint [OBO:defs]"""
141,
142"""[Typedef]
143id: instance_of
144name: instance_of
145range: OBO:TERM
146domain: OBO:INSTANCE
147definition: Indicates the type of an instance [OBO:defs]"""
148,
149"""[Typedef]
150id: inverse_of
151name: inverse_of
152range: OBO:TYPE
153domain: OBO:TYPE
154definition: Indicates that one relationship type is the inverse of another [OBO:defs]"""
155,
156"""[Typedef]
157id: union_of
158name: union_of
159range: OBO:TERM
160domain: OBO:TERM
161definition: Indicates that a term is the union of several others [OBO:defs]"""
162,
163"""[Typedef]
164id: intersection_of
165name: intersection_of
166range: OBO:TERM
167domain: OBO:TERM
168definition: Indicates that a term is the intersection of several others [OBO:defs]"""]
169
170
171class OBOObject(object):
172    """Represents a generic OBO object (e.g. Term, Typedef, Instance, ...)
173    Example:
174
175    >>> OBOObject(r"[Term]\nid: FOO:001\nname: bar", ontology)
176
177    """
178    _INTERN_TAGS = ["id", "name", "namespace", "alt_id", "is_a"]
179
180    def __init__(self, stanza=None, ontology=None):
181        self.ontology = ontology
182        self._lines = []
183        self.values = {}
184        self.related = set()
185        self.relatedTo = set()
186        if stanza:
187            self.ParseStanza(stanza)
188
189    def ParseStanza(self, stanza):
190        intern_tags = set(self._INTERN_TAGS)
191        for line in stanza.splitlines():
192            if ":" not in line:
193                continue
194            tag, rest = line.split(":", 1)
195            value, modifiers, comment = "", "", ""
196            if "!" in rest:
197                rest, comment = rest.split("!")
198            if "{" in rest:
199                value, modifiers = rest.split("{", 1)
200                modifiers = modifiers.strip("}")
201            else:
202                value = rest
203            tag = intern(tag)
204            value = value.strip()
205            comment = comment.strip()
206            if tag in intern_tags:
207                value, comment = intern(value), intern(comment)
208            self._lines.append((tag, value, modifiers, comment))
209            if tag in multipleTagSet:
210                self.values.setdefault(tag, []).append(value)
211            else:
212                self.values[tag] = value
213        self.related = set(self.GetRelatedObjects())
214        self.__dict__.update(self.values)
215        if "def" in self.__dict__:
216            self.__dict__["def_"] = self.__dict__["def"]
217
218    def GetRelatedObjects(self):
219        """Return a list of tuple pairs where the first element is relationship
220        typeId and the second id of object to whom the relationship applies to.
221
222        """
223        # TODO: add other defined Typedef ids
224        typeIds = [intern("is_a")]
225        result = [(typeId, id) for typeId in typeIds
226                  for id in self.values.get(typeId, [])]
227        result = result + [tuple(map(intern, r.split(None, 1)))
228                           for r in self.values.get("relationship", [])]
229        return result
230
231    def __repr__(self):
232        """ Return a string representation of the object in OBO format
233        """
234        repr = "[%s]\n" % type(self).__name__
235        for tag, value, modifiers, comment in self._lines:
236            repr = repr + tag + ": " + value
237            if modifiers:
238                repr = repr + "{ " + modifiers + " }"
239            if comment:
240                repr = repr + " ! " + comment
241            repr = repr + "\n"
242        return repr
243
244    def __str__(self):
245        """ Return the OBO object id entry
246        """
247        return "%s: %s" % (self.id, self.name)
248
249    def __iter__(self):
250        """ Iterates over sub terms
251        """
252        for typeId, id in self.relatedTo:
253            yield (typeId, self.ontology[id])
254
255
256class Term(OBOObject):
257    pass
258
259
260class Typedef(OBOObject):
261    pass
262
263
264class Instance(OBOObject):
265    pass
266
267
268class Ontology(object):
269    """ Ontology is the main class representing a gene ontology.
270
271    Example::
272        >>> ontology = Ontology("my_gene_ontology.obo")
273
274    """
275    version = 1
276
277    def __init__(self, file=None, progressCallback=None, rev=None):
278        """ Initialize the ontology from file (if `None` the default gene
279        ontology will be loaded). The optional `progressCallback` will be
280        called with a single argument to report on the progress.
281
282        """
283        self.terms = {}
284        self.typedefs = {}
285        self.instances = {}
286        self.slimsSubset = set()
287        if isinstance(file, basestring):
288            self.ParseFile(file, progressCallback)
289
290        elif rev is not None:
291            if not _CVS_REVISION_RE.match(rev):
292                raise ValueError("Invalid revision format.")
293            if rev.startswith("rev"):
294                rev = rev[3:]
295            pc = lambda v: progressCallback(v / 2.0) \
296                 if progressCallback else None
297            filename = os.path.join(default_database_path,
298                                    "gene_ontology_edit@rev%s.obo" % rev)
299            if not os.path.exists(filename):
300                self.DownloadOntologyAtRev(rev, filename, pc)
301            self.ParseFile(filename, lambda v: progressCallback(v / 2.0 + 50) \
302                                     if progressCallback else None)
303        else:
304            fool = self.Load(progressCallback)
305            # A fool and his attributes are soon parted
306            self.__dict__ = fool.__dict__
307
308    @classmethod
309    def Load(cls, progressCallback=None):
310        """ A class method that tries to load the ontology file from
311        default_database_path. It looks for a filename starting with
312        'gene_ontology'. If not found it will download it.
313        """
314        filename = os.path.join(default_database_path,
315                                "gene_ontology_edit.obo.tar.gz")
316        if not os.path.isfile(filename) and not os.path.isdir(filename):
317            orngServerFiles.download("GO", "gene_ontology_edit.obo.tar.gz")
318        try:
319            return cls(filename, progressCallback=progressCallback)
320        except (IOError, OSError), ex:
321            raise Exception("Could not locate ontology file")
322
323    def ParseFile(self, file, progressCallback=None):
324        """ Parse the file. file can be a filename string or an open filelike
325        object. The optional progressCallback will be called with a single
326        argument to report on the progress.
327        """
328        if type(file) == str:
329            if os.path.isfile(file) and tarfile.is_tarfile(file):
330                f = tarfile.open(file).extractfile("gene_ontology_edit.obo")
331            elif os.path.isfile(file):
332                f = open(file)
333            else:
334                f = open(os.path.join(file, "gene_ontology_edit.obo"))
335        else:
336            f = file
337
338        data = f.readlines()
339        data = "".join([line for line in data if not line.startswith("!")])
340        self.header = data[: data.index("[Term]")]
341        c = re.compile("\[.+?\].*?\n\n", re.DOTALL)
342        data = c.findall(data)
343
344        milestones = orngMisc.progressBarMilestones(len(data), 90)
345        for i, block in enumerate(builtinOBOObjects + data):
346            if block.startswith("[Term]"):
347                term = Term(block, self)
348                self.terms[term.id] = term
349            elif block.startswith("[Typedef]"):
350                typedef = Typedef(block, self)
351                self.typedefs[typedef.id] = typedef
352            elif block.startswith("[Instance]"):
353                instance = Instance(block, self)
354                self.instances[instance.id] = instance
355            if progressCallback and i in milestones:
356                progressCallback(90.0 * i / len(data))
357
358        self.aliasMapper = {}
359        self.reverseAliasMapper = defaultdict(set)
360        milestones = orngMisc.progressBarMilestones(len(self.terms), 10)
361        for i, (id, term) in enumerate(self.terms.iteritems()):
362            for typeId, parent in term.related:
363                self.terms[parent].relatedTo.add((typeId, id))
364            try:
365                self.aliasMapper.update([(alt_id, id)
366                                         for alt_id in term.alt_id])
367                self.reverseAliasMapper[id].union_update(term.alt_id)
368            except AttributeError:
369                pass
370            if progressCallback and i in milestones:
371                progressCallback(90.0 + 10.0 * i / len(self.terms))
372
373    def GetDefinedSlimsSubsets(self):
374        """ Return a list of defined subsets.
375        """
376        return [line.split()[1] for line in self.header.splitlines()
377                if line.startswith("subsetdef:")]
378
379    def SetSlimsSubset(self, subset):
380        """ Set the slims term subset to subset. If subset is a string it
381        must equal one of the defined subsetdef.
382        """
383        if type(subset) == str:
384            self.slimsSubset = [id for id, term in self.terms.items()
385                                if subset in getattr(term, "subset", set())]
386        else:
387            self.slimsSubset = set(subset)
388
389    def GetSlimsSubset(self, subset):
390        return [id for id, term in self.terms.items()
391                if subset in getattr(term, "subset", set())]
392
393    def GetSlimTerms(self, termId):
394        """ Return a list of slim terms for termId.
395        """
396        queue = set([termId])
397        visited = set()
398        slims = set()
399        while queue:
400            term = queue.pop()
401            visited.add(term)
402            if term in self.slimsSubset:
403                slims.add(term)
404            else:
405                queue.update(set(id for typeId, id in self[term].related) -
406                             visited)
407        return slims
408
409    def ExtractSuperGraph(self, terms):
410        """ Return all super terms of terms up to the most general one.
411        """
412        terms = [terms] if type(terms) == str else terms
413        visited = set()
414        queue = set(terms)
415        while queue:
416            term = queue.pop()
417            visited.add(term)
418            queue.update(set(id for typeId, id in self[term].related) -
419                         visited)
420        return visited
421
422    def ExtractSubGraph(self, terms):
423        """ Return all sub terms of terms.
424        """
425        terms = [terms] if type(terms) == str else terms
426        visited = set()
427        queue = set(terms)
428        while queue:
429            term = queue.pop()
430            visited.add(term)
431            queue.update(set(id for typeId, id in self[term].relatedTo) -
432                         visited)
433        return visited
434
435    def GetTermDepth(self, term, cache_={}):
436        """ Return the minimum depth of a term (length of the shortest
437        path to this term from the top level term).
438        """
439        if term not in cache_:
440            cache_[term] = min([self.GetTermDepth(parent) + 1
441                                for typeId, parent in self[term].related] or
442                               [1])
443        return cache_[term]
444
445    def __getitem__(self, id):
446        """
447        Return a object with id.
448        """
449        if id in self.terms:
450            return self.terms[id]
451        elif id in self.aliasMapper:
452            return self.terms[self.aliasMapper[id]]
453        else:
454            raise KeyError(id)
455
456    def __iter__(self):
457        """
458        Iterate over all ids in ontology.
459        """
460        return iter(self.terms)
461
462    def __len__(self):
463        """
464        Return number of objects in ontology.
465        """
466        return len(self.terms)
467
468    def __contains__(self, id):
469        """
470        Return `True` if a term with `id` is present in the ontology.
471        """
472        return id in self.terms or id in self.aliasMapper
473
474    @staticmethod
475    def DownloadOntology(file, progressCallback=None):
476        tFile = tarfile.open(file, "w:gz") if isinstance(file, basestring) \
477                else file
478        tmpDir = os.path.join(orngEnviron.bufferDir, "tmp_go/")
479        try:
480            os.mkdir(tmpDir)
481        except Exception:
482            pass
483        urlretrieve("http://www.geneontology.org/ontology/gene_ontology_edit.obo",
484                    os.path.join(tmpDir, "gene_ontology_edit.obo"),
485                    progressCallback and __progressCallbackWrapper(progressCallback))
486        tFile.add(os.path.join(tmpDir, "gene_ontology_edit.obo"),
487                  "gene_ontology_edit.obo")
488        tFile.close()
489        os.remove(os.path.join(tmpDir, "gene_ontology_edit.obo"))
490
491    @staticmethod
492    def DownloadOntologyAtRev(rev, filename=None, progressCallback=None):
493        url = "http://cvsweb.geneontology.org/cgi-bin/cvsweb.cgi/~checkout~/go/ontology/gene_ontology_edit.obo?rev=%s" % rev
494        url += ";content-type=text%2Fplain"
495        if filename is None:
496            filename = os.path.join(default_database_path,
497                                    "gene_ontology_edit@rev%s.obo" % rev)
498        r = urllib2.urlopen(url)
499
500        with open(filename + ".part", "wb") as f:
501            shutil.copyfileobj(r, f)
502
503        os.rename(filename + ".part", filename)
504
505
506_re_obj_name_ = re.compile("([a-zA-z0-9-_]+)")
507
508
509class AnnotationRecord(object):
510    """Holds the data for an annotation record read from the annotation file.
511    Fields can be accessed with the names: DB, DB_Object_ID, DB_Object_Symbol,
512    Qualifier, GO_ID, DB_Reference, Evidence_code, With_or_From, Aspect,
513    DB_Object_Name, DB_Object_Synonym, DB_Object_Type, taxon, Date,
514    Assigned_by (e.g. rec.GO_ID) or by supplying the original name of the
515    field (see http://geneontology.org/GO.annotation.shtml#file) to the get
516    method (e.g. rec.get("GO ID")) The object also provides the following
517    data members for quicker access: geneName, GOId, evidence, aspect and
518    alias(a list of aliases)
519
520    """
521    __slots__ = annotationFields + ["geneName", "GOId", "evidence",
522                                    "aspect", "alias", "additionalAliases"]
523
524    def __init__(self, fullText):
525        """\
526        :param fulText: A single line from the annotation file.
527
528        """
529        for slot, val in zip(annotationFields, fullText.split("\t")):
530            setattr(self, slot, intern(val))
531        self.geneName = self.DB_Object_Symbol
532        self.GOId = self.GO_ID
533        self.evidence = self.Evidence_Code
534        self.aspect = self.Aspect
535        self.alias = list(map(intern, self.DB_Object_Synonym.split("|")))
536
537        self.additionalAliases = []
538        if ":" in self.DB_Object_Name:
539            self.additionalAliases = []  # _re_obj_name_.findall(self.DB_Object_Name.split(":")[0])
540
541    def __getattr__(self, name):
542        if name in annotationFieldsDict:
543            return getattr(self, self.__slots__[annotationFieldsDict[name]])
544        else:
545            raise AttributeError(name)
546
547
548class Annotations(object):
549    """Annotations object holds the annotations.
550    """
551    version = 2
552
553    def __init__(self, file=None, ontology=None, genematcher=None,
554                 progressCallback=None, rev=None):
555        """Initialize the annotations from file by calling ParseFile on it.
556        The ontology must be an instance of Ontology class.
557        The optional progressCallback will be called with a single argument
558        to report on the progress.
559
560        """
561        self.file = file
562        self.ontology = ontology
563        self.allAnnotations = defaultdict(list)
564        self.geneAnnotations = defaultdict(list)
565        self.termAnnotations = defaultdict(list)
566        self._geneNames = None
567        self._geneNamesDict = None
568        self._aliasMapper = None
569        self.additionalAliases = {}
570        self.annotations = []
571        self.header = ""
572        self.genematcher = genematcher
573        self.taxid = None
574        if type(file) in [list, set, dict, Annotations]:
575            for ann in file:
576                self.AddAnnotation(ann)
577            if type(file, Annotations):
578                taxid = file.taxid
579        elif isinstance(file, basestring) and os.path.exists(file):
580            self.ParseFile(file, progressCallback)
581            try:
582                self.taxid = to_taxid(os.path.basename(file).split(".")[1]).pop()
583            except IOError:
584                pass
585        elif file is not None:
586            # Organism code
587            if rev is not None:
588                if not _CVS_REVISION_RE.match(rev):
589                    raise ValueError("Invalid revision format")
590
591                if rev.startswith("rev"):
592                    rev = rev[3:]
593                code = self.organism_name_search(file)
594                filename = os.path.join(default_database_path,
595                                        "gene_association.%s@rev%s.tar.gz" %
596                                        (code, rev))
597
598                if not os.path.exists(filename):
599                    self.DownloadAnnotationsAtRev(code, rev, filename,
600                                                  progressCallback
601                                                  )
602
603                self.ParseFile(filename, progressCallback)
604                self.taxid = to_taxid(code).pop()
605            else:
606                a = self.Load(file, ontology, genematcher, progressCallback)
607                self.__dict__ = a.__dict__
608                self.taxid = to_taxid(organism_name_search(file)).pop()
609        if not self.genematcher and self.taxid:
610            self.genematcher = obiGene.matcher(
611                [obiGene.GMGO(self.taxid)] +
612                ([obiGene.GMDicty(),
613                  [obiGene.GMGO(self.taxid), obiGene.GMDicty()]]
614                 if self.taxid == "352472" else [])
615            )
616        if self.genematcher:
617            self.genematcher.set_targets(self.geneNames)
618
619    @classmethod
620    def organism_name_search(cls, org):
621        ids = to_taxid(org)
622        from . import obiTaxonomy as tax
623        if not ids:
624            ids = [org] if org in Taxonomy().common_org_map.keys() + \
625                  Taxonomy().code_map.keys() else []
626        if not ids:
627            ids = tax.to_taxid(org, mapTo=Taxonomy().keys())
628        if not ids:
629            ids = tax.search(org, exact=True)
630            ids = set(ids).intersection(Taxonomy().keys())
631        if not ids:
632            ids = tax.search(org)
633            ids = set(ids).intersection(Taxonomy().keys())
634        codes = set([from_taxid(id) for id in ids])
635        if len(codes) > 1:
636            raise tax.MultipleSpeciesException(
637                ", ".join(["%s: %s" % (str(from_taxid(id)), tax.name(id))
638                           for id in ids]))
639        elif len(codes) == 0:
640            raise tax.UnknownSpeciesIdentifier(org)
641        return codes.pop()
642
643    @classmethod
644    def organism_version(cls, name):
645        name = organism_name_search(name)
646        orngServerFiles.localpath_download(
647            "GO", "gene_association.%s.tar.gz" % name)
648        return ("v%i." % cls.version) + orngServerFiles.info("GO",
649                        "gene_association.%s.tar.gz" % name)["datetime"]
650
651    def SetOntology(self, ontology):
652        """ Set the ontology to use in the annotations mapping.
653        """
654        self.allAnnotations = defaultdict(list)
655        self._ontology = ontology
656
657    def GetOntology(self):
658        return self._ontology
659
660    ontology = property(GetOntology, SetOntology,
661                        doc="Ontology object for annotations")
662
663    @classmethod
664    def Load(cls, org, ontology=None, genematcher=None, progressCallback=None):
665        """A class method that tries to load the association file for the
666        given organism from default_database_path.
667        """
668        code = organism_name_search(org)
669
670        file = "gene_association.%s.tar.gz" % code
671
672        path = os.path.join(orngServerFiles.localpath("GO"), file)
673
674        if not os.path.exists(path):
675            sf = orngServerFiles.ServerFiles()
676            available = sf.listfiles("GO")
677            if file not in available:
678                from . import obiKEGG
679                raise obiKEGG.OrganismNotFoundError(org + str(code))
680            orngServerFiles.download("GO", file)
681
682        return cls(path, ontology=ontology, genematcher=genematcher,
683                   progressCallback=progressCallback)
684
685    def ParseFile(self, file, progressCallback=None):
686        """ Parse and load the annotations from file. Report progress
687        with progressCallback.
688        File can be:
689            - a tarball containing the association file named gene_association
690            - a directory name containing the association file named gene_association
691            - a path to the actual association file
692            - an open file-like object of the association file
693        """
694        if type(file) == str:
695            if os.path.isfile(file) and tarfile.is_tarfile(file):
696                f = tarfile.open(file).extractfile("gene_association")
697            elif os.path.isfile(file) and file.endswith(".gz"):
698                f = gzip.open(file)
699            elif os.path.isfile(file):
700                f = open(file)
701            else:
702                f = open(os.path.join(file, "gene_association"))
703        else:
704            f = file
705        lines = [line for line in f.read().splitlines() if line.strip()]
706
707        milestones = orngMisc.progressBarMilestones(len(lines), 100)
708        for i, line in enumerate(lines):
709            if line.startswith("!"):
710                self.header = self.header + line + "\n"
711                continue
712
713            a = AnnotationRecord(line)
714            self.AddAnnotation(a)
715#            self.annotations.append(a)
716            if progressCallback and i in milestones:
717                progressCallback(100.0 * i / len(lines))
718
719    def AddAnnotation(self, a):
720        """ Add a single `AnotationRecord` instance to this `Annotations`
721        object.
722        """
723        if not isinstance(a, AnnotationRecord):
724            a = AnnotationRecord(a)
725        if not a.geneName or not a.GOId or a.Qualifier == "NOT":
726            return
727
728        self.geneAnnotations[a.geneName].append(a)
729        self.annotations.append(a)
730        self.termAnnotations[a.GOId].append(a)
731        self.allAnnotations = defaultdict(list)
732
733        self._geneNames = None
734        self._geneNamesDict = None
735        self._aliasMapper = None
736
737    @property
738    def geneNamesDict(self):
739        if getattr(self, "_geneNamesDict", None) is None:
740            self._geneNamesDict = defaultdict(set)
741            for alias, name in self.aliasMapper.iteritems():
742                self._geneNamesDict[name].add(alias)
743        return self._geneNamesDict
744
745    @property
746    def geneNames(self):
747        if getattr(self, "_geneNames", None) is None:
748            self._geneNames = set([ann.geneName for ann in self.annotations])
749        return self._geneNames
750
751    @property
752    def aliasMapper(self):
753        if getattr(self, "_aliasMapper", None) is None:
754            self._aliasMapper = {}
755            for ann in self.annotations:
756                self._aliasMapper.update([(alias, ann.geneName)
757                                          for alias in ann.alias +
758                                           [ann.geneName, ann.DB_Object_ID]])
759        return self._aliasMapper
760
761    def GetGeneNamesTranslator(self, genes):
762        """ Return a dictionary mapping canonical names (DB_Object_Symbol)
763        to `genes`.
764
765        """
766        def alias(gene):
767            if self.genematcher:
768                return self.genematcher.umatch(gene)
769            else:
770                return gene if gene in self.geneNames else \
771                        self.aliasMapper.get(gene,
772                             self.additionalAliases.get(gene, None))
773
774        return dict([(alias(gene), gene) for gene in genes if alias(gene)])
775
776    def _CollectAnnotations(self, id, visited):
777        """ Recursive function collects and caches all annotations for id
778        """
779        if id not in self.allAnnotations and id not in visited:
780            if id in self.ontology.reverseAliasMapper:
781                annotations = [self.termAnnotations.get(alt_id, [])
782                               for alt_id in
783                               self.ontology.reverseAliasMapper[id]] + \
784                              [self.termAnnotations[id]]
785            else:
786                ## annotations for this term alone
787                annotations = [self.termAnnotations[id]]
788            visited.add(id)
789            for typeId, child in self.ontology[id].relatedTo:
790                aa = self._CollectAnnotations(child, visited)
791                if type(aa) == set:
792                    ## if it was already reduced in GetAllAnnotations
793                    annotations.append(aa)
794                else:
795                    annotations.extend(aa)
796            self.allAnnotations[id] = annotations
797        return self.allAnnotations[id]
798
799    def GetAllAnnotations(self, id):
800        """ Return a set of all annotations (instances if `AnnotationRectord`)
801        for GO term `id` and all it's subterms.
802
803        :param id: GO term id
804        :type id: str
805
806        """
807        visited = set()
808        id = self.ontology.aliasMapper.get(id, id)
809        if id not in self.allAnnotations or \
810                type(self.allAnnotations[id]) == list:
811            annot_set = set()
812            for annots in self._CollectAnnotations(id, set()):
813                annot_set.update(annots)
814            self.allAnnotations[id] = annot_set
815        return self.allAnnotations[id]
816
817    def GetAllGenes(self, id, evidenceCodes=None):
818        """ Return a list of genes annotated by specified `evidenceCodes`
819        to GO term 'id' and all it's subterms."
820
821        :param id: GO term id
822        :type id: str
823
824        :param evidneceCodes: List of evidence codes to consider when
825                              matching annotations to terms.
826        :type evidenceCodes: list-of-strings
827        """
828        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
829        annotations = self.GetAllAnnotations(id)
830        return list(set([ann.geneName for ann in annotations
831                         if ann.Evidence_Code in evidenceCodes]))
832
833    def GetEnrichedTerms(self, genes, reference=None, evidenceCodes=None,
834                         slimsOnly=False, aspect=None, prob=obiProb.Binomial(),
835                         useFDR=True, progressCallback=None):
836        """ Return a dictionary of enriched terms, with tuples of
837        (list_of_genes, p_value, reference_count) for items and term
838        ids as keys. P-Values are FDR adjusted if useFDR is True (default).
839
840        :param genes: List of genes
841        :param reference: list of genes (if None all genes included in the
842                          annotations will be used).
843        :param evidenceCodes: List of evidence codes to consider.
844        :param slimsOnly: If `True` return only slim terms
845        :param aspect: Which aspects to use. Use all by default. "P", "F", "C"
846            or a set containing these elements.
847        """
848        revGenesDict = self.GetGeneNamesTranslator(genes)
849        genes = set(revGenesDict.keys())
850        if reference:
851            refGenesDict = self.GetGeneNamesTranslator(reference)
852            reference = set(refGenesDict.keys())
853        else:
854            reference = self.geneNames
855
856        if aspect == None:
857            aspects_set = set(["P", "C", "F"])
858        elif isinstance(aspect, basestring):
859            aspects_set = set([aspect])
860        else:
861            aspects_set = aspect
862
863        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
864        annotations = [ann
865                       for gene in genes for ann in self.geneAnnotations[gene]
866                       if ann.Evidence_Code in evidenceCodes and
867                       ann.Aspect in aspects_set]
868
869        refAnnotations = set(
870            [ann
871             for gene in reference for ann in self.geneAnnotations[gene]
872             if ann.Evidence_Code in evidenceCodes and
873             ann.Aspect in aspects_set]
874        )
875
876        annotationsDict = defaultdict(set)
877        for ann in annotations:
878            annotationsDict[ann.GO_ID].add(ann)
879
880        if slimsOnly and not self.ontology.slimsSubset:
881            warnings.warn("Unspecified slims subset in the ontology! "
882                          "Using 'goslim_generic' subset", UserWarning)
883            self.ontology.SetSlimsSubset("goslim_generic")
884
885        terms = annotationsDict.keys()
886        filteredTerms = [term for term in terms if term in self.ontology]
887
888        if len(terms) != len(filteredTerms):
889            termDiff = set(terms) - set(filteredTerms)
890            warnings.warn("%s terms in the annotations were not found in the "
891                          "ontology." % ",".join(map(repr, termDiff)),
892                          UserWarning)
893
894        terms = self.ontology.ExtractSuperGraph(filteredTerms)
895        res = {}
896
897        milestones = orngMisc.progressBarMilestones(len(terms), 100)
898        for i, term in enumerate(terms):
899            if slimsOnly and term not in self.ontology.slimsSubset:
900                continue
901            allAnnotations = self.GetAllAnnotations(term).intersection(refAnnotations)
902##            allAnnotations.intersection_update(refAnnotations)
903            allAnnotatedGenes = set([ann.geneName for ann in allAnnotations])
904            mappedGenes = genes.intersection(allAnnotatedGenes)
905##            if not mappedGenes:
906##                print >> sys.stderr, term, sorted(genes)
907##                print >> sys.stderr, sorted(allAnnotatedGenes)
908##                return
909            if len(reference) > len(allAnnotatedGenes):
910                mappedReferenceGenes = reference.intersection(allAnnotatedGenes)
911            else:
912                mappedReferenceGenes = allAnnotatedGenes.intersection(reference)
913            res[term] = ([revGenesDict[g] for g in mappedGenes],
914                         prob.p_value(len(mappedGenes), len(reference),
915                                      len(mappedReferenceGenes), len(genes)),
916                         len(mappedReferenceGenes))
917            if progressCallback and i in milestones:
918                progressCallback(100.0 * i / len(terms))
919        if useFDR:
920            res = sorted(res.items(), key=lambda (_1, (_2, p, _3)): p)
921            res = dict([(id, (genes, p, ref))
922                        for (id, (genes, _, ref)), p in
923                        zip(res, obiProb.FDR([p for _, (_, p, _) in res]))])
924        return res
925
926    def GetAnnotatedTerms(self, genes, directAnnotationOnly=False,
927                          evidenceCodes=None, progressCallback=None):
928        """ Return all terms that are annotated by genes with evidenceCodes.
929        """
930        genes = [genes] if type(genes) == str else genes
931        revGenesDict = self.GetGeneNamesTranslator(genes)
932        genes = set(revGenesDict.keys())
933        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
934        annotations = [ann for gene in genes for ann in self.geneAnnotations[gene]
935                       if ann.Evidence_Code in evidenceCodes]
936        dd = defaultdict(set)
937        for ann in annotations:
938            dd[ann.GO_ID].add(revGenesDict.get(ann.geneName, ann.geneName))
939        if not directAnnotationOnly:
940            terms = dd.keys()
941            filteredTerms = [term for term in terms if term in self.ontology]
942            if len(terms) != len(filteredTerms):
943                termDiff = set(terms) - set(filteredTerms)
944                warnings.warn(
945                    "%s terms in the annotations were not found in the "
946                    "ontology." % ",".join(map(repr, termDiff)),
947                    UserWarning)
948
949            terms = self.ontology.ExtractSuperGraph(filteredTerms)
950            for i, term in enumerate(terms):
951                termAnnots = self.GetAllAnnotations(term).intersection(annotations)
952##                termAnnots.intersection_update(annotations)
953                dd[term].update([revGenesDict.get(ann.geneName, ann.geneName)
954                                 for ann in termAnnots])
955        return dict(dd)
956
957    def DrawEnrichmentGraph(self, terms, clusterSize, refSize=None,
958                            file="graph.png", width=None, height=None,
959                            precison=3):
960        refSize = len(self.geneNames) if refSize == None else refSize
961        sortedterms = sorted(terms.items(), key=lambda term: term[1][1])
962        fdr = dict(zip([t[0] for t in sortedterms],
963                       obiProb.FDR([t[1][1] for t in sortedterms])))
964        termsList = [(term,
965                      ((float(len(terms[term][0])) / clusterSize) /
966                       (float(terms[term][2]) / refSize)),
967                      len(terms[term][0]),
968                      terms[term][2],
969                      terms[term][1],
970                      fdr[term],
971                      terms[term][0])
972                     for term in terms]
973
974        drawEnrichmentGraph(termsList, file, width, height,
975                            ontology=self.ontology, precison=precison)
976
977    def __add__(self, iterable):
978        """ Return a new Annotations object with combined annotations
979        """
980        return Annotations([a for a in self] + [a for a in iterable],
981                           ontology=self.ontology)
982
983    def __iadd__(self, iterable):
984        """ Add annotations to this instance
985        """
986        self.extend(iterable)
987        return self
988
989    def __contains__(self, item):
990        return item in self.annotations
991
992    def __iter__(self):
993        """ Iterate over all AnnotationRecord objects in annotations
994        """
995        return iter(self.annotations)
996
997    def __len__(self):
998        """ Return the number of annotations
999        """
1000        return len(self.annotations)
1001
1002    def __getitem__(self, index):
1003        """ Return the i-th annotation record
1004        """
1005        return self.annotations[index]
1006
1007    def __getslice__(self, *args):
1008        return self.annotations.__getslice__(*args)
1009
1010    def add(self, line):
1011        """ Add one annotation
1012        """
1013        self.AddAnnotation(line)
1014
1015    def append(self, line):
1016        """ Add one annotation
1017        """
1018        self.AddAnnotation(line)
1019
1020    def extend(self, lines):
1021        """ Add multiple annotations
1022        """
1023        for line in lines:
1024            self.AddAnnotation(line)
1025
1026    def RemapGenes(self, map):
1027        """
1028        """
1029        for gene in map:
1030            annotations = self.geneAnnotations[gene]
1031            for ann in annotations:
1032                for name in map[gene]:
1033                    ann1 = copy.copy(ann)
1034                    ann1.geneName = name
1035                    self.add(ann1)
1036        self.genematcher = obiGene.GMDirect()
1037        try:
1038            del self._geneNames
1039        except Exception:
1040            pass
1041        self.genematcher.set_targets(self.geneNames)
1042
1043    @staticmethod
1044    def DownloadAnnotations(org, file, progressCallback=None):
1045        tFile = tarfile.open(file, "w:gz") if type(file) == str else file
1046        tmpDir = os.path.join(orngEnviron.bufferDir, "tmp_go/")
1047        try:
1048            os.mkdir(tmpDir)
1049        except Exception:
1050            pass
1051        fileName = "gene_association." + org + ".gz"
1052        urlretrieve(("http://www.geneontology.org/gene-associations/" +
1053                     fileName),
1054                    os.path.join(tmpDir, fileName),
1055                    progressCallback and __progressCallbackWrapper(progressCallback))
1056        gzFile = GzipFile(os.path.join(tmpDir, fileName), "r")
1057        file = open(os.path.join(tmpDir, "gene_association." + org), "w")
1058        file.writelines(gzFile.readlines())
1059        file.flush()
1060        file.close()
1061
1062        tFile.add(os.path.join(tmpDir, "gene_association." + org),
1063                  "gene_association")
1064        annotation = Annotations(os.path.join(tmpDir, "gene_association." + org),
1065                    genematcher=obiGene.GMDirect(), progressCallback=progressCallback)
1066        cPickle.dump(annotation.geneNames, open(os.path.join(tmpDir, "gene_names.pickle"), "wb"))
1067        tFile.add(os.path.join(tmpDir, "gene_names.pickle"), "gene_names.pickle")
1068        tFile.close()
1069        os.remove(os.path.join(tmpDir, "gene_association." + org))
1070        os.remove(os.path.join(tmpDir, "gene_names.pickle"))
1071
1072    @staticmethod
1073    def DownloadAnnotationsAtRev(org, rev, filename=None,
1074                                 progressCallback=None):
1075        if filename is None:
1076            filename = os.path.join(default_database_path,
1077                                    "gene_association.%s@rev%s.tar.gz" %
1078                                    (code, rev))
1079        url = ("http://cvsweb.geneontology.org/cgi-bin/cvsweb.cgi/~checkout~/go/gene-associations/gene_association.%s.gz?rev=%s" %
1080               (org, rev))
1081        url += ";content-type=application%2Fx-gzip"
1082        r = urllib2.urlopen(url)
1083
1084        with open(filename + ".part", "wb") as f:
1085            shutil.copyfileobj(r, f)
1086
1087        os.rename(filename + ".part", filename)
1088
1089from .obiTaxonomy import pickled_cache
1090
1091
1092@pickled_cache(None, [("GO", "taxonomy.pickle"),
1093                      ("Taxonomy", "ncbi_taxonomy.tar.gz")])
1094def organism_name_search(name):
1095    return Annotations.organism_name_search(name)
1096
1097
1098def filterByPValue(terms, maxPValue=0.1):
1099    """ Filters the terms by the p-value. Asumes terms is a dict with
1100    the same structure as returned from GetEnrichedTerms.
1101
1102    """
1103    return dict(filter(lambda (k, e): e[1] <= maxPValue, terms.items()))
1104
1105
1106def filterByFrequency(terms, minF=2):
1107    """ Filters the terms by the cluster frequency. Asumes terms is
1108    a dict with the same structure as returned from GetEnrichedTerms.
1109
1110    """
1111    return dict(filter(lambda (k, e): len(e[0]) >= minF, terms.items()))
1112
1113
1114def filterByRefFrequency(terms, minF=4):
1115    """ Filters the terms by the reference frequency. Asumes terms is
1116    a dict with the same structure as returned from GetEnrichedTerms.
1117
1118    """
1119    return dict(filter(lambda (k, e): e[2] >= minF, terms.items()))
1120
1121
1122def drawEnrichmentGraph(enriched, file="graph.png", width=None, height=None,
1123                        header=None, ontology=None, precison=3):
1124    file = open(file, "wb") if type(file) == str else file
1125    drawEnrichmentGraph_tostreamMk2(enriched, file, width, height, header,
1126                                    ontology, precison)
1127
1128
1129def drawEnrichmentGraph_tostreamMk2(enriched, fh, width, height, header=None,
1130                                    ontology=None, precison=4):
1131    ontology = ontology if ontology else Ontology()
1132    header = header if header else ["List", "Total", "p-value", "FDR",
1133                                    "Names", "Genes"]
1134    GOTerms = dict([(t[0], t) for t in enriched if t[0] in ontology])
1135
1136    def getParents(term):
1137        parents = ontology.ExtractSuperGraph([term])
1138        parents = [id for id in parents if id in GOTerms and id != term]
1139        c = reduce(set.union, [set(ontology.ExtractSuperGraph([id])) -
1140                               set([id]) for id in parents], set())
1141        parents = [t for t in parents if t not in c]
1142        return parents
1143    parents = dict([(term, getParents(term)) for term in GOTerms])
1144    # print "Parentes", parents
1145
1146    def getChildren(term):
1147        return [id for id in GOTerms if term in parents[id]]
1148    topLevelTerms = [id for id in parents if not parents[id]]
1149    # print "Top level terms", topLevelTerms
1150    termsList = []
1151    fmt = "%" + ".%if" % precison
1152
1153    def collect(term, parent):
1154        termsList.append(GOTerms[term][1:4] + \
1155                         (fmt % GOTerms[term][4],
1156                          fmt % GOTerms[term][5],
1157                          ontology[term].name,
1158                          ", ".join(GOTerms[term][6])) + (parent,))
1159        parent = len(termsList) - 1
1160        for c in getChildren(term):
1161            collect(c, parent)
1162
1163    for topTerm in topLevelTerms:
1164        collect(topTerm, None)
1165    for entry in enriched:
1166        if entry[0] not in ontology:
1167            termsList.append(entry[1:4] + \
1168                             (fmt % entry[4],
1169                              fmt % entry[5],
1170                              entry[0],
1171                              ", ".join(entry[6])) + (None,))
1172
1173    drawEnrichmentGraphPIL_tostream(termsList, header, fh, width, height)
1174
1175
1176def drawEnrichmentGraphPIL_tostream(termsList, headers, fh, width=None, height=None):
1177    from PIL import Image, ImageDraw, ImageFont
1178    backgroundColor = (255, 255, 255)
1179    textColor = (0, 0, 0)
1180    graphColor = (0, 0, 255)
1181    fontSize = height == None and 12 or (height - 60) / len(termsList)
1182    font = ImageFont.load_default()
1183    try:
1184        font = ImageFont.truetype("arial.ttf", fontSize)
1185    except:
1186        pass
1187    getMaxTextHeightHint = lambda l: max([font.getsize(t)[1] for t in l])
1188    getMaxTextWidthHint = lambda l: max([font.getsize(t)[0] for t in l])
1189    maxFoldWidth = width != None and min(150, width / 6) or 150
1190    maxFoldEnrichment = max([t[0] for t in termsList])
1191    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
1192    foldWidths = [int(foldNormalizationFactor * term[0]) for term in termsList]
1193    treeStep = 10
1194    treeWidth = {}
1195    for i, term in enumerate(termsList):
1196        treeWidth[i] = (term[-1] == None and 1 or treeWidth[term[-1]] + 1)
1197    treeStep = width != None and min(treeStep, width / (6 * max(treeWidth.values())) or 2) or treeStep
1198    treeWidth = [w * treeStep + foldWidths[i] for i, w in treeWidth.items()]
1199    treeWidth = max(treeWidth) - maxFoldWidth
1200    verticalMargin = 10
1201    horizontalMargin = 10
1202##    print verticalMargin, maxFoldWidth, treeWidth
1203##    treeWidth = 100
1204    firstColumnStart = verticalMargin + maxFoldWidth + treeWidth + 10
1205    secondColumnStart = firstColumnStart + getMaxTextWidthHint([str(t[1]) for t in termsList] + [headers[0]]) + 2
1206    thirdColumnStart = secondColumnStart + getMaxTextWidthHint([str(t[2]) for t in termsList] + [headers[1]]) + 2
1207    fourthColumnStart = thirdColumnStart + getMaxTextWidthHint([str(t[3]) for t in termsList] + [headers[2]]) + 2
1208    fifthColumnStart = fourthColumnStart + getMaxTextWidthHint([str(t[4]) for t in termsList] + [headers[3]]) + 4
1209##    maxAnnotationTextWidth = width==None and getMaxTextWidthHint([str(t[4]) for t in termsList]+["Annotation"]) or (width - fourthColumnStart - verticalMargin) * 2 / 3
1210    maxAnnotationTextWidth = width == None and getMaxTextWidthHint([str(t[5]) for t in termsList] + [headers[4]]) or max((width - fifthColumnStart - verticalMargin) * 2 / 3, getMaxTextWidthHint([t[5] for t in termsList] + [headers[4]]))
1211    sixthColumnStart = fifthColumnStart + maxAnnotationTextWidth + 4
1212    maxGenesTextWidth = width == None and getMaxTextWidthHint([str(t[6]) for t in termsList] + [headers[5]]) or (width - fifthColumnStart - verticalMargin) / 3
1213
1214    legendHeight = font.getsize("1234567890")[1] * 2
1215    termHeight = font.getsize("A")[1]
1216##    print fourthColumnStart, maxAnnotationTextWidth, verticalMargin
1217    width = sixthColumnStart + maxGenesTextWidth + verticalMargin
1218    height = len(termsList) * termHeight + 2 * (legendHeight + horizontalMargin)
1219
1220    image = Image.new("RGB", (width, height), backgroundColor)
1221    draw = ImageDraw.Draw(image)
1222
1223    def truncText(text, maxWidth, append=""):
1224        # print getMaxTextWidthHint([text]), maxAnnotationTextWidth
1225        if getMaxTextWidthHint([text]) > maxWidth:
1226            while getMaxTextWidthHint([text + "..." + append]) > maxWidth and text:
1227                text = text[:-1]
1228            if text:
1229                text = text + "..." + append
1230            else:
1231                text = append
1232        return text
1233    currentY = horizontalMargin + legendHeight
1234    connectAtX = {}
1235    for i, term in enumerate(termsList):
1236        draw.line([(verticalMargin, currentY + termHeight / 2), (verticalMargin + foldWidths[i], currentY + termHeight / 2)], width=termHeight - 2, fill=graphColor)
1237        draw.text((firstColumnStart, currentY), str(term[1]), font=font, fill=textColor)
1238        draw.text((secondColumnStart, currentY), str(term[2]), font=font, fill=textColor)
1239        draw.text((thirdColumnStart, currentY), str(term[3]), font=font, fill=textColor)
1240        draw.text((fourthColumnStart, currentY), str(term[4]), font=font, fill=textColor)
1241##        annotText = width!=None and truncText(str(term[5]), maxAnnotationTextWidth, str(term[5])) or str(term[4])
1242        annotText = width != None and truncText(str(term[5]), maxAnnotationTextWidth)
1243        draw.text((fifthColumnStart, currentY), annotText, font=font, fill=textColor)
1244        genesText = width != None and truncText(str(term[6]), maxGenesTextWidth) or str(term[6])
1245        draw.text((sixthColumnStart, currentY), genesText, font=font, fill=textColor)
1246        lineEnd = term[-1] == None and firstColumnStart - 10 or connectAtX[term[-1]]
1247        draw.line([(verticalMargin + foldWidths[i] + 1, currentY + termHeight / 2), (lineEnd, currentY + termHeight / 2)], width=1, fill=textColor)
1248        if term[-1] != None:
1249            draw.line([(lineEnd, currentY + termHeight / 2), (lineEnd, currentY + termHeight / 2 - termHeight * (i - term[-1]))], width=1, fill=textColor)
1250        connectAtX[i] = lineEnd - treeStep
1251        currentY += termHeight
1252
1253    currentY = horizontalMargin
1254    draw.text((firstColumnStart, currentY), headers[0], font=font, fill=textColor)
1255    draw.text((secondColumnStart, currentY), headers[1], font=font, fill=textColor)
1256    draw.text((thirdColumnStart, currentY), headers[2], font=font, fill=textColor)
1257    draw.text((fourthColumnStart, currentY), headers[3], font=font, fill=textColor)
1258    draw.text((fifthColumnStart, currentY), headers[4], font=font, fill=textColor)
1259    draw.text((sixthColumnStart, currentY), headers[5], font=font, fill=textColor)
1260
1261    horizontalMargin = 0
1262    # draw.line([(verticalMargin, height - horizontalMargin - legendHeight), (verticalMargin + maxFoldWidth, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1263    draw.line([(verticalMargin, horizontalMargin + legendHeight), (verticalMargin + maxFoldWidth, horizontalMargin + legendHeight)], width=1, fill=textColor)
1264    maxLabelWidth = getMaxTextWidthHint([" " + str(i) for i in range(int(maxFoldEnrichment + 1))])
1265    numOfLegendLabels = max(int(maxFoldWidth / maxLabelWidth), 2)
1266    for i in range(numOfLegendLabels + 1):
1267        # draw.line([(verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight/2), (verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1268        # draw.text((verticalMargin + i*maxFoldWidth/10 - font.getsize(str(i))[0]/2, height - horizontalMargin - legendHeight/2), str(i), font=font, fill=textColor)
1269
1270        label = str(int(i * maxFoldEnrichment / numOfLegendLabels))
1271        draw.line([(verticalMargin + i * maxFoldWidth / numOfLegendLabels, horizontalMargin + legendHeight / 2), (verticalMargin + i * maxFoldWidth / numOfLegendLabels, horizontalMargin + legendHeight)], width=1, fill=textColor)
1272        draw.text((verticalMargin + i * maxFoldWidth / numOfLegendLabels - font.getsize(label)[0] / 2, horizontalMargin), label, font=font, fill=textColor)
1273
1274    image.save(fh)
1275
1276
1277def drawEnrichmentGraphPylab_tostream(termsList, headers, fh, width=None, height=None, show=True):
1278    from matplotlib import pyplot as plt
1279    from matplotlib.patches import Rectangle
1280
1281    maxFoldWidth = width != None and min(150, width / 6) or 150
1282    maxFoldEnrichment = max([t[0] for t in termsList])
1283    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
1284##    foldWidths = [int(foldNormalizationFactor*term[0]) for term in termsList]
1285    foldWidths = [term[0] for term in termsList]
1286    treeStep = maxFoldEnrichment * 0.05
1287    treeWidth = {}
1288
1289    for i, term in enumerate(termsList):
1290        treeWidth[i] = (term[-1] == None and treeStep or treeWidth[term[-1]] + treeStep)
1291    maxTreeWidth = max(treeWidth)
1292
1293    connectAt = {}
1294    cellText = []
1295    axes1 = plt.axes([0.1, 0.1, 0.2, 0.8])
1296    for i, line in enumerate(termsList):
1297        enrichment, n, m, p_val, fdr_val, name, genes, parent = line
1298        r = Rectangle((0, len(termsList) - i - 0.4), enrichment, 0.8)
1299        plt.gca().add_patch(r)
1300        plt.plot([enrichment, connectAt.get(parent, maxFoldEnrichment + maxTreeWidth)], [len(termsList) - i, len(termsList) - i], color="black")
1301        connectAt[i] = connectAt.get(parent, maxFoldEnrichment + maxTreeWidth) - treeStep
1302        if parent != None:
1303            plt.plot([connectAt.get(parent)] * 2, [len(termsList) - i, len(termsList) - parent], color="black")
1304        cellText.append((str(n), str(m), p_val, fdr_val, name, genes))
1305
1306##    from Orange.orng.orngClustering import TableTextLayout
1307##    text = TableTextLayout((maxFoldEnrichment*1.1, len(termsList)), cellText)
1308    from Orange.orng.orngClustering import TablePlot
1309    if True:
1310        axes2 = plt.axes([0.3, 0.1, 0.6, 0.8], sharey=axes1)
1311        axes2.set_axis_off()
1312        table = TablePlot((0, len(termsList)), axes=plt.gca())
1313        for i, line in enumerate(cellText):
1314            for j, text in enumerate(line):
1315                table.add_cell(i, j, width=len(text), height=1, text=text, loc="left", edgecolor="w", facecolor="w")
1316
1317        table.set_figure(plt.gcf())
1318        plt.gca().add_artist(table)
1319        plt.gca()._set_artist_props(table)
1320##    plt.text(3, 3, "\n".join(["\t".join(text) for text in cellText]))
1321
1322##    table = plt.table(cellText=cellText, colLabels=headers, loc="right")
1323##    table.set_transform(plt.gca().transData)
1324##
1325##    table.set_xy(20,20)
1326    plt.show()
1327
1328
1329class Taxonomy(object):
1330    """Maps NCBI taxonomy ids to corresponding GO organism codes
1331    """
1332    common_org_map = {"297284": "9913", "30523": "9913",  # Bos taurus
1333                      "5782": "352472", "44689": "352472", "366501": "352472",  # Dictyostelium discoideum
1334                      "83333": "562",  # Escherichia coli
1335                      "52545": "4530", "4532": "4530", "65489": "4530", "4533": "4530", "77588": "4530", "29689": "4530",
1336                      "4538": "4530", "40148": "4530", "29690": "4530", "110450": "4530", "4534": "4530", "83309": "4530",
1337                      "4528": "4530", "127571": "4530", "40149": "4530", "83307": "4530", "63629": "4530", "4536": "4530",
1338                      "4535": "4530", "4537": "4530", "65491": "4530", "83308": "4530", "4529": "4530", "4530": "4530",
1339                      "39946": "4530", "39947": "4530", "110451": "4530", "364100": "4530", "364099": "4530", "4539": "4530",
1340                      }
1341    code_map = {"3702": "tair",  # Arabidopsis thaliana
1342                "9913": "goa_cow",  # Bos taurus
1343                "6239": "wb",  # Caenorhabditis elegans
1344                "3055": None,  # Chlamydomonas reinhardtii
1345                "7955": "zfin",  # Danio rerio (zebrafish)
1346                "352472": "dictyBase",  # Dictyostelium discoideum
1347                "7227": "fb",  # Drosophila melanogaster
1348                "562": "ecocyc",  # Escherichia coli
1349                "11103": None,  # Hepatitis C virus
1350                "9606": "goa_human",  # Homo sapiens
1351                "10090": "mgi",  # Mus musculus
1352                "2104": None,  # Mycoplasma pneumoniae
1353                "4530": "gramene_oryza",  # Oryza sativa
1354                "5833": "GeneDB_Pfalciparum",  # Plasmodium falciparum
1355                "4754": None,  # Pneumocystis carinii
1356                "10116": "rgd",  # Rattus norvegicus
1357                "4932": "sgd",  # Saccharomyces cerevisiae
1358                "4896": "GeneDB_Spombe",  # Schizosaccharomyces pombe
1359                "31033": None,  # Takifugu rubripes
1360                "8355": None,  # Xenopus laevis
1361                "4577": None  # Zea mays
1362                }
1363    version = 1
1364    __shared_state = {"tax": None}
1365
1366    def __init__(self):
1367        self.__dict__ = self.__shared_state
1368        if not self.tax:
1369            path = orngServerFiles.localpath_download("GO", "taxonomy.pickle")
1370            if os.path.isfile(path):
1371                self.tax = cPickle.load(open(path, "rb"))
1372            else:
1373                orngServerFiles.download("GO", "taxonomy.pickle")
1374                self.tax = cPickle.load(open(path, "rb"))
1375
1376    def __getitem__(self, key):
1377        key = self.common_org_map.get(key, key)
1378        return self.code_map[key]
1379
1380    def keys(self):
1381        return list(set(self.common_org_map.keys() + self.code_map.keys()))
1382
1383
1384def from_taxid(id):
1385    """ Return a set of GO organism codes that correspond to NCBI taxonomy id
1386    """
1387    return Taxonomy()[id]
1388
1389
1390def to_taxid(db_code):
1391    """ Return a set of NCBI taxonomy ids from db_code GO organism annotations
1392    """
1393    r = [key for key, val in Taxonomy().code_map.items() if db_code == val]
1394    return set(r)
1395
1396
1397class __progressCallbackWrapper:
1398    def __init__(self, callback):
1399        self.callback = callback
1400
1401    def __call__(self, bCount, bSize, fSize):
1402        fSize = 10000000 if fSize == -1 else fSize
1403        self.callback(100 * bCount * bSize / fSize)
1404
1405from .obiGenomicsUpdate import Update as UpdateBase
1406
1407
1408class Update(UpdateBase):
1409    def __init__(self, local_database_path=None, progressCallback=None):
1410        UpdateBase.__init__(self, local_database_path or getDataDir(), progressCallback)
1411
1412    def CheckModified(self, addr, date=None):
1413        return date < self.GetLastModified(addr) if date else True
1414
1415    def CheckModifiedOrg(self, org):
1416        return self.CheckModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz", self.LastModifiedOrg(org))
1417
1418    def LastModifiedOrg(self, org):
1419        return self.shelve.get((Update.UpdateAnnotation, (org,)), None)
1420
1421    def GetLastModified(self, addr):
1422        stream = urllib2.urlopen(addr)
1423        return datetime.strptime(stream.headers.get("Last-Modified"), "%a, %d %b %Y %H:%M:%S %Z")
1424##        return stream.headers.get("Last-Modified")
1425
1426    def GetAvailableOrganisms(self):
1427        source = urllib2.urlopen("http://www.geneontology.org/gene-associations/").read()
1428        return [s.split(".")[1] for s in sorted(set(re.findall("gene_association\.[a-zA-z0-9_]+?\.gz", source)))]
1429
1430    def GetDownloadedOrganisms(self):
1431        return [name.split(".")[1] for name in os.listdir(self.local_database_path) if name.startswith("gene_association")]
1432
1433    def IsUpdatable(self, func, args):
1434        if func == Update.UpdateOntology:
1435            return self.CheckModified("http://www.geneontology.org/ontology/gene_ontology.obo", self.shelve.get((Update.UpdateOntology, ()), None))
1436        elif func == Update.UpdateAnnotation:
1437            return self.CheckModifiedOrg(args[0])
1438
1439    def GetDownloadable(self):
1440        orgs = set(self.GetAvailableOrganisms()) - set(self.GetDownloadedOrganisms())
1441        ret = []
1442        if (Update.UpdateOntology, ()) not in self.shelve:
1443            ret.append((Update.UpdateOntology, ()))
1444        if orgs:
1445            ret.extend([(Update.UpdateAnnotation, (org,)) for org in orgs])
1446        return ret
1447
1448    def UpdateOntology(self):
1449        Ontology.DownloadOntology(os.path.join(self.local_database_path, "gene_ontology_edit.obo.tar.gz"), self.progressCallback)
1450        self._update(Update.UpdateOntology, (), self.GetLastModified("http://www.geneontology.org/ontology/gene_ontology.obo"))
1451
1452    def UpdateAnnotation(self, org):
1453        Annotations.DownloadAnnotations(org, os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"), self.progressCallback)
1454        self._update(Update.UpdateAnnotation, (org,), self.GetLastModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz"))
1455
1456    def UpdateTaxonomy(self, org):
1457        exclude = ["goa_uniprot", "goa_pdb", "GeneDB_tsetse", "reactome", "goa_zebrafish", "goa_rat", "goa_mouse"]
1458
1459        orgs = self.GetAvailableOrganisms()
1460        tax = defaultdict(set)
1461
1462        for org in orgs:
1463            if org in exclude:
1464                continue
1465            try:
1466                a = obiGO.Annotations(os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"))
1467                taxons = set(ann.taxon for ann in a.annotations)
1468                for taxId in [t.split(":")[-1] for t in taxons if "|" not in t]:  # exclude taxons with cardinality 2
1469                    tax[taxId].add(org)
1470            except Exception, ex:
1471                print ex
1472
1473        cPickle.dump(dict(tax), open(os.path.join(path, "taxonomy.pickle"), "wb"))
1474
1475
1476def _test1():
1477##    Ontology.DownloadOntology("ontology_arch.tar.gz")
1478##    Annotations.DownloadAnnotations("sgd", "annotations_arch.tar.gz")
1479    def _print(f):
1480        print f
1481    o = Ontology("ontology_arch.tar.gz")
1482    a = Annotations("annotations_arch.tar.gz", ontology=o)
1483
1484    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
1485##    profile.runctx("a.GetEnrichedTerms(sorted(a.geneNames)[:100])", {"a":a}, {})
1486    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
1487    d1 = a.GetEnrichedTerms(sorted(a.geneNames)[:1000])  # , progressCallback=_print)
1488
1489##    print a.GetEnrichedTerms(sorted(a.geneNames)[:100])#, progressCallback=_print)
1490
1491
1492def _test2():
1493    o = Ontology()
1494    a = Annotations("human", ontology=o)
1495    clusterGenes = sorted(a.geneNames)[:100]
1496    for i in range(10):
1497        terms = a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["P"])
1498        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["C"])
1499        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["F"])
1500        print i
1501#    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1502
1503##    drawEnrichmentGraph([("bal", 1.0, 5, 6, 0.1, 0.4, ["vv"]),
1504##                        ("GO:0019079", 0.5, 5, 6, 0.1, 0.4, ["cc", "bb"]),
1505##                        ("GO:0022415", 0.4, 5, 7, 0.11, 0.4, ["cc1", "bb"])], open("graph.png", "wb"), None, None)
1506
1507
1508def _test3():
1509    o = Ontology()
1510    a = Annotations("sgd", ontology=o)
1511##    a = Annotations(list(a)[3:len(a)/3], ontology=o)
1512    clusterGenes = sorted(a.geneNames)[:1] + sorted(a.geneNames)[-1:]
1513##    clusterGenes = [g + "_" + str(i%5) for g in sorted(a.geneNames)[:2]]
1514    exonMap = dict([(gene, [gene + "_E%i" % i for i in range(10)]) for gene in a.geneNames])
1515    a.RemapGenes(exonMap)
1516##    o.reverseAliasMapper = o.aliasMapper = {}
1517    terms = a.GetEnrichedTerms(exonMap.values()[0][:2] + exonMap.values()[-1][2:])
1518##    terms = a.GetEnrichedTerms(clusterGenes)
1519    print terms
1520##    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1521    a.DrawEnrichmentGraph(filterByPValue(terms, maxPValue=0.1), len(clusterGenes), len(a.geneNames))
1522
1523if __name__ == "__main__":
1524    _test2()
Note: See TracBrowser for help on using the repository browser.