source: orange-bioinformatics/obiGO.py @ 1422:2b381e58061b

Revision 1422:2b381e58061b, 58.8 KB checked in by markotoplak, 3 years ago (diff)

obiGO: GetEnrichedTerms uses all aspects by default.

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