source: orange-bioinformatics/obiGO.py @ 1334:3c8eb7c872d1

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