source: orange-bioinformatics/obiGO.py @ 1610:f6c73952d74f

Revision 1610:f6c73952d74f, 61.9 KB checked in by Ales Erjavec <ales.erjavec@…>, 2 years ago (diff)

Save download file to a temporary '.part' file.

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