source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1716:4c4266f7c5a5

Revision 1716:4c4266f7c5a5, 62.3 KB checked in by markotoplak, 20 months ago (diff)

Removed the old version of obiKEGG. Renamed obiKEGG2 -> obiKEGG.

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.splitlines():
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.splitlines() 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
627        if not os.path.exists(path):
628            sf = orngServerFiles.ServerFiles()
629            available = sf.listfiles("GO")
630            if file not in available:
631                from . import obiKEGG
632                raise obiKEGG.OrganismNotFoundError(org + str(code))
633            orngServerFiles.download("GO", file)
634
635        return cls(path, ontology=ontology, genematcher=genematcher, progressCallback=progressCallback)
636   
637    def ParseFile(self, file, progressCallback=None):
638        """ Parse and load the annotations from file. Report progress
639        with progressCallback.
640        File can be:
641            - a tarball containing the association file named gene_association
642            - a directory name containing the association file named gene_association
643            - a path to the actual association file
644            - an open file-like object of the association file
645        """
646        if type(file) == str:
647            if os.path.isfile(file) and tarfile.is_tarfile(file):
648                f = tarfile.open(file).extractfile("gene_association")
649            elif os.path.isfile(file) and file.endswith(".gz"):
650                f = gzip.open(file) 
651            elif os.path.isfile(file):
652                f = open(file)
653            else:
654                f = open(os.path.join(file, "gene_association"))
655        else:
656            f = file
657        lines = [line for line in f.read().splitlines() if line.strip()]
658        from Orange.orng import orngMisc
659        milestones = orngMisc.progressBarMilestones(len(lines), 100)
660        for i,line in enumerate(lines):
661            if line.startswith("!"):
662                self.header = self.header + line + "\n"
663                continue
664           
665            a=AnnotationRecord(line)
666            self.AddAnnotation(a)
667#            self.annotations.append(a)
668            if progressCallback and i in milestones:
669                progressCallback(100.0*i/len(lines))
670
671    def AddAnnotation(self, a):
672        """ Add a single `AnotationRecord` instance to this `Annotations`
673        object.
674        """
675        if not isinstance(a, AnnotationRecord):
676            a = AnnotationRecord(a)
677        if not a.geneName or not a.GOId or a.Qualifier == "NOT":
678            return
679       
680        self.geneAnnotations[a.geneName].append(a)
681        self.annotations.append(a)
682        self.termAnnotations[a.GOId].append(a)
683        self.allAnnotations = defaultdict(list)
684       
685        self._geneNames = None
686        self._geneNamesDict = None
687        self._aliasMapper = None
688
689    @property
690    def geneNamesDict(self):
691        if getattr(self, "_geneNamesDict", None) is None:
692            self._geneNamesDict = defaultdict(set)
693            for alias, name in self.aliasMapper.iteritems():
694                self._geneNamesDict[name].add(alias)
695        return self._geneNamesDict
696
697    @property
698    def geneNames(self):
699        if getattr(self, "_geneNames", None) is None:
700            self._geneNames = set([ann.geneName for ann in self.annotations])
701        return self._geneNames
702
703    @property
704    def aliasMapper(self):
705        if getattr(self, "_aliasMapper", None) is None:
706            self._aliasMapper = {}
707            for ann in self.annotations:
708                self._aliasMapper.update([(alias, ann.geneName) for alias in ann.alias + \
709                                          [ann.geneName, ann.DB_Object_ID]])
710        return self._aliasMapper
711
712    def GetGeneNamesTranslator(self, genes):
713        """ Return a dictionary mapping canonical names (DB_Object_Symbol)
714        to `genes`.
715         
716        """
717        def alias(gene):
718            if self.genematcher:
719                return self.genematcher.umatch(gene)
720            else:
721                return gene if gene in self.geneNames else \
722                        self.aliasMapper.get(gene,
723                             self.additionalAliases.get(gene, None))
724                       
725        return dict([(alias(gene), gene) for gene in genes if alias(gene)])
726
727    def _CollectAnnotations(self, id, visited):
728        """ Recursive function collects and caches all annotations for id
729        """
730        if id not in self.allAnnotations and id not in visited:
731            if id in self.ontology.reverseAliasMapper:
732                annotations = [self.termAnnotations.get(alt_id, []) \
733                               for alt_id in self.ontology.reverseAliasMapper[id]] + \
734                                             [self.termAnnotations[id]]
735            else:
736                ## annotations for this term alone
737                annotations = [self.termAnnotations[id]] 
738            visited.add(id)
739            for typeId, child in self.ontology[id].relatedTo:
740                aa = self._CollectAnnotations(child, visited)
741                if type(aa) == set:
742                    ## if it was already reduced in GetAllAnnotations
743                    annotations.append(aa)
744                else:
745                    annotations.extend(aa)
746            self.allAnnotations[id] = annotations
747        return self.allAnnotations[id]
748
749    def GetAllAnnotations(self, id):
750        """ Return a set of all annotations (instances if `AnnotationRectord`)
751        for GO term `id` and all it's subterms.
752       
753        :param id: GO term id
754        :type id: str
755       
756        """
757        visited = set()
758        id = self.ontology.aliasMapper.get(id, id)
759        if id not in self.allAnnotations or type(self.allAnnotations[id]) == list:
760            annot_set = set()
761            for annots in self._CollectAnnotations(id, set()):
762                annot_set.update(annots)
763            self.allAnnotations[id] = annot_set
764        return self.allAnnotations[id]
765
766    def GetAllGenes(self, id, evidenceCodes = None):
767        """ Return a list of genes annotated by specified `evidenceCodes`
768        to GO term 'id' and all it's subterms."
769       
770        :param id: GO term id
771        :type id: str
772       
773        :param evidneceCodes: List of evidence codes to consider when
774                              matching annotations to terms.
775        :type evidenceCodes: list-of-strings
776        """
777        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
778        annotations = self.GetAllAnnotations(id)
779        return list(set([ann.geneName for ann in annotations if ann.Evidence_Code in evidenceCodes]))
780
781    def GetEnrichedTerms(self, genes, reference=None, evidenceCodes=None, 
782                         slimsOnly=False, aspect=None, prob=obiProb.Binomial(),
783                         useFDR=True, progressCallback=None):
784        """ Return a dictionary of enriched terms, with tuples of
785        (list_of_genes, p_value, reference_count) for items and term
786        ids as keys. P-Values are FDR adjusted if useFDR is True (default).
787       
788        :param genes: List of genes
789        :param reference: list of genes (if None all genes included in the
790                          annotations will be used).
791        :param evidenceCodes: List of evidence codes to consider.
792        :param slimsOnly: If `True` return only slim terms
793        :param aspect: Which aspects to use. Use all by default. "P", "F", "C"
794            or a set containing these elements.
795        """
796        revGenesDict = self.GetGeneNamesTranslator(genes)
797        genes = set(revGenesDict.keys())
798        if reference:
799            refGenesDict = self.GetGeneNamesTranslator(reference)
800            reference = set(refGenesDict.keys())
801        else:
802            reference = self.geneNames
803
804        if aspect == None:
805            aspects_set = set(["P", "C", "F"])
806        else:
807            aspects_set = set([aspect]) if isinstance(aspect, basestring) else aspect
808
809        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
810        annotations = [ann for gene in genes for ann in self.geneAnnotations[gene] \
811                       if ann.Evidence_Code in evidenceCodes and ann.Aspect in aspects_set]
812        refAnnotations = set([ann for gene in reference for ann in self.geneAnnotations[gene] \
813                              if ann.Evidence_Code in evidenceCodes and ann.Aspect in aspects_set])
814        annotationsDict = defaultdict(set)
815        for ann in annotations:
816            annotationsDict[ann.GO_ID].add(ann)
817           
818        if slimsOnly and not self.ontology.slimsSubset:
819            import warnings
820            warnings.warn("Unspecified slims subset in the ontology! Using 'goslim_generic' subset", UserWarning)
821            self.ontology.SetSlimsSubset("goslim_generic")
822           
823        terms = self.ontology.ExtractSuperGraph(annotationsDict.keys())
824        res = {}
825        from Orange.orng import orngMisc
826        milestones = orngMisc.progressBarMilestones(len(terms), 100)
827        for i, term in enumerate(terms):
828            if slimsOnly and term not in self.ontology.slimsSubset:
829                continue
830            allAnnotations = self.GetAllAnnotations(term).intersection(refAnnotations)
831##            allAnnotations.intersection_update(refAnnotations)
832            allAnnotatedGenes = set([ann.geneName for ann in allAnnotations])
833            mappedGenes = genes.intersection(allAnnotatedGenes)
834##            if not mappedGenes:
835##                print >> sys.stderr, term, sorted(genes)
836##                print >> sys.stderr, sorted(allAnnotatedGenes)
837##                return
838            if len(reference) > len(allAnnotatedGenes):
839                mappedReferenceGenes = reference.intersection(allAnnotatedGenes)
840            else:
841                mappedReferenceGenes = allAnnotatedGenes.intersection(reference)
842            res[term] = ([revGenesDict[g] for g in mappedGenes],
843                         prob.p_value(len(mappedGenes), len(reference),
844                                      len(mappedReferenceGenes), len(genes)),
845                         len(mappedReferenceGenes))
846            if progressCallback and i in milestones:
847                progressCallback(100.0 * i / len(terms))
848        if useFDR:
849            res = sorted(res.items(), key = lambda (_1, (_2, p, _3)): p)
850            res = dict([(id, (genes, p, ref)) \
851                        for (id, (genes, _, ref)), p in zip(res, obiProb.FDR([p for _, (_, p, _) in res]))])
852        return res
853
854    def GetAnnotatedTerms(self, genes, directAnnotationOnly=False, evidenceCodes=None, progressCallback=None):
855        """ Return all terms that are annotated by genes with evidenceCodes.
856        """
857        genes = [genes] if type(genes) == str else genes
858        revGenesDict = self.GetGeneNamesTranslator(genes)
859        genes = set(revGenesDict.keys())
860        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
861        annotations = [ann for gene in genes for ann in self.geneAnnotations[gene] \
862                       if ann.Evidence_Code in evidenceCodes]
863        dd = defaultdict(set)
864        for ann in annotations:
865            dd[ann.GO_ID].add(revGenesDict.get(ann.geneName, ann.geneName))
866        if not directAnnotationOnly:
867            terms = self.ontology.ExtractSuperGraph(dd.keys())
868            for i, term in enumerate(terms):
869                termAnnots = self.GetAllAnnotations(term).intersection(annotations)
870##                termAnnots.intersection_update(annotations)
871                dd[term].update([revGenesDict.get(ann.geneName, ann.geneName) for ann in termAnnots])
872        return dict(dd)
873
874    def DrawEnrichmentGraph(self, terms, clusterSize, refSize=None, file="graph.png", width=None, height=None, precison=3):
875        refSize = len(self.geneNames) if refSize == None else refSize
876        sortedterms = sorted(terms.items(), key=lambda term:term[1][1])
877        fdr = dict(zip([t[0] for t in sortedterms], obiProb.FDR([t[1][1] for t in sortedterms])))
878        termsList = [(term, (float(len(terms[term][0]))/clusterSize) / (float(terms[term][2])/refSize),
879                      len(terms[term][0]), terms[term][2], terms[term][1],
880                      fdr[term], terms[term][0]) for term in terms]
881                         
882        drawEnrichmentGraph(termsList, file, width, height, ontology=self.ontology, precison=precison)
883
884    def __add__(self, iterable):
885        """ Return a new Annotations object with combined annotations
886        """
887        return Annotations([a for a in self] + [a for a in iterable], ontology=self.ontology)
888
889    def __iadd__(self, iterable):
890        """ Add annotations to this instance
891        """
892        self.extend(iterable)
893        return self
894
895    def __contains__(self, item):
896        return item in self.annotations
897           
898    def __iter__(self):
899        """ Iterate over all AnnotationRecord objects in annotations
900        """
901        return iter(self.annotations)
902
903    def __len__(self):
904        """ Return the number of annotations
905        """
906        return len(self.annotations)
907
908    def __getitem__(self, index):
909        """ Return the i-th annotation record
910        """
911        return self.annotations[index]
912
913    def __getslice__(self, *args):
914        return self.annotations.__getslice__(*args)
915
916    def add(self, line):
917        """ Add one annotation
918        """
919        self.AddAnnotation(line)
920
921    def append(self, line):
922        """ Add one annotation
923        """
924        self.AddAnnotation(line)
925
926    def extend(self, lines):
927        """ Add multiple annotations
928        """
929        for line in lines:
930            self.AddAnnotation(line)
931
932    def RemapGenes(self, map):
933        """
934        """
935        from copy import copy
936        for gene in map:
937            annotations = self.geneAnnotations[gene]
938            for ann in annotations:
939                for name in map[gene]:
940                    ann1 = copy(ann)
941                    ann1.geneName = name
942                    self.add(ann1)
943        self.genematcher = obiGene.GMDirect()
944        try:
945            del self._geneNames
946        except Exception:
947            pass
948        self.genematcher.set_targets(self.geneNames)
949   
950    @staticmethod
951    def DownloadAnnotations(org, file, progressCallback=None):
952        tFile = tarfile.open(file, "w:gz") if type(file) == str else file
953        tmpDir = os.path.join(orngEnviron.bufferDir, "tmp_go/")
954        try:
955            os.mkdir(tmpDir)
956        except Exception:
957            pass
958        fileName = "gene_association." + org + ".gz"
959        urlretrieve("http://www.geneontology.org/gene-associations/" + fileName,
960                    os.path.join(tmpDir, fileName),
961                    progressCallback and __progressCallbackWraper(progressCallback))
962        gzFile = GzipFile(os.path.join(tmpDir, fileName), "r")
963        file = open(os.path.join(tmpDir, "gene_association." + org), "w")
964        file.writelines(gzFile.readlines())
965        file.flush()
966        file.close()
967##        tFile = tarfile.open(os.path.join(tmpDir, "gene_association." + org + ".tar.gz"), "w:gz")
968        tFile.add(os.path.join(tmpDir, "gene_association." + org), "gene_association")
969        annotation = Annotations(os.path.join(tmpDir, "gene_association." + org),
970                    genematcher=obiGene.GMDirect(), progressCallback=progressCallback)
971        cPickle.dump(annotation.geneNames, open(os.path.join(tmpDir, "gene_names.pickle"), "wb"))
972        tFile.add(os.path.join(tmpDir, "gene_names.pickle"), "gene_names.pickle")
973        tFile.close()
974        os.remove(os.path.join(tmpDir, "gene_association." + org))
975        os.remove(os.path.join(tmpDir, "gene_names.pickle"))
976       
977    @staticmethod
978    def DownloadAnnotationsAtRev(org, rev, filename=None, progressCallback=None):
979        import urllib2
980        import shutil
981        if filename is None:
982            filename = os.path.join(default_database_path,
983                                    "gene_association.%s@rev%s.tar.gz" % (code, rev))
984        url = "http://cvsweb.geneontology.org/cgi-bin/cvsweb.cgi/~checkout~/go/gene-associations/gene_association.%s.gz?rev=%s" % (org, rev)
985        url += ";content-type=application%2Fx-gzip"
986        r = urllib2.urlopen(url)
987       
988        with open(filename + ".part", "wb") as f:
989            shutil.copyfileobj(r, f)
990       
991        os.rename(filename + ".part", filename)
992
993from .obiTaxonomy import pickled_cache
994
995@pickled_cache(None, [("GO", "taxonomy.pickle"), ("Taxonomy", "ncbi_taxonomy.tar.gz")])
996def organism_name_search(name):
997    return Annotations.organism_name_search(name)
998
999def filterByPValue(terms, maxPValue=0.1):
1000    """ Filters the terms by the p-value. Asumes terms is a dict with
1001    the same structure as returned from GetEnrichedTerms.
1002   
1003    """
1004    return dict(filter(lambda (k,e): e[1]<=maxPValue, terms.items()))
1005
1006def filterByFrequency(terms, minF=2):
1007    """ Filters the terms by the cluster frequency. Asumes terms is
1008    a dict with the same structure as returned from GetEnrichedTerms.
1009   
1010    """
1011    return dict(filter(lambda (k,e): len(e[0])>=minF, terms.items()))
1012
1013def filterByRefFrequency(terms, minF=4):
1014    """ Filters the terms by the reference frequency. Asumes terms is
1015    a dict with the same structure as returned from GetEnrichedTerms.
1016   
1017    """
1018    return dict(filter(lambda (k,e): e[2]>=minF, terms.items()))
1019
1020##def drawEnrichmentGraph(termsList, clusterSize, refSize, filename="graph.png", width=None, height=None):
1021##    if type(termsList) == dict:
1022##        termsList = [(term, (float(len(termsList[term][0]))/clusterSize) / (float(termsList[term][2])/refSize),
1023##                      len(termsList[term][0]), termsList[term][2], termsList[term][1], 1.0, termsList[term][0]) for term in termsList]
1024##                     
1025##                     
1026##                             
1027##    drawEnrichmentGraph_tostreamMk2(termsList, open(filename, "wb"), width, height)
1028
1029def drawEnrichmentGraph_tostream(GOTerms, clusterSize, refSize, fh, width=None, height=None):
1030    def getParents(term):
1031        parents = extractGODAG([term])
1032        parents = filter(lambda t: t.id in GOTerms and t.id!=term, parents)
1033        c = []
1034        map(c.extend, [getParents(t.id) for t in parents])
1035        parents = filter(lambda t: t not in c, parents)
1036        return parents
1037    parents = dict([(term, getParents(term)) for term in GOTerms])
1038    #print "Parentes", parents
1039    def getChildren(term):
1040        return filter(lambda t: term in [p.id for p in parents[t]], GOTerms.keys())
1041    topLevelTerms = filter(lambda t: not parents[t], parents.keys())
1042    #print "Top level terms", topLevelTerms
1043    termsList=[]
1044    def collect(term, parent):
1045        termsList.append(
1046            ((float(len(GOTerms[term][0]))/clusterSize) / (float(GOTerms[term][2])/refSize),
1047            len(GOTerms[term][0]),
1048            GOTerms[term][2],
1049            "%.4f" % GOTerms[term][1],
1050            loadedGO.termDict[term].name,
1051            loadedGO.termDict[term].id,
1052            ", ".join(GOTerms[term][0]),
1053            parent)
1054            )
1055##        print float(len(GOTerms[term][0])), float(GOTerms[term][2]), clusterSize, refSize
1056        parent = len(termsList)-1
1057        for c in getChildren(term):
1058            collect(c, parent)
1059                         
1060    for topTerm in topLevelTerms:
1061        collect(topTerm, None)
1062
1063    drawEnrichmentGraphPIL_tostream(termsList, fh, width, height)
1064
1065def drawEnrichmentGraph(enriched, file="graph.png", width=None, height=None, header=None, ontology = None, precison=3):
1066    file = open(file, "wb") if type(file) == str else file
1067    drawEnrichmentGraph_tostreamMk2(enriched, file,  width, height, header, ontology, precison)
1068   
1069def drawEnrichmentGraph_tostreamMk2(enriched, fh, width, height, header=None, ontology = None, precison=4):
1070    ontology = ontology if ontology else Ontology()
1071    header = header if header else ["List", "Total", "p-value", "FDR", "Names", "Genes"]
1072    GOTerms = dict([(t[0], t) for t in enriched if t[0] in ontology])
1073    def getParents(term):
1074        parents = ontology.ExtractSuperGraph([term])
1075        parents = [id for id in parents if id in GOTerms and id != term]
1076        c = reduce(set.union, [set(ontology.ExtractSuperGraph([id])) - set([id]) for id in parents], set())
1077        parents = [t for t in parents if t not in c]
1078        return parents
1079    parents = dict([(term, getParents(term)) for term in GOTerms])
1080    #print "Parentes", parents
1081    def getChildren(term):
1082        return [id for id in GOTerms if term in parents[id]]
1083    topLevelTerms = [id for id in parents if not parents[id]]
1084    #print "Top level terms", topLevelTerms
1085    termsList=[]
1086    fmt = "%" + ".%if" % precison
1087    def collect(term, parent):
1088##        termsList.append(
1089##            ((float(len(GOTerms[term][0]))/clusterSize) / (float(GOTerms[term][2])/refSize),
1090##            len(GOTerms[term][0]),
1091##            GOTerms[term][2],
1092##            "%.4f" % GOTerms[term][1],
1093##            loadedGO.termDict[term].name,
1094##            loadedGO.termDict[term].id,
1095##            ", ".join(GOTerms[term][0]),
1096##            parent)
1097##            )
1098        termsList.append(GOTerms[term][1:4] + \
1099                         (fmt % GOTerms[term][4],
1100                          fmt % GOTerms[term][5],
1101                          ontology[term].name,
1102                          ", ".join(GOTerms[term][6])) + (parent,))
1103##        print float(len(GOTerms[term][0])), float(GOTerms[term][2]), clusterSize, refSize
1104        parent = len(termsList)-1
1105        for c in getChildren(term):
1106            collect(c, parent)
1107                         
1108    for topTerm in topLevelTerms:
1109        collect(topTerm, None)
1110    for entry in enriched:
1111        if entry[0] not in ontology:
1112            termsList.append(entry[1:4] + \
1113                             (fmt % entry[4],
1114                              fmt % entry[5],
1115                              entry[0],
1116                              ", ".join(entry[6])) + (None,))
1117
1118    drawEnrichmentGraphPIL_tostream(termsList, header, fh, width, height)
1119##    drawEnrichmentGraphPylab_tostream(termsList, header, fh, width, height)
1120   
1121def drawEnrichmentGraphPIL_tostream(termsList, headers, fh, width=None, height=None):
1122    from PIL import Image, ImageDraw, ImageFont
1123    backgroundColor = (255, 255, 255)
1124    textColor = (0, 0, 0)
1125    graphColor = (0, 0, 255)
1126    fontSize = height==None and 12 or (height-60)/len(termsList)
1127    font = ImageFont.load_default()
1128    try:
1129        font = ImageFont.truetype("arial.ttf", fontSize)
1130    except:
1131        pass
1132    getMaxTextHeightHint = lambda l: max([font.getsize(t)[1] for t in l])
1133    getMaxTextWidthHint = lambda l: max([font.getsize(t)[0] for t in l])
1134    maxFoldWidth = width!=None and min(150, width/6) or 150
1135    maxFoldEnrichment = max([t[0] for t in termsList])
1136    foldNormalizationFactor = float(maxFoldWidth)/maxFoldEnrichment
1137    foldWidths = [int(foldNormalizationFactor*term[0]) for term in termsList]
1138    treeStep = 10
1139    treeWidth = {}
1140    for i, term in enumerate(termsList):
1141        treeWidth[i] = (term[-1]==None and 1 or treeWidth[term[-1]]+1)
1142    treeStep = width!=None and min(treeStep, width/(6*max(treeWidth.values())) or 2) or treeStep
1143    treeWidth = [w*treeStep + foldWidths[i] for i, w in treeWidth.items()]
1144    treeWidth = max(treeWidth) - maxFoldWidth
1145    verticalMargin = 10
1146    horizontalMargin = 10
1147##    print verticalMargin, maxFoldWidth, treeWidth
1148##    treeWidth = 100
1149    firstColumnStart = verticalMargin + maxFoldWidth + treeWidth + 10
1150    secondColumnStart = firstColumnStart + getMaxTextWidthHint([str(t[1]) for t in termsList]+[headers[0]]) + 2
1151    thirdColumnStart = secondColumnStart + getMaxTextWidthHint([str(t[2]) for t in termsList]+[headers[1]]) + 2
1152    fourthColumnStart = thirdColumnStart + getMaxTextWidthHint([str(t[3]) for t in termsList]+[headers[2]]) + 2
1153    fifthColumnStart = fourthColumnStart + getMaxTextWidthHint([str(t[4]) for t in termsList]+[headers[3]]) + 4
1154##    maxAnnotationTextWidth = width==None and getMaxTextWidthHint([str(t[4]) for t in termsList]+["Annotation"]) or (width - fourthColumnStart - verticalMargin) * 2 / 3
1155    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]]))
1156    sixthColumnStart  = fifthColumnStart + maxAnnotationTextWidth + 4
1157    maxGenesTextWidth = width==None and getMaxTextWidthHint([str(t[6]) for t in termsList]+[headers[5]]) or (width - fifthColumnStart - verticalMargin) / 3
1158   
1159    legendHeight = font.getsize("1234567890")[1]*2
1160    termHeight = font.getsize("A")[1]
1161##    print fourthColumnStart, maxAnnotationTextWidth, verticalMargin
1162    width = sixthColumnStart + maxGenesTextWidth + verticalMargin
1163    height = len(termsList)*termHeight+2*(legendHeight+horizontalMargin)
1164
1165    image = Image.new("RGB", (width, height), backgroundColor)
1166    draw = ImageDraw.Draw(image)
1167
1168    def truncText(text, maxWidth, append=""):
1169        #print getMaxTextWidthHint([text]), maxAnnotationTextWidth
1170        if getMaxTextWidthHint([text])>maxWidth:
1171            while getMaxTextWidthHint([text+"..."+append])>maxWidth and text:
1172                text = text[:-1]
1173            if text:
1174                text = text+"..."+append
1175            else:
1176                text = append
1177        return text
1178    currentY = horizontalMargin + legendHeight
1179    connectAtX = {}
1180    for i, term in enumerate(termsList):
1181        draw.line([(verticalMargin, currentY+termHeight/2), (verticalMargin + foldWidths[i], currentY+termHeight/2)], width=termHeight-2, fill=graphColor)
1182        draw.text((firstColumnStart, currentY), str(term[1]), font=font, fill=textColor)
1183        draw.text((secondColumnStart, currentY), str(term[2]), font=font, fill=textColor)
1184        draw.text((thirdColumnStart, currentY), str(term[3]), font=font, fill=textColor)
1185        draw.text((fourthColumnStart, currentY), str(term[4]), font=font, fill=textColor)
1186##        annotText = width!=None and truncText(str(term[5]), maxAnnotationTextWidth, str(term[5])) or str(term[4])
1187        annotText = width!=None and truncText(str(term[5]), maxAnnotationTextWidth)
1188        draw.text((fifthColumnStart, currentY), annotText, font=font, fill=textColor)
1189        genesText = width!=None and truncText(str(term[6]), maxGenesTextWidth) or str(term[6])
1190        draw.text((sixthColumnStart, currentY), genesText, font=font, fill=textColor)
1191        lineEnd = term[-1]==None and firstColumnStart-10 or connectAtX[term[-1]]
1192        draw.line([(verticalMargin+foldWidths[i]+1, currentY+termHeight/2), (lineEnd, currentY+termHeight/2)], width=1, fill=textColor)
1193        if term[-1]!=None:
1194            draw.line([(lineEnd, currentY+termHeight/2), (lineEnd, currentY+termHeight/2 - termHeight*(i-term[-1]))], width=1, fill=textColor)
1195        connectAtX[i] = lineEnd - treeStep
1196        currentY+=termHeight
1197
1198    currentY = horizontalMargin
1199    draw.text((firstColumnStart, currentY), headers[0], font=font, fill=textColor)
1200    draw.text((secondColumnStart, currentY), headers[1], font=font, fill=textColor)
1201    draw.text((thirdColumnStart, currentY), headers[2], font=font, fill=textColor)
1202    draw.text((fourthColumnStart, currentY), headers[3], font=font, fill=textColor)
1203    draw.text((fifthColumnStart, currentY), headers[4], font=font, fill=textColor)
1204    draw.text((sixthColumnStart, currentY), headers[5], font=font, fill=textColor)
1205
1206    horizontalMargin = 0
1207    #draw.line([(verticalMargin, height - horizontalMargin - legendHeight), (verticalMargin + maxFoldWidth, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1208    draw.line([(verticalMargin, horizontalMargin + legendHeight), (verticalMargin + maxFoldWidth, horizontalMargin + legendHeight)], width=1, fill=textColor)
1209    maxLabelWidth = getMaxTextWidthHint([" "+str(i) for i in range(int(maxFoldEnrichment+1))])
1210    numOfLegendLabels = max(int(maxFoldWidth/maxLabelWidth), 2)
1211    for i in range(numOfLegendLabels+1):
1212        #draw.line([(verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight/2), (verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1213        #draw.text((verticalMargin + i*maxFoldWidth/10 - font.getsize(str(i))[0]/2, height - horizontalMargin - legendHeight/2), str(i), font=font, fill=textColor)
1214
1215        label = str(int(i*maxFoldEnrichment/numOfLegendLabels))
1216        draw.line([(verticalMargin + i*maxFoldWidth/numOfLegendLabels, horizontalMargin + legendHeight/2), (verticalMargin + i*maxFoldWidth/numOfLegendLabels, horizontalMargin + legendHeight)], width=1, fill=textColor)
1217        draw.text((verticalMargin + i*maxFoldWidth/numOfLegendLabels - font.getsize(label)[0]/2, horizontalMargin), label, font=font, fill=textColor)
1218       
1219    image.save(fh)
1220
1221def drawEnrichmentGraphPylab_tostream(termsList, headers, fh, width=None, height=None, show=True):
1222    from matplotlib import pyplot as plt
1223    from matplotlib.patches import Rectangle
1224   
1225    maxFoldWidth = width!=None and min(150, width/6) or 150
1226    maxFoldEnrichment = max([t[0] for t in termsList])
1227    foldNormalizationFactor = float(maxFoldWidth)/maxFoldEnrichment
1228##    foldWidths = [int(foldNormalizationFactor*term[0]) for term in termsList]
1229    foldWidths = [term[0] for term in termsList]
1230    treeStep = maxFoldEnrichment*0.05
1231    treeWidth = {}
1232
1233    for i, term in enumerate(termsList):
1234        treeWidth[i] = (term[-1]==None and treeStep or treeWidth[term[-1]] + treeStep)
1235    maxTreeWidth = max(treeWidth)
1236
1237    connectAt = {}
1238    cellText = []
1239    axes1 = plt.axes([0.1, 0.1, 0.2, 0.8])
1240    for i, line in enumerate(termsList):
1241        enrichment, n, m, p_val, fdr_val, name, genes, parent = line
1242        r = Rectangle((0, len(termsList) - i - 0.4), enrichment, 0.8)
1243        plt.gca().add_patch(r)
1244        plt.plot([enrichment, connectAt.get(parent, maxFoldEnrichment + maxTreeWidth)], [len(termsList) - i, len(termsList) - i], color="black")
1245        connectAt[i] = connectAt.get(parent, maxFoldEnrichment + maxTreeWidth) - treeStep
1246        if parent != None:
1247            plt.plot([connectAt.get(parent)]*2, [len(termsList) - i, len(termsList) - parent], color="black")
1248        cellText.append((str(n), str(m), p_val, fdr_val, name, genes))
1249
1250##    from Orange.orng.orngClustering import TableTextLayout
1251##    text = TableTextLayout((maxFoldEnrichment*1.1, len(termsList)), cellText)
1252    from Orange.orng.orngClustering import TablePlot
1253    if True:
1254        axes2 = plt.axes([0.3, 0.1, 0.6, 0.8], sharey=axes1)
1255        axes2.set_axis_off()
1256        table = TablePlot((0, len(termsList)), axes=plt.gca())
1257        for i, line in enumerate(cellText):
1258            for j, text in enumerate(line):
1259                table.add_cell(i, j,width=len(text), height=1, text=text, loc="left", edgecolor="w", facecolor="w")
1260
1261        table.set_figure(plt.gcf())
1262        plt.gca().add_artist(table)
1263        plt.gca()._set_artist_props(table)
1264##    plt.text(3, 3, "\n".join(["\t".join(text) for text in cellText]))
1265
1266##    table = plt.table(cellText=cellText, colLabels=headers, loc="right")
1267##    table.set_transform(plt.gca().transData)
1268##   
1269##    table.set_xy(20,20)
1270    plt.show()
1271   
1272class Taxonomy(object):
1273    """Maps NCBI taxonomy ids to coresponding GO organism codes
1274    """
1275    common_org_map = {"297284":"9913", "30523":"9913", # Bos taurus
1276                      "5782":"352472", "44689":"352472", "366501":"352472", # Dictyostelium discoideum
1277                      "83333": "562", # Escherichia coli
1278                      "52545":"4530", "4532":"4530", "65489":"4530", "4533":"4530", "77588":"4530", "29689":"4530",
1279                      "4538":"4530", "40148":"4530", "29690":"4530", "110450":"4530", "4534":"4530", "83309":"4530",
1280                      "4528":"4530", "127571":"4530", "40149":"4530", "83307":"4530", "63629":"4530", "4536": "4530",
1281                      "4535":"4530", "4537":"4530", "65491":"4530", "83308":"4530", "4529":"4530", "4530":"4530",
1282                      "39946":"4530", "39947":"4530", "110451":"4530", "364100":"4530", "364099":"4530", "4539":"4530",
1283                      }
1284    code_map = {"3702":"tair",  # Arabidopsis thaliana
1285                "9913":"goa_cow", # Bos taurus
1286                "6239":"wb", # Caenorhabditis elegans
1287                "3055":None, # Chlamydomonas reinhardtii
1288                "7955":"zfin", # Danio rerio (zebrafish)
1289                "352472":"dictyBase", # Dictyostelium discoideum
1290                "7227": "fb", # Drosophila melanogaster
1291                "562": "ecocyc", # Escherichia coli
1292                "11103": None, # Hepatitis C virus
1293                "9606": "goa_human",  # Homo sapiens
1294                "10090": "mgi", # Mus musculus
1295                "2104": None,  # Mycoplasma pneumoniae
1296                "4530": "gramene_oryza",  # Oryza sativa
1297                "5833": "GeneDB_Pfalciparum",  # Plasmodium falciparum
1298                "4754": None,  # Pneumocystis carinii
1299                "10116": "rgd", # Rattus norvegicus
1300                "4932": "sgd",  # Saccharomyces cerevisiae
1301                "4896": "GeneDB_Spombe", # Schizosaccharomyces pombe
1302                "31033": None, # Takifugu rubripes
1303                "8355": None,  # Xenopus laevis
1304                "4577": None # Zea mays
1305                }
1306    version = 1
1307    __shared_state = {"tax": None}
1308    def __init__(self):
1309        self.__dict__ = self.__shared_state
1310        if not self.tax:
1311            from Orange.orng import orngServerFiles
1312            path = orngServerFiles.localpath_download("GO", "taxonomy.pickle")
1313            if os.path.isfile(path):
1314                self.tax = cPickle.load(open(path, "rb"))
1315            else:
1316                orngServerFiles.download("GO", "taxonomy.pickle")
1317                self.tax = cPickle.load(open(path, "rb"))
1318               
1319    def __getitem__(self, key):
1320        key = self.common_org_map.get(key, key)
1321        return self.code_map[key]
1322   
1323    def keys(self):
1324        return list(set(self.common_org_map.keys() + self.code_map.keys()))
1325   
1326#    @classmethod
1327#    def get_taxonomy(cls):
1328#        import urllib2 as url
1329#        import sgmllib
1330#        organisms
1331#        class MyParser(sgmllib.SGMLParser):
1332#            inTable = False
1333#            def start_table(self, attributes):
1334#                self.inTable = dict(attributes).get("summary", False)
1335#            def end_table(self):
1336#                self.inTable = False
1337#            def start
1338   
1339def from_taxid(id):
1340    """ Return a set of GO organism codes that correspond to NCBI taxonomy id
1341    """
1342    return Taxonomy()[id]
1343
1344def to_taxid(db_code):
1345    """ Return a set of NCBI taxonomy ids from db_code GO organism annotations
1346    """
1347    r = [key for key, val in Taxonomy().code_map.items() if db_code == val]
1348    return set(r)
1349   
1350
1351class __progressCallbackWrapper:
1352    def __init__(self, callback):
1353        self.callback = callback
1354    def __call__(self, bCount, bSize, fSize):
1355        fSize = 10000000 if fSize == -1 else fSize
1356        self.callback(100*bCount*bSize/fSize)
1357       
1358from .obiGenomicsUpdate import Update as UpdateBase
1359
1360import urllib2
1361
1362class Update(UpdateBase):
1363    def __init__(self, local_database_path=None, progressCallback=None):
1364        UpdateBase.__init__(self, local_database_path or getDataDir(), progressCallback)
1365    def CheckModified(self, addr, date=None):
1366        return date < self.GetLastModified(addr) if date else True
1367       
1368    def CheckModifiedOrg(self, org):
1369        return self.CheckModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz", self.LastModifiedOrg(org))
1370   
1371    def LastModifiedOrg(self, org):
1372        return self.shelve.get((Update.UpdateAnnotation, (org,)), None)
1373
1374    def GetLastModified(self, addr):
1375        stream = urllib2.urlopen(addr)
1376        return datetime.strptime(stream.headers.get("Last-Modified"), "%a, %d %b %Y %H:%M:%S %Z")
1377##        return stream.headers.get("Last-Modified")
1378
1379    def GetAvailableOrganisms(self):
1380        source = urllib2.urlopen("http://www.geneontology.org/gene-associations/").read()
1381        return [s.split(".")[1] for s in sorted(set(re.findall("gene_association\.[a-zA-z0-9_]+?\.gz", source)))]
1382
1383    def GetDownloadedOrganisms(self):
1384        return [name.split(".")[1] for name in os.listdir(self.local_database_path) if name.startswith("gene_association")]
1385
1386    def IsUpdatable(self, func, args):
1387        if func == Update.UpdateOntology:
1388            return self.CheckModified("http://www.geneontology.org/ontology/gene_ontology.obo", self.shelve.get((Update.UpdateOntology, ()), None))
1389        elif func == Update.UpdateAnnotation:
1390            return self.CheckModifiedOrg(args[0])
1391           
1392    def GetDownloadable(self):
1393        orgs = set(self.GetAvailableOrganisms()) - set(self.GetDownloadedOrganisms())
1394        ret = []
1395        if (Update.UpdateOntology, ()) not in self.shelve:
1396            ret.append((Update.UpdateOntology, ()))
1397        if orgs:
1398            ret.extend([(Update.UpdateAnnotation, (org,)) for org in orgs])
1399        return ret
1400
1401    def UpdateOntology(self):
1402        Ontology.DownloadOntology(os.path.join(self.local_database_path, "gene_ontology_edit.obo.tar.gz"), self.progressCallback)
1403        self._update(Update.UpdateOntology, (), self.GetLastModified("http://www.geneontology.org/ontology/gene_ontology.obo"))
1404
1405    def UpdateAnnotation(self, org):
1406        Annotations.DownloadAnnotations(org, os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"), self.progressCallback)
1407        self._update(Update.UpdateAnnotation, (org,), self.GetLastModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz"))
1408       
1409    def UpdateTaxonomy(self, org):
1410        exclude = ["goa_uniprot", "goa_pdb", "GeneDB_tsetse", "reactome", "goa_zebrafish", "goa_rat", "goa_mouse"]
1411
1412        orgs = self.GetAvailableOrganisms()
1413        tax = defaultdict(set)
1414
1415        for org in orgs:
1416            if org in exclude:
1417                continue
1418            try:
1419                a = obiGO.Annotations(os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"))
1420                taxons = set(ann.taxon for ann in a.annotations)
1421                for taxId in [t.split(":")[-1] for t in taxons if "|" not in t]: ## exclude taxons with cardinality 2
1422                    tax[taxId].add(org)
1423            except Exception, ex:
1424                print ex
1425               
1426        cPickle.dump(dict(tax), open(os.path.join(path, "taxonomy.pickle"), "wb"))
1427           
1428
1429def _test1():
1430##    Ontology.DownloadOntology("ontology_arch.tar.gz")
1431##    Annotations.DownloadAnnotations("sgd", "annotations_arch.tar.gz")
1432    def _print(f):
1433        print f
1434    o = Ontology("ontology_arch.tar.gz")
1435    a = Annotations("annotations_arch.tar.gz", ontology=o)
1436   
1437    a.GetEnrichedTerms(sorted(a.geneNames)[:100])#, progressCallback=_print)
1438##    profile.runctx("a.GetEnrichedTerms(sorted(a.geneNames)[:100])", {"a":a}, {})
1439    a.GetEnrichedTerms(sorted(a.geneNames)[:100])#, progressCallback=_print)
1440    d1 = a.GetEnrichedTerms(sorted(a.geneNames)[:1000])#, progressCallback=_print)
1441   
1442##    print a.GetEnrichedTerms(sorted(a.geneNames)[:100])#, progressCallback=_print)
1443
1444def _test2():
1445    o = Ontology()
1446    a = Annotations("human", ontology=o)
1447    clusterGenes = sorted(a.geneNames)[:100]
1448    for i in range(10):
1449        terms = a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["P"])
1450        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["C"])
1451        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["F"])
1452        print i
1453#    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1454             
1455##    drawEnrichmentGraph([("bal", 1.0, 5, 6, 0.1, 0.4, ["vv"]),
1456##                        ("GO:0019079", 0.5, 5, 6, 0.1, 0.4, ["cc", "bb"]),
1457##                        ("GO:0022415", 0.4, 5, 7, 0.11, 0.4, ["cc1", "bb"])], open("graph.png", "wb"), None, None)
1458
1459def _test3():
1460    o = Ontology()
1461    a = Annotations("sgd", ontology=o)
1462##    a = Annotations(list(a)[3:len(a)/3], ontology=o)
1463    clusterGenes = sorted(a.geneNames)[:1] + sorted(a.geneNames)[-1:]
1464##    clusterGenes = [g + "_" + str(i%5) for g in sorted(a.geneNames)[:2]]
1465    exonMap = dict([(gene, [gene+"_E%i" %i for i in range(10)]) for gene in a.geneNames])
1466    a.RemapGenes(exonMap)
1467##    o.reverseAliasMapper = o.aliasMapper = {}
1468    terms = a.GetEnrichedTerms(exonMap.values()[0][:2] + exonMap.values()[-1][2:])
1469##    terms = a.GetEnrichedTerms(clusterGenes)
1470    print terms
1471##    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1472    a.DrawEnrichmentGraph(filterByPValue(terms, maxPValue=0.1), len(clusterGenes), len(a.geneNames))
1473   
1474if __name__ == "__main__":
1475    _test2()
Note: See TracBrowser for help on using the repository browser.