source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1767:214b1db4e31a

Revision 1767:214b1db4e31a, 61.2 KB checked in by Ales Erjavec <ales.erjavec@…>, 12 months ago (diff)

obiGO code style fixes.

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