source: orange-bioinformatics/obiGO.py @ 1256:84a64c254d74

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