source: orange-bioinformatics/Orange/bioinformatics/obiGO.py @ 1632:9cf919d0f343

Revision 1632:9cf919d0f343, 62.1 KB checked in by mitar, 2 years ago (diff)

Fixing imports.

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