source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1768:9cb38945684f

Revision 1768:9cb38945684f, 59.3 KB checked in by Ales Erjavec <ales.erjavec@…>, 12 months ago (diff)

Removed deprecated and nonfunctioning 'drawEnrichmentGraph_tostream' function.

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        """ Return object with id (same as ontology.terms[id]
447        """
448        return self.terms.get(id, self.terms.get(self.aliasMapper.get(id, id)))
449
450    def __iter__(self):
451        """ Iterate over all ids in ontology
452        """
453        return iter(self.terms)
454
455    def __len__(self):
456        """ Return number of objects in ontology
457        """
458        return len(self.terms)
459
460    def __contains__(self, id):
461        return id in self.terms or id in self.aliasMapper
462
463    @staticmethod
464    def DownloadOntology(file, progressCallback=None):
465        tFile = tarfile.open(file, "w:gz") if isinstance(file, basestring) \
466                else file
467        tmpDir = os.path.join(orngEnviron.bufferDir, "tmp_go/")
468        try:
469            os.mkdir(tmpDir)
470        except Exception:
471            pass
472        urlretrieve("http://www.geneontology.org/ontology/gene_ontology_edit.obo",
473                    os.path.join(tmpDir, "gene_ontology_edit.obo"),
474                    progressCallback and __progressCallbackWrapper(progressCallback))
475        tFile.add(os.path.join(tmpDir, "gene_ontology_edit.obo"),
476                  "gene_ontology_edit.obo")
477        tFile.close()
478        os.remove(os.path.join(tmpDir, "gene_ontology_edit.obo"))
479
480    @staticmethod
481    def DownloadOntologyAtRev(rev, filename=None, progressCallback=None):
482        url = "http://cvsweb.geneontology.org/cgi-bin/cvsweb.cgi/~checkout~/go/ontology/gene_ontology_edit.obo?rev=%s" % rev
483        url += ";content-type=text%2Fplain"
484        if filename is None:
485            filename = os.path.join(default_database_path,
486                                    "gene_ontology_edit@rev%s.obo" % rev)
487        r = urllib2.urlopen(url)
488
489        with open(filename + ".part", "wb") as f:
490            shutil.copyfileobj(r, f)
491
492        os.rename(filename + ".part", filename)
493
494
495_re_obj_name_ = re.compile("([a-zA-z0-9-_]+)")
496
497
498class AnnotationRecord(object):
499    """Holds the data for an annotation record read from the annotation file.
500    Fields can be accessed with the names: DB, DB_Object_ID, DB_Object_Symbol,
501    Qualifier, GO_ID, DB_Reference, Evidence_code, With_or_From, Aspect,
502    DB_Object_Name, DB_Object_Synonym, DB_Object_Type, taxon, Date,
503    Assigned_by (e.g. rec.GO_ID) or by supplying the original name of the
504    field (see http://geneontology.org/GO.annotation.shtml#file) to the get
505    method (e.g. rec.get("GO ID")) The object also provides the following
506    data members for quicker access: geneName, GOId, evidence, aspect and
507    alias(a list of aliases)
508
509    """
510    __slots__ = annotationFields + ["geneName", "GOId", "evidence",
511                                    "aspect", "alias", "additionalAliases"]
512
513    def __init__(self, fullText):
514        """\
515        :param fulText: A single line from the annotation file.
516
517        """
518        for slot, val in zip(annotationFields, fullText.split("\t")):
519            setattr(self, slot, intern(val))
520        self.geneName = self.DB_Object_Symbol
521        self.GOId = self.GO_ID
522        self.evidence = self.Evidence_Code
523        self.aspect = self.Aspect
524        self.alias = list(map(intern, self.DB_Object_Synonym.split("|")))
525
526        self.additionalAliases = []
527        if ":" in self.DB_Object_Name:
528            self.additionalAliases = []  # _re_obj_name_.findall(self.DB_Object_Name.split(":")[0])
529
530    def __getattr__(self, name):
531        if name in annotationFieldsDict:
532            return getattr(self, self.__slots__[annotationFieldsDict[name]])
533        else:
534            raise AttributeError(name)
535
536
537class Annotations(object):
538    """Annotations object holds the annotations.
539    """
540    version = 2
541
542    def __init__(self, file=None, ontology=None, genematcher=None,
543                 progressCallback=None, rev=None):
544        """Initialize the annotations from file by calling ParseFile on it.
545        The ontology must be an instance of Ontology class.
546        The optional progressCallback will be called with a single argument
547        to report on the progress.
548
549        """
550        self.file = file
551        self.ontology = ontology
552        self.allAnnotations = defaultdict(list)
553        self.geneAnnotations = defaultdict(list)
554        self.termAnnotations = defaultdict(list)
555        self._geneNames = None
556        self._geneNamesDict = None
557        self._aliasMapper = None
558        self.additionalAliases = {}
559        self.annotations = []
560        self.header = ""
561        self.genematcher = genematcher
562        self.taxid = None
563        if type(file) in [list, set, dict, Annotations]:
564            for ann in file:
565                self.AddAnnotation(ann)
566            if type(file, Annotations):
567                taxid = file.taxid
568        elif isinstance(file, basestring) and os.path.exists(file):
569            self.ParseFile(file, progressCallback)
570            try:
571                self.taxid = to_taxid(os.path.basename(file).split(".")[1]).pop()
572            except IOError:
573                pass
574        elif file is not None:
575            # Organism code
576            if rev is not None:
577                if not _CVS_REVISION_RE.match(rev):
578                    raise ValueError("Invalid revision format")
579
580                if rev.startswith("rev"):
581                    rev = rev[3:]
582                code = self.organism_name_search(file)
583                filename = os.path.join(default_database_path,
584                                        "gene_association.%s@rev%s.tar.gz" %
585                                        (code, rev))
586
587                if not os.path.exists(filename):
588                    self.DownloadAnnotationsAtRev(code, rev, filename,
589                                                  progressCallback
590                                                  )
591
592                self.ParseFile(filename, progressCallback)
593                self.taxid = to_taxid(code).pop()
594            else:
595                a = self.Load(file, ontology, genematcher, progressCallback)
596                self.__dict__ = a.__dict__
597                self.taxid = to_taxid(organism_name_search(file)).pop()
598        if not self.genematcher and self.taxid:
599            self.genematcher = obiGene.matcher(
600                [obiGene.GMGO(self.taxid)] +
601                ([obiGene.GMDicty(),
602                  [obiGene.GMGO(self.taxid), obiGene.GMDicty()]]
603                 if self.taxid == "352472" else [])
604            )
605        if self.genematcher:
606            self.genematcher.set_targets(self.geneNames)
607
608    @classmethod
609    def organism_name_search(cls, org):
610        ids = to_taxid(org)
611        from . import obiTaxonomy as tax
612        if not ids:
613            ids = [org] if org in Taxonomy().common_org_map.keys() + \
614                  Taxonomy().code_map.keys() else []
615        if not ids:
616            ids = tax.to_taxid(org, mapTo=Taxonomy().keys())
617        if not ids:
618            ids = tax.search(org, exact=True)
619            ids = set(ids).intersection(Taxonomy().keys())
620        if not ids:
621            ids = tax.search(org)
622            ids = set(ids).intersection(Taxonomy().keys())
623        codes = set([from_taxid(id) for id in ids])
624        if len(codes) > 1:
625            raise tax.MultipleSpeciesException(
626                ", ".join(["%s: %s" % (str(from_taxid(id)), tax.name(id))
627                           for id in ids]))
628        elif len(codes) == 0:
629            raise tax.UnknownSpeciesIdentifier(org)
630        return codes.pop()
631
632    @classmethod
633    def organism_version(cls, name):
634        name = organism_name_search(name)
635        orngServerFiles.localpath_download(
636            "GO", "gene_association.%s.tar.gz" % name)
637        return ("v%i." % cls.version) + orngServerFiles.info("GO",
638                        "gene_association.%s.tar.gz" % name)["datetime"]
639
640    def SetOntology(self, ontology):
641        """ Set the ontology to use in the annotations mapping.
642        """
643        self.allAnnotations = defaultdict(list)
644        self._ontology = ontology
645
646    def GetOntology(self):
647        return self._ontology
648
649    ontology = property(GetOntology, SetOntology,
650                        doc="Ontology object for annotations")
651
652    @classmethod
653    def Load(cls, org, ontology=None, genematcher=None, progressCallback=None):
654        """A class method that tries to load the association file for the
655        given organism from default_database_path.
656        """
657        code = organism_name_search(org)
658
659        file = "gene_association.%s.tar.gz" % code
660
661        path = os.path.join(orngServerFiles.localpath("GO"), file)
662
663        if not os.path.exists(path):
664            sf = orngServerFiles.ServerFiles()
665            available = sf.listfiles("GO")
666            if file not in available:
667                from . import obiKEGG
668                raise obiKEGG.OrganismNotFoundError(org + str(code))
669            orngServerFiles.download("GO", file)
670
671        return cls(path, ontology=ontology, genematcher=genematcher,
672                   progressCallback=progressCallback)
673
674    def ParseFile(self, file, progressCallback=None):
675        """ Parse and load the annotations from file. Report progress
676        with progressCallback.
677        File can be:
678            - a tarball containing the association file named gene_association
679            - a directory name containing the association file named gene_association
680            - a path to the actual association file
681            - an open file-like object of the association file
682        """
683        if type(file) == str:
684            if os.path.isfile(file) and tarfile.is_tarfile(file):
685                f = tarfile.open(file).extractfile("gene_association")
686            elif os.path.isfile(file) and file.endswith(".gz"):
687                f = gzip.open(file)
688            elif os.path.isfile(file):
689                f = open(file)
690            else:
691                f = open(os.path.join(file, "gene_association"))
692        else:
693            f = file
694        lines = [line for line in f.read().splitlines() if line.strip()]
695
696        milestones = orngMisc.progressBarMilestones(len(lines), 100)
697        for i, line in enumerate(lines):
698            if line.startswith("!"):
699                self.header = self.header + line + "\n"
700                continue
701
702            a = AnnotationRecord(line)
703            self.AddAnnotation(a)
704#            self.annotations.append(a)
705            if progressCallback and i in milestones:
706                progressCallback(100.0 * i / len(lines))
707
708    def AddAnnotation(self, a):
709        """ Add a single `AnotationRecord` instance to this `Annotations`
710        object.
711        """
712        if not isinstance(a, AnnotationRecord):
713            a = AnnotationRecord(a)
714        if not a.geneName or not a.GOId or a.Qualifier == "NOT":
715            return
716
717        self.geneAnnotations[a.geneName].append(a)
718        self.annotations.append(a)
719        self.termAnnotations[a.GOId].append(a)
720        self.allAnnotations = defaultdict(list)
721
722        self._geneNames = None
723        self._geneNamesDict = None
724        self._aliasMapper = None
725
726    @property
727    def geneNamesDict(self):
728        if getattr(self, "_geneNamesDict", None) is None:
729            self._geneNamesDict = defaultdict(set)
730            for alias, name in self.aliasMapper.iteritems():
731                self._geneNamesDict[name].add(alias)
732        return self._geneNamesDict
733
734    @property
735    def geneNames(self):
736        if getattr(self, "_geneNames", None) is None:
737            self._geneNames = set([ann.geneName for ann in self.annotations])
738        return self._geneNames
739
740    @property
741    def aliasMapper(self):
742        if getattr(self, "_aliasMapper", None) is None:
743            self._aliasMapper = {}
744            for ann in self.annotations:
745                self._aliasMapper.update([(alias, ann.geneName)
746                                          for alias in ann.alias +
747                                           [ann.geneName, ann.DB_Object_ID]])
748        return self._aliasMapper
749
750    def GetGeneNamesTranslator(self, genes):
751        """ Return a dictionary mapping canonical names (DB_Object_Symbol)
752        to `genes`.
753
754        """
755        def alias(gene):
756            if self.genematcher:
757                return self.genematcher.umatch(gene)
758            else:
759                return gene if gene in self.geneNames else \
760                        self.aliasMapper.get(gene,
761                             self.additionalAliases.get(gene, None))
762
763        return dict([(alias(gene), gene) for gene in genes if alias(gene)])
764
765    def _CollectAnnotations(self, id, visited):
766        """ Recursive function collects and caches all annotations for id
767        """
768        if id not in self.allAnnotations and id not in visited:
769            if id in self.ontology.reverseAliasMapper:
770                annotations = [self.termAnnotations.get(alt_id, [])
771                               for alt_id in
772                               self.ontology.reverseAliasMapper[id]] + \
773                              [self.termAnnotations[id]]
774            else:
775                ## annotations for this term alone
776                annotations = [self.termAnnotations[id]]
777            visited.add(id)
778            for typeId, child in self.ontology[id].relatedTo:
779                aa = self._CollectAnnotations(child, visited)
780                if type(aa) == set:
781                    ## if it was already reduced in GetAllAnnotations
782                    annotations.append(aa)
783                else:
784                    annotations.extend(aa)
785            self.allAnnotations[id] = annotations
786        return self.allAnnotations[id]
787
788    def GetAllAnnotations(self, id):
789        """ Return a set of all annotations (instances if `AnnotationRectord`)
790        for GO term `id` and all it's subterms.
791
792        :param id: GO term id
793        :type id: str
794
795        """
796        visited = set()
797        id = self.ontology.aliasMapper.get(id, id)
798        if id not in self.allAnnotations or \
799                type(self.allAnnotations[id]) == list:
800            annot_set = set()
801            for annots in self._CollectAnnotations(id, set()):
802                annot_set.update(annots)
803            self.allAnnotations[id] = annot_set
804        return self.allAnnotations[id]
805
806    def GetAllGenes(self, id, evidenceCodes=None):
807        """ Return a list of genes annotated by specified `evidenceCodes`
808        to GO term 'id' and all it's subterms."
809
810        :param id: GO term id
811        :type id: str
812
813        :param evidneceCodes: List of evidence codes to consider when
814                              matching annotations to terms.
815        :type evidenceCodes: list-of-strings
816        """
817        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
818        annotations = self.GetAllAnnotations(id)
819        return list(set([ann.geneName for ann in annotations
820                         if ann.Evidence_Code in evidenceCodes]))
821
822    def GetEnrichedTerms(self, genes, reference=None, evidenceCodes=None,
823                         slimsOnly=False, aspect=None, prob=obiProb.Binomial(),
824                         useFDR=True, progressCallback=None):
825        """ Return a dictionary of enriched terms, with tuples of
826        (list_of_genes, p_value, reference_count) for items and term
827        ids as keys. P-Values are FDR adjusted if useFDR is True (default).
828
829        :param genes: List of genes
830        :param reference: list of genes (if None all genes included in the
831                          annotations will be used).
832        :param evidenceCodes: List of evidence codes to consider.
833        :param slimsOnly: If `True` return only slim terms
834        :param aspect: Which aspects to use. Use all by default. "P", "F", "C"
835            or a set containing these elements.
836        """
837        revGenesDict = self.GetGeneNamesTranslator(genes)
838        genes = set(revGenesDict.keys())
839        if reference:
840            refGenesDict = self.GetGeneNamesTranslator(reference)
841            reference = set(refGenesDict.keys())
842        else:
843            reference = self.geneNames
844
845        if aspect == None:
846            aspects_set = set(["P", "C", "F"])
847        elif isinstance(aspect, basestring):
848            aspects_set = set([aspect])
849        else:
850            aspects_set = aspect
851
852        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
853        annotations = [ann
854                       for gene in genes for ann in self.geneAnnotations[gene]
855                       if ann.Evidence_Code in evidenceCodes and
856                       ann.Aspect in aspects_set]
857
858        refAnnotations = set(
859            [ann
860             for gene in reference for ann in self.geneAnnotations[gene]
861             if ann.Evidence_Code in evidenceCodes and
862             ann.Aspect in aspects_set]
863        )
864
865        annotationsDict = defaultdict(set)
866        for ann in annotations:
867            annotationsDict[ann.GO_ID].add(ann)
868
869        if slimsOnly and not self.ontology.slimsSubset:
870            warnings.warn("Unspecified slims subset in the ontology! "
871                          "Using 'goslim_generic' subset", UserWarning)
872            self.ontology.SetSlimsSubset("goslim_generic")
873
874        terms = self.ontology.ExtractSuperGraph(annotationsDict.keys())
875        res = {}
876
877        milestones = orngMisc.progressBarMilestones(len(terms), 100)
878        for i, term in enumerate(terms):
879            if slimsOnly and term not in self.ontology.slimsSubset:
880                continue
881            allAnnotations = self.GetAllAnnotations(term).intersection(refAnnotations)
882##            allAnnotations.intersection_update(refAnnotations)
883            allAnnotatedGenes = set([ann.geneName for ann in allAnnotations])
884            mappedGenes = genes.intersection(allAnnotatedGenes)
885##            if not mappedGenes:
886##                print >> sys.stderr, term, sorted(genes)
887##                print >> sys.stderr, sorted(allAnnotatedGenes)
888##                return
889            if len(reference) > len(allAnnotatedGenes):
890                mappedReferenceGenes = reference.intersection(allAnnotatedGenes)
891            else:
892                mappedReferenceGenes = allAnnotatedGenes.intersection(reference)
893            res[term] = ([revGenesDict[g] for g in mappedGenes],
894                         prob.p_value(len(mappedGenes), len(reference),
895                                      len(mappedReferenceGenes), len(genes)),
896                         len(mappedReferenceGenes))
897            if progressCallback and i in milestones:
898                progressCallback(100.0 * i / len(terms))
899        if useFDR:
900            res = sorted(res.items(), key=lambda (_1, (_2, p, _3)): p)
901            res = dict([(id, (genes, p, ref))
902                        for (id, (genes, _, ref)), p in
903                        zip(res, obiProb.FDR([p for _, (_, p, _) in res]))])
904        return res
905
906    def GetAnnotatedTerms(self, genes, directAnnotationOnly=False,
907                          evidenceCodes=None, progressCallback=None):
908        """ Return all terms that are annotated by genes with evidenceCodes.
909        """
910        genes = [genes] if type(genes) == str else genes
911        revGenesDict = self.GetGeneNamesTranslator(genes)
912        genes = set(revGenesDict.keys())
913        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
914        annotations = [ann for gene in genes for ann in self.geneAnnotations[gene]
915                       if ann.Evidence_Code in evidenceCodes]
916        dd = defaultdict(set)
917        for ann in annotations:
918            dd[ann.GO_ID].add(revGenesDict.get(ann.geneName, ann.geneName))
919        if not directAnnotationOnly:
920            terms = self.ontology.ExtractSuperGraph(dd.keys())
921            for i, term in enumerate(terms):
922                termAnnots = self.GetAllAnnotations(term).intersection(annotations)
923##                termAnnots.intersection_update(annotations)
924                dd[term].update([revGenesDict.get(ann.geneName, ann.geneName)
925                                 for ann in termAnnots])
926        return dict(dd)
927
928    def DrawEnrichmentGraph(self, terms, clusterSize, refSize=None,
929                            file="graph.png", width=None, height=None,
930                            precison=3):
931        refSize = len(self.geneNames) if refSize == None else refSize
932        sortedterms = sorted(terms.items(), key=lambda term: term[1][1])
933        fdr = dict(zip([t[0] for t in sortedterms],
934                       obiProb.FDR([t[1][1] for t in sortedterms])))
935        termsList = [(term,
936                      ((float(len(terms[term][0])) / clusterSize) /
937                       (float(terms[term][2]) / refSize)),
938                      len(terms[term][0]),
939                      terms[term][2],
940                      terms[term][1],
941                      fdr[term],
942                      terms[term][0])
943                     for term in terms]
944
945        drawEnrichmentGraph(termsList, file, width, height,
946                            ontology=self.ontology, precison=precison)
947
948    def __add__(self, iterable):
949        """ Return a new Annotations object with combined annotations
950        """
951        return Annotations([a for a in self] + [a for a in iterable],
952                           ontology=self.ontology)
953
954    def __iadd__(self, iterable):
955        """ Add annotations to this instance
956        """
957        self.extend(iterable)
958        return self
959
960    def __contains__(self, item):
961        return item in self.annotations
962
963    def __iter__(self):
964        """ Iterate over all AnnotationRecord objects in annotations
965        """
966        return iter(self.annotations)
967
968    def __len__(self):
969        """ Return the number of annotations
970        """
971        return len(self.annotations)
972
973    def __getitem__(self, index):
974        """ Return the i-th annotation record
975        """
976        return self.annotations[index]
977
978    def __getslice__(self, *args):
979        return self.annotations.__getslice__(*args)
980
981    def add(self, line):
982        """ Add one annotation
983        """
984        self.AddAnnotation(line)
985
986    def append(self, line):
987        """ Add one annotation
988        """
989        self.AddAnnotation(line)
990
991    def extend(self, lines):
992        """ Add multiple annotations
993        """
994        for line in lines:
995            self.AddAnnotation(line)
996
997    def RemapGenes(self, map):
998        """
999        """
1000        for gene in map:
1001            annotations = self.geneAnnotations[gene]
1002            for ann in annotations:
1003                for name in map[gene]:
1004                    ann1 = copy.copy(ann)
1005                    ann1.geneName = name
1006                    self.add(ann1)
1007        self.genematcher = obiGene.GMDirect()
1008        try:
1009            del self._geneNames
1010        except Exception:
1011            pass
1012        self.genematcher.set_targets(self.geneNames)
1013
1014    @staticmethod
1015    def DownloadAnnotations(org, file, progressCallback=None):
1016        tFile = tarfile.open(file, "w:gz") if type(file) == str else file
1017        tmpDir = os.path.join(orngEnviron.bufferDir, "tmp_go/")
1018        try:
1019            os.mkdir(tmpDir)
1020        except Exception:
1021            pass
1022        fileName = "gene_association." + org + ".gz"
1023        urlretrieve(("http://www.geneontology.org/gene-associations/" +
1024                     fileName),
1025                    os.path.join(tmpDir, fileName),
1026                    progressCallback and __progressCallbackWrapper(progressCallback))
1027        gzFile = GzipFile(os.path.join(tmpDir, fileName), "r")
1028        file = open(os.path.join(tmpDir, "gene_association." + org), "w")
1029        file.writelines(gzFile.readlines())
1030        file.flush()
1031        file.close()
1032
1033        tFile.add(os.path.join(tmpDir, "gene_association." + org),
1034                  "gene_association")
1035        annotation = Annotations(os.path.join(tmpDir, "gene_association." + org),
1036                    genematcher=obiGene.GMDirect(), progressCallback=progressCallback)
1037        cPickle.dump(annotation.geneNames, open(os.path.join(tmpDir, "gene_names.pickle"), "wb"))
1038        tFile.add(os.path.join(tmpDir, "gene_names.pickle"), "gene_names.pickle")
1039        tFile.close()
1040        os.remove(os.path.join(tmpDir, "gene_association." + org))
1041        os.remove(os.path.join(tmpDir, "gene_names.pickle"))
1042
1043    @staticmethod
1044    def DownloadAnnotationsAtRev(org, rev, filename=None,
1045                                 progressCallback=None):
1046        if filename is None:
1047            filename = os.path.join(default_database_path,
1048                                    "gene_association.%s@rev%s.tar.gz" %
1049                                    (code, rev))
1050        url = ("http://cvsweb.geneontology.org/cgi-bin/cvsweb.cgi/~checkout~/go/gene-associations/gene_association.%s.gz?rev=%s" %
1051               (org, rev))
1052        url += ";content-type=application%2Fx-gzip"
1053        r = urllib2.urlopen(url)
1054
1055        with open(filename + ".part", "wb") as f:
1056            shutil.copyfileobj(r, f)
1057
1058        os.rename(filename + ".part", filename)
1059
1060from .obiTaxonomy import pickled_cache
1061
1062
1063@pickled_cache(None, [("GO", "taxonomy.pickle"),
1064                      ("Taxonomy", "ncbi_taxonomy.tar.gz")])
1065def organism_name_search(name):
1066    return Annotations.organism_name_search(name)
1067
1068
1069def filterByPValue(terms, maxPValue=0.1):
1070    """ Filters the terms by the p-value. Asumes terms is a dict with
1071    the same structure as returned from GetEnrichedTerms.
1072
1073    """
1074    return dict(filter(lambda (k, e): e[1] <= maxPValue, terms.items()))
1075
1076
1077def filterByFrequency(terms, minF=2):
1078    """ Filters the terms by the cluster frequency. Asumes terms is
1079    a dict with the same structure as returned from GetEnrichedTerms.
1080
1081    """
1082    return dict(filter(lambda (k, e): len(e[0]) >= minF, terms.items()))
1083
1084
1085def filterByRefFrequency(terms, minF=4):
1086    """ Filters the terms by the reference frequency. Asumes terms is
1087    a dict with the same structure as returned from GetEnrichedTerms.
1088
1089    """
1090    return dict(filter(lambda (k, e): e[2] >= minF, terms.items()))
1091
1092
1093def drawEnrichmentGraph(enriched, file="graph.png", width=None, height=None,
1094                        header=None, ontology=None, precison=3):
1095    file = open(file, "wb") if type(file) == str else file
1096    drawEnrichmentGraph_tostreamMk2(enriched, file, width, height, header,
1097                                    ontology, precison)
1098
1099
1100def drawEnrichmentGraph_tostreamMk2(enriched, fh, width, height, header=None,
1101                                    ontology=None, precison=4):
1102    ontology = ontology if ontology else Ontology()
1103    header = header if header else ["List", "Total", "p-value", "FDR",
1104                                    "Names", "Genes"]
1105    GOTerms = dict([(t[0], t) for t in enriched if t[0] in ontology])
1106
1107    def getParents(term):
1108        parents = ontology.ExtractSuperGraph([term])
1109        parents = [id for id in parents if id in GOTerms and id != term]
1110        c = reduce(set.union, [set(ontology.ExtractSuperGraph([id])) -
1111                               set([id]) for id in parents], set())
1112        parents = [t for t in parents if t not in c]
1113        return parents
1114    parents = dict([(term, getParents(term)) for term in GOTerms])
1115    # print "Parentes", parents
1116
1117    def getChildren(term):
1118        return [id for id in GOTerms if term in parents[id]]
1119    topLevelTerms = [id for id in parents if not parents[id]]
1120    # print "Top level terms", topLevelTerms
1121    termsList = []
1122    fmt = "%" + ".%if" % precison
1123
1124    def collect(term, parent):
1125        termsList.append(GOTerms[term][1:4] + \
1126                         (fmt % GOTerms[term][4],
1127                          fmt % GOTerms[term][5],
1128                          ontology[term].name,
1129                          ", ".join(GOTerms[term][6])) + (parent,))
1130        parent = len(termsList) - 1
1131        for c in getChildren(term):
1132            collect(c, parent)
1133
1134    for topTerm in topLevelTerms:
1135        collect(topTerm, None)
1136    for entry in enriched:
1137        if entry[0] not in ontology:
1138            termsList.append(entry[1:4] + \
1139                             (fmt % entry[4],
1140                              fmt % entry[5],
1141                              entry[0],
1142                              ", ".join(entry[6])) + (None,))
1143
1144    drawEnrichmentGraphPIL_tostream(termsList, header, fh, width, height)
1145
1146
1147def drawEnrichmentGraphPIL_tostream(termsList, headers, fh, width=None, height=None):
1148    from PIL import Image, ImageDraw, ImageFont
1149    backgroundColor = (255, 255, 255)
1150    textColor = (0, 0, 0)
1151    graphColor = (0, 0, 255)
1152    fontSize = height == None and 12 or (height - 60) / len(termsList)
1153    font = ImageFont.load_default()
1154    try:
1155        font = ImageFont.truetype("arial.ttf", fontSize)
1156    except:
1157        pass
1158    getMaxTextHeightHint = lambda l: max([font.getsize(t)[1] for t in l])
1159    getMaxTextWidthHint = lambda l: max([font.getsize(t)[0] for t in l])
1160    maxFoldWidth = width != None and min(150, width / 6) or 150
1161    maxFoldEnrichment = max([t[0] for t in termsList])
1162    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
1163    foldWidths = [int(foldNormalizationFactor * term[0]) for term in termsList]
1164    treeStep = 10
1165    treeWidth = {}
1166    for i, term in enumerate(termsList):
1167        treeWidth[i] = (term[-1] == None and 1 or treeWidth[term[-1]] + 1)
1168    treeStep = width != None and min(treeStep, width / (6 * max(treeWidth.values())) or 2) or treeStep
1169    treeWidth = [w * treeStep + foldWidths[i] for i, w in treeWidth.items()]
1170    treeWidth = max(treeWidth) - maxFoldWidth
1171    verticalMargin = 10
1172    horizontalMargin = 10
1173##    print verticalMargin, maxFoldWidth, treeWidth
1174##    treeWidth = 100
1175    firstColumnStart = verticalMargin + maxFoldWidth + treeWidth + 10
1176    secondColumnStart = firstColumnStart + getMaxTextWidthHint([str(t[1]) for t in termsList] + [headers[0]]) + 2
1177    thirdColumnStart = secondColumnStart + getMaxTextWidthHint([str(t[2]) for t in termsList] + [headers[1]]) + 2
1178    fourthColumnStart = thirdColumnStart + getMaxTextWidthHint([str(t[3]) for t in termsList] + [headers[2]]) + 2
1179    fifthColumnStart = fourthColumnStart + getMaxTextWidthHint([str(t[4]) for t in termsList] + [headers[3]]) + 4
1180##    maxAnnotationTextWidth = width==None and getMaxTextWidthHint([str(t[4]) for t in termsList]+["Annotation"]) or (width - fourthColumnStart - verticalMargin) * 2 / 3
1181    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]]))
1182    sixthColumnStart = fifthColumnStart + maxAnnotationTextWidth + 4
1183    maxGenesTextWidth = width == None and getMaxTextWidthHint([str(t[6]) for t in termsList] + [headers[5]]) or (width - fifthColumnStart - verticalMargin) / 3
1184
1185    legendHeight = font.getsize("1234567890")[1] * 2
1186    termHeight = font.getsize("A")[1]
1187##    print fourthColumnStart, maxAnnotationTextWidth, verticalMargin
1188    width = sixthColumnStart + maxGenesTextWidth + verticalMargin
1189    height = len(termsList) * termHeight + 2 * (legendHeight + horizontalMargin)
1190
1191    image = Image.new("RGB", (width, height), backgroundColor)
1192    draw = ImageDraw.Draw(image)
1193
1194    def truncText(text, maxWidth, append=""):
1195        # print getMaxTextWidthHint([text]), maxAnnotationTextWidth
1196        if getMaxTextWidthHint([text]) > maxWidth:
1197            while getMaxTextWidthHint([text + "..." + append]) > maxWidth and text:
1198                text = text[:-1]
1199            if text:
1200                text = text + "..." + append
1201            else:
1202                text = append
1203        return text
1204    currentY = horizontalMargin + legendHeight
1205    connectAtX = {}
1206    for i, term in enumerate(termsList):
1207        draw.line([(verticalMargin, currentY + termHeight / 2), (verticalMargin + foldWidths[i], currentY + termHeight / 2)], width=termHeight - 2, fill=graphColor)
1208        draw.text((firstColumnStart, currentY), str(term[1]), font=font, fill=textColor)
1209        draw.text((secondColumnStart, currentY), str(term[2]), font=font, fill=textColor)
1210        draw.text((thirdColumnStart, currentY), str(term[3]), font=font, fill=textColor)
1211        draw.text((fourthColumnStart, currentY), str(term[4]), font=font, fill=textColor)
1212##        annotText = width!=None and truncText(str(term[5]), maxAnnotationTextWidth, str(term[5])) or str(term[4])
1213        annotText = width != None and truncText(str(term[5]), maxAnnotationTextWidth)
1214        draw.text((fifthColumnStart, currentY), annotText, font=font, fill=textColor)
1215        genesText = width != None and truncText(str(term[6]), maxGenesTextWidth) or str(term[6])
1216        draw.text((sixthColumnStart, currentY), genesText, font=font, fill=textColor)
1217        lineEnd = term[-1] == None and firstColumnStart - 10 or connectAtX[term[-1]]
1218        draw.line([(verticalMargin + foldWidths[i] + 1, currentY + termHeight / 2), (lineEnd, currentY + termHeight / 2)], width=1, fill=textColor)
1219        if term[-1] != None:
1220            draw.line([(lineEnd, currentY + termHeight / 2), (lineEnd, currentY + termHeight / 2 - termHeight * (i - term[-1]))], width=1, fill=textColor)
1221        connectAtX[i] = lineEnd - treeStep
1222        currentY += termHeight
1223
1224    currentY = horizontalMargin
1225    draw.text((firstColumnStart, currentY), headers[0], font=font, fill=textColor)
1226    draw.text((secondColumnStart, currentY), headers[1], font=font, fill=textColor)
1227    draw.text((thirdColumnStart, currentY), headers[2], font=font, fill=textColor)
1228    draw.text((fourthColumnStart, currentY), headers[3], font=font, fill=textColor)
1229    draw.text((fifthColumnStart, currentY), headers[4], font=font, fill=textColor)
1230    draw.text((sixthColumnStart, currentY), headers[5], font=font, fill=textColor)
1231
1232    horizontalMargin = 0
1233    # draw.line([(verticalMargin, height - horizontalMargin - legendHeight), (verticalMargin + maxFoldWidth, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1234    draw.line([(verticalMargin, horizontalMargin + legendHeight), (verticalMargin + maxFoldWidth, horizontalMargin + legendHeight)], width=1, fill=textColor)
1235    maxLabelWidth = getMaxTextWidthHint([" " + str(i) for i in range(int(maxFoldEnrichment + 1))])
1236    numOfLegendLabels = max(int(maxFoldWidth / maxLabelWidth), 2)
1237    for i in range(numOfLegendLabels + 1):
1238        # draw.line([(verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight/2), (verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1239        # draw.text((verticalMargin + i*maxFoldWidth/10 - font.getsize(str(i))[0]/2, height - horizontalMargin - legendHeight/2), str(i), font=font, fill=textColor)
1240
1241        label = str(int(i * maxFoldEnrichment / numOfLegendLabels))
1242        draw.line([(verticalMargin + i * maxFoldWidth / numOfLegendLabels, horizontalMargin + legendHeight / 2), (verticalMargin + i * maxFoldWidth / numOfLegendLabels, horizontalMargin + legendHeight)], width=1, fill=textColor)
1243        draw.text((verticalMargin + i * maxFoldWidth / numOfLegendLabels - font.getsize(label)[0] / 2, horizontalMargin), label, font=font, fill=textColor)
1244
1245    image.save(fh)
1246
1247
1248def drawEnrichmentGraphPylab_tostream(termsList, headers, fh, width=None, height=None, show=True):
1249    from matplotlib import pyplot as plt
1250    from matplotlib.patches import Rectangle
1251
1252    maxFoldWidth = width != None and min(150, width / 6) or 150
1253    maxFoldEnrichment = max([t[0] for t in termsList])
1254    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
1255##    foldWidths = [int(foldNormalizationFactor*term[0]) for term in termsList]
1256    foldWidths = [term[0] for term in termsList]
1257    treeStep = maxFoldEnrichment * 0.05
1258    treeWidth = {}
1259
1260    for i, term in enumerate(termsList):
1261        treeWidth[i] = (term[-1] == None and treeStep or treeWidth[term[-1]] + treeStep)
1262    maxTreeWidth = max(treeWidth)
1263
1264    connectAt = {}
1265    cellText = []
1266    axes1 = plt.axes([0.1, 0.1, 0.2, 0.8])
1267    for i, line in enumerate(termsList):
1268        enrichment, n, m, p_val, fdr_val, name, genes, parent = line
1269        r = Rectangle((0, len(termsList) - i - 0.4), enrichment, 0.8)
1270        plt.gca().add_patch(r)
1271        plt.plot([enrichment, connectAt.get(parent, maxFoldEnrichment + maxTreeWidth)], [len(termsList) - i, len(termsList) - i], color="black")
1272        connectAt[i] = connectAt.get(parent, maxFoldEnrichment + maxTreeWidth) - treeStep
1273        if parent != None:
1274            plt.plot([connectAt.get(parent)] * 2, [len(termsList) - i, len(termsList) - parent], color="black")
1275        cellText.append((str(n), str(m), p_val, fdr_val, name, genes))
1276
1277##    from Orange.orng.orngClustering import TableTextLayout
1278##    text = TableTextLayout((maxFoldEnrichment*1.1, len(termsList)), cellText)
1279    from Orange.orng.orngClustering import TablePlot
1280    if True:
1281        axes2 = plt.axes([0.3, 0.1, 0.6, 0.8], sharey=axes1)
1282        axes2.set_axis_off()
1283        table = TablePlot((0, len(termsList)), axes=plt.gca())
1284        for i, line in enumerate(cellText):
1285            for j, text in enumerate(line):
1286                table.add_cell(i, j, width=len(text), height=1, text=text, loc="left", edgecolor="w", facecolor="w")
1287
1288        table.set_figure(plt.gcf())
1289        plt.gca().add_artist(table)
1290        plt.gca()._set_artist_props(table)
1291##    plt.text(3, 3, "\n".join(["\t".join(text) for text in cellText]))
1292
1293##    table = plt.table(cellText=cellText, colLabels=headers, loc="right")
1294##    table.set_transform(plt.gca().transData)
1295##
1296##    table.set_xy(20,20)
1297    plt.show()
1298
1299
1300class Taxonomy(object):
1301    """Maps NCBI taxonomy ids to corresponding GO organism codes
1302    """
1303    common_org_map = {"297284": "9913", "30523": "9913",  # Bos taurus
1304                      "5782": "352472", "44689": "352472", "366501": "352472",  # Dictyostelium discoideum
1305                      "83333": "562",  # Escherichia coli
1306                      "52545": "4530", "4532": "4530", "65489": "4530", "4533": "4530", "77588": "4530", "29689": "4530",
1307                      "4538": "4530", "40148": "4530", "29690": "4530", "110450": "4530", "4534": "4530", "83309": "4530",
1308                      "4528": "4530", "127571": "4530", "40149": "4530", "83307": "4530", "63629": "4530", "4536": "4530",
1309                      "4535": "4530", "4537": "4530", "65491": "4530", "83308": "4530", "4529": "4530", "4530": "4530",
1310                      "39946": "4530", "39947": "4530", "110451": "4530", "364100": "4530", "364099": "4530", "4539": "4530",
1311                      }
1312    code_map = {"3702": "tair",  # Arabidopsis thaliana
1313                "9913": "goa_cow",  # Bos taurus
1314                "6239": "wb",  # Caenorhabditis elegans
1315                "3055": None,  # Chlamydomonas reinhardtii
1316                "7955": "zfin",  # Danio rerio (zebrafish)
1317                "352472": "dictyBase",  # Dictyostelium discoideum
1318                "7227": "fb",  # Drosophila melanogaster
1319                "562": "ecocyc",  # Escherichia coli
1320                "11103": None,  # Hepatitis C virus
1321                "9606": "goa_human",  # Homo sapiens
1322                "10090": "mgi",  # Mus musculus
1323                "2104": None,  # Mycoplasma pneumoniae
1324                "4530": "gramene_oryza",  # Oryza sativa
1325                "5833": "GeneDB_Pfalciparum",  # Plasmodium falciparum
1326                "4754": None,  # Pneumocystis carinii
1327                "10116": "rgd",  # Rattus norvegicus
1328                "4932": "sgd",  # Saccharomyces cerevisiae
1329                "4896": "GeneDB_Spombe",  # Schizosaccharomyces pombe
1330                "31033": None,  # Takifugu rubripes
1331                "8355": None,  # Xenopus laevis
1332                "4577": None  # Zea mays
1333                }
1334    version = 1
1335    __shared_state = {"tax": None}
1336
1337    def __init__(self):
1338        self.__dict__ = self.__shared_state
1339        if not self.tax:
1340            path = orngServerFiles.localpath_download("GO", "taxonomy.pickle")
1341            if os.path.isfile(path):
1342                self.tax = cPickle.load(open(path, "rb"))
1343            else:
1344                orngServerFiles.download("GO", "taxonomy.pickle")
1345                self.tax = cPickle.load(open(path, "rb"))
1346
1347    def __getitem__(self, key):
1348        key = self.common_org_map.get(key, key)
1349        return self.code_map[key]
1350
1351    def keys(self):
1352        return list(set(self.common_org_map.keys() + self.code_map.keys()))
1353
1354
1355def from_taxid(id):
1356    """ Return a set of GO organism codes that correspond to NCBI taxonomy id
1357    """
1358    return Taxonomy()[id]
1359
1360
1361def to_taxid(db_code):
1362    """ Return a set of NCBI taxonomy ids from db_code GO organism annotations
1363    """
1364    r = [key for key, val in Taxonomy().code_map.items() if db_code == val]
1365    return set(r)
1366
1367
1368class __progressCallbackWrapper:
1369    def __init__(self, callback):
1370        self.callback = callback
1371
1372    def __call__(self, bCount, bSize, fSize):
1373        fSize = 10000000 if fSize == -1 else fSize
1374        self.callback(100 * bCount * bSize / fSize)
1375
1376from .obiGenomicsUpdate import Update as UpdateBase
1377
1378
1379class Update(UpdateBase):
1380    def __init__(self, local_database_path=None, progressCallback=None):
1381        UpdateBase.__init__(self, local_database_path or getDataDir(), progressCallback)
1382
1383    def CheckModified(self, addr, date=None):
1384        return date < self.GetLastModified(addr) if date else True
1385
1386    def CheckModifiedOrg(self, org):
1387        return self.CheckModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz", self.LastModifiedOrg(org))
1388
1389    def LastModifiedOrg(self, org):
1390        return self.shelve.get((Update.UpdateAnnotation, (org,)), None)
1391
1392    def GetLastModified(self, addr):
1393        stream = urllib2.urlopen(addr)
1394        return datetime.strptime(stream.headers.get("Last-Modified"), "%a, %d %b %Y %H:%M:%S %Z")
1395##        return stream.headers.get("Last-Modified")
1396
1397    def GetAvailableOrganisms(self):
1398        source = urllib2.urlopen("http://www.geneontology.org/gene-associations/").read()
1399        return [s.split(".")[1] for s in sorted(set(re.findall("gene_association\.[a-zA-z0-9_]+?\.gz", source)))]
1400
1401    def GetDownloadedOrganisms(self):
1402        return [name.split(".")[1] for name in os.listdir(self.local_database_path) if name.startswith("gene_association")]
1403
1404    def IsUpdatable(self, func, args):
1405        if func == Update.UpdateOntology:
1406            return self.CheckModified("http://www.geneontology.org/ontology/gene_ontology.obo", self.shelve.get((Update.UpdateOntology, ()), None))
1407        elif func == Update.UpdateAnnotation:
1408            return self.CheckModifiedOrg(args[0])
1409
1410    def GetDownloadable(self):
1411        orgs = set(self.GetAvailableOrganisms()) - set(self.GetDownloadedOrganisms())
1412        ret = []
1413        if (Update.UpdateOntology, ()) not in self.shelve:
1414            ret.append((Update.UpdateOntology, ()))
1415        if orgs:
1416            ret.extend([(Update.UpdateAnnotation, (org,)) for org in orgs])
1417        return ret
1418
1419    def UpdateOntology(self):
1420        Ontology.DownloadOntology(os.path.join(self.local_database_path, "gene_ontology_edit.obo.tar.gz"), self.progressCallback)
1421        self._update(Update.UpdateOntology, (), self.GetLastModified("http://www.geneontology.org/ontology/gene_ontology.obo"))
1422
1423    def UpdateAnnotation(self, org):
1424        Annotations.DownloadAnnotations(org, os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"), self.progressCallback)
1425        self._update(Update.UpdateAnnotation, (org,), self.GetLastModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz"))
1426
1427    def UpdateTaxonomy(self, org):
1428        exclude = ["goa_uniprot", "goa_pdb", "GeneDB_tsetse", "reactome", "goa_zebrafish", "goa_rat", "goa_mouse"]
1429
1430        orgs = self.GetAvailableOrganisms()
1431        tax = defaultdict(set)
1432
1433        for org in orgs:
1434            if org in exclude:
1435                continue
1436            try:
1437                a = obiGO.Annotations(os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"))
1438                taxons = set(ann.taxon for ann in a.annotations)
1439                for taxId in [t.split(":")[-1] for t in taxons if "|" not in t]:  # exclude taxons with cardinality 2
1440                    tax[taxId].add(org)
1441            except Exception, ex:
1442                print ex
1443
1444        cPickle.dump(dict(tax), open(os.path.join(path, "taxonomy.pickle"), "wb"))
1445
1446
1447def _test1():
1448##    Ontology.DownloadOntology("ontology_arch.tar.gz")
1449##    Annotations.DownloadAnnotations("sgd", "annotations_arch.tar.gz")
1450    def _print(f):
1451        print f
1452    o = Ontology("ontology_arch.tar.gz")
1453    a = Annotations("annotations_arch.tar.gz", ontology=o)
1454
1455    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
1456##    profile.runctx("a.GetEnrichedTerms(sorted(a.geneNames)[:100])", {"a":a}, {})
1457    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
1458    d1 = a.GetEnrichedTerms(sorted(a.geneNames)[:1000])  # , progressCallback=_print)
1459
1460##    print a.GetEnrichedTerms(sorted(a.geneNames)[:100])#, progressCallback=_print)
1461
1462
1463def _test2():
1464    o = Ontology()
1465    a = Annotations("human", ontology=o)
1466    clusterGenes = sorted(a.geneNames)[:100]
1467    for i in range(10):
1468        terms = a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["P"])
1469        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["C"])
1470        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["F"])
1471        print i
1472#    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1473
1474##    drawEnrichmentGraph([("bal", 1.0, 5, 6, 0.1, 0.4, ["vv"]),
1475##                        ("GO:0019079", 0.5, 5, 6, 0.1, 0.4, ["cc", "bb"]),
1476##                        ("GO:0022415", 0.4, 5, 7, 0.11, 0.4, ["cc1", "bb"])], open("graph.png", "wb"), None, None)
1477
1478
1479def _test3():
1480    o = Ontology()
1481    a = Annotations("sgd", ontology=o)
1482##    a = Annotations(list(a)[3:len(a)/3], ontology=o)
1483    clusterGenes = sorted(a.geneNames)[:1] + sorted(a.geneNames)[-1:]
1484##    clusterGenes = [g + "_" + str(i%5) for g in sorted(a.geneNames)[:2]]
1485    exonMap = dict([(gene, [gene + "_E%i" % i for i in range(10)]) for gene in a.geneNames])
1486    a.RemapGenes(exonMap)
1487##    o.reverseAliasMapper = o.aliasMapper = {}
1488    terms = a.GetEnrichedTerms(exonMap.values()[0][:2] + exonMap.values()[-1][2:])
1489##    terms = a.GetEnrichedTerms(clusterGenes)
1490    print terms
1491##    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1492    a.DrawEnrichmentGraph(filterByPValue(terms, maxPValue=0.1), len(clusterGenes), len(a.geneNames))
1493
1494if __name__ == "__main__":
1495    _test2()
Note: See TracBrowser for help on using the repository browser.