source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1769:c0c6b64879cf

Revision 1769:c0c6b64879cf, 60.1 KB checked in by Ales Erjavec <ales.erjavec@…>, 12 months ago (diff)

Filter out terms from the annotations that are not present in the ontology.

If the annotations are newer then the ontology they can contain terms not yet
present in the ontology. This leads to an error when extracting the term's
parents.

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 = annotationsDict.keys()
875        filteredTerms = [term for term in terms if term in self.ontology]
876
877        if len(terms) != len(filteredTerms):
878            termDiff = set(terms) - set(filteredTerms)
879            warnings.warn("%s terms in the annotations were not found in the "
880                          "ontology." % ",".join(map(repr, termDiff)),
881                          UserWarning)
882
883        terms = self.ontology.ExtractSuperGraph(filteredTerms)
884        res = {}
885
886        milestones = orngMisc.progressBarMilestones(len(terms), 100)
887        for i, term in enumerate(terms):
888            if slimsOnly and term not in self.ontology.slimsSubset:
889                continue
890            allAnnotations = self.GetAllAnnotations(term).intersection(refAnnotations)
891##            allAnnotations.intersection_update(refAnnotations)
892            allAnnotatedGenes = set([ann.geneName for ann in allAnnotations])
893            mappedGenes = genes.intersection(allAnnotatedGenes)
894##            if not mappedGenes:
895##                print >> sys.stderr, term, sorted(genes)
896##                print >> sys.stderr, sorted(allAnnotatedGenes)
897##                return
898            if len(reference) > len(allAnnotatedGenes):
899                mappedReferenceGenes = reference.intersection(allAnnotatedGenes)
900            else:
901                mappedReferenceGenes = allAnnotatedGenes.intersection(reference)
902            res[term] = ([revGenesDict[g] for g in mappedGenes],
903                         prob.p_value(len(mappedGenes), len(reference),
904                                      len(mappedReferenceGenes), len(genes)),
905                         len(mappedReferenceGenes))
906            if progressCallback and i in milestones:
907                progressCallback(100.0 * i / len(terms))
908        if useFDR:
909            res = sorted(res.items(), key=lambda (_1, (_2, p, _3)): p)
910            res = dict([(id, (genes, p, ref))
911                        for (id, (genes, _, ref)), p in
912                        zip(res, obiProb.FDR([p for _, (_, p, _) in res]))])
913        return res
914
915    def GetAnnotatedTerms(self, genes, directAnnotationOnly=False,
916                          evidenceCodes=None, progressCallback=None):
917        """ Return all terms that are annotated by genes with evidenceCodes.
918        """
919        genes = [genes] if type(genes) == str else genes
920        revGenesDict = self.GetGeneNamesTranslator(genes)
921        genes = set(revGenesDict.keys())
922        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
923        annotations = [ann for gene in genes for ann in self.geneAnnotations[gene]
924                       if ann.Evidence_Code in evidenceCodes]
925        dd = defaultdict(set)
926        for ann in annotations:
927            dd[ann.GO_ID].add(revGenesDict.get(ann.geneName, ann.geneName))
928        if not directAnnotationOnly:
929            terms = dd.keys()
930            filteredTerms = [term for term in terms if term in self.ontology]
931            if len(terms) != len(filteredTerms):
932                termDiff = set(terms) - set(filteredTerms)
933                warnings.warn(
934                    "%s terms in the annotations were not found in the "
935                    "ontology." % ",".join(map(repr, termDiff)),
936                    UserWarning)
937
938            terms = self.ontology.ExtractSuperGraph(filteredTerms)
939            for i, term in enumerate(terms):
940                termAnnots = self.GetAllAnnotations(term).intersection(annotations)
941##                termAnnots.intersection_update(annotations)
942                dd[term].update([revGenesDict.get(ann.geneName, ann.geneName)
943                                 for ann in termAnnots])
944        return dict(dd)
945
946    def DrawEnrichmentGraph(self, terms, clusterSize, refSize=None,
947                            file="graph.png", width=None, height=None,
948                            precison=3):
949        refSize = len(self.geneNames) if refSize == None else refSize
950        sortedterms = sorted(terms.items(), key=lambda term: term[1][1])
951        fdr = dict(zip([t[0] for t in sortedterms],
952                       obiProb.FDR([t[1][1] for t in sortedterms])))
953        termsList = [(term,
954                      ((float(len(terms[term][0])) / clusterSize) /
955                       (float(terms[term][2]) / refSize)),
956                      len(terms[term][0]),
957                      terms[term][2],
958                      terms[term][1],
959                      fdr[term],
960                      terms[term][0])
961                     for term in terms]
962
963        drawEnrichmentGraph(termsList, file, width, height,
964                            ontology=self.ontology, precison=precison)
965
966    def __add__(self, iterable):
967        """ Return a new Annotations object with combined annotations
968        """
969        return Annotations([a for a in self] + [a for a in iterable],
970                           ontology=self.ontology)
971
972    def __iadd__(self, iterable):
973        """ Add annotations to this instance
974        """
975        self.extend(iterable)
976        return self
977
978    def __contains__(self, item):
979        return item in self.annotations
980
981    def __iter__(self):
982        """ Iterate over all AnnotationRecord objects in annotations
983        """
984        return iter(self.annotations)
985
986    def __len__(self):
987        """ Return the number of annotations
988        """
989        return len(self.annotations)
990
991    def __getitem__(self, index):
992        """ Return the i-th annotation record
993        """
994        return self.annotations[index]
995
996    def __getslice__(self, *args):
997        return self.annotations.__getslice__(*args)
998
999    def add(self, line):
1000        """ Add one annotation
1001        """
1002        self.AddAnnotation(line)
1003
1004    def append(self, line):
1005        """ Add one annotation
1006        """
1007        self.AddAnnotation(line)
1008
1009    def extend(self, lines):
1010        """ Add multiple annotations
1011        """
1012        for line in lines:
1013            self.AddAnnotation(line)
1014
1015    def RemapGenes(self, map):
1016        """
1017        """
1018        for gene in map:
1019            annotations = self.geneAnnotations[gene]
1020            for ann in annotations:
1021                for name in map[gene]:
1022                    ann1 = copy.copy(ann)
1023                    ann1.geneName = name
1024                    self.add(ann1)
1025        self.genematcher = obiGene.GMDirect()
1026        try:
1027            del self._geneNames
1028        except Exception:
1029            pass
1030        self.genematcher.set_targets(self.geneNames)
1031
1032    @staticmethod
1033    def DownloadAnnotations(org, file, progressCallback=None):
1034        tFile = tarfile.open(file, "w:gz") if type(file) == str else file
1035        tmpDir = os.path.join(orngEnviron.bufferDir, "tmp_go/")
1036        try:
1037            os.mkdir(tmpDir)
1038        except Exception:
1039            pass
1040        fileName = "gene_association." + org + ".gz"
1041        urlretrieve(("http://www.geneontology.org/gene-associations/" +
1042                     fileName),
1043                    os.path.join(tmpDir, fileName),
1044                    progressCallback and __progressCallbackWrapper(progressCallback))
1045        gzFile = GzipFile(os.path.join(tmpDir, fileName), "r")
1046        file = open(os.path.join(tmpDir, "gene_association." + org), "w")
1047        file.writelines(gzFile.readlines())
1048        file.flush()
1049        file.close()
1050
1051        tFile.add(os.path.join(tmpDir, "gene_association." + org),
1052                  "gene_association")
1053        annotation = Annotations(os.path.join(tmpDir, "gene_association." + org),
1054                    genematcher=obiGene.GMDirect(), progressCallback=progressCallback)
1055        cPickle.dump(annotation.geneNames, open(os.path.join(tmpDir, "gene_names.pickle"), "wb"))
1056        tFile.add(os.path.join(tmpDir, "gene_names.pickle"), "gene_names.pickle")
1057        tFile.close()
1058        os.remove(os.path.join(tmpDir, "gene_association." + org))
1059        os.remove(os.path.join(tmpDir, "gene_names.pickle"))
1060
1061    @staticmethod
1062    def DownloadAnnotationsAtRev(org, rev, filename=None,
1063                                 progressCallback=None):
1064        if filename is None:
1065            filename = os.path.join(default_database_path,
1066                                    "gene_association.%s@rev%s.tar.gz" %
1067                                    (code, rev))
1068        url = ("http://cvsweb.geneontology.org/cgi-bin/cvsweb.cgi/~checkout~/go/gene-associations/gene_association.%s.gz?rev=%s" %
1069               (org, rev))
1070        url += ";content-type=application%2Fx-gzip"
1071        r = urllib2.urlopen(url)
1072
1073        with open(filename + ".part", "wb") as f:
1074            shutil.copyfileobj(r, f)
1075
1076        os.rename(filename + ".part", filename)
1077
1078from .obiTaxonomy import pickled_cache
1079
1080
1081@pickled_cache(None, [("GO", "taxonomy.pickle"),
1082                      ("Taxonomy", "ncbi_taxonomy.tar.gz")])
1083def organism_name_search(name):
1084    return Annotations.organism_name_search(name)
1085
1086
1087def filterByPValue(terms, maxPValue=0.1):
1088    """ Filters the terms by the p-value. Asumes terms is a dict with
1089    the same structure as returned from GetEnrichedTerms.
1090
1091    """
1092    return dict(filter(lambda (k, e): e[1] <= maxPValue, terms.items()))
1093
1094
1095def filterByFrequency(terms, minF=2):
1096    """ Filters the terms by the cluster frequency. Asumes terms is
1097    a dict with the same structure as returned from GetEnrichedTerms.
1098
1099    """
1100    return dict(filter(lambda (k, e): len(e[0]) >= minF, terms.items()))
1101
1102
1103def filterByRefFrequency(terms, minF=4):
1104    """ Filters the terms by the reference frequency. Asumes terms is
1105    a dict with the same structure as returned from GetEnrichedTerms.
1106
1107    """
1108    return dict(filter(lambda (k, e): e[2] >= minF, terms.items()))
1109
1110
1111def drawEnrichmentGraph(enriched, file="graph.png", width=None, height=None,
1112                        header=None, ontology=None, precison=3):
1113    file = open(file, "wb") if type(file) == str else file
1114    drawEnrichmentGraph_tostreamMk2(enriched, file, width, height, header,
1115                                    ontology, precison)
1116
1117
1118def drawEnrichmentGraph_tostreamMk2(enriched, fh, width, height, header=None,
1119                                    ontology=None, precison=4):
1120    ontology = ontology if ontology else Ontology()
1121    header = header if header else ["List", "Total", "p-value", "FDR",
1122                                    "Names", "Genes"]
1123    GOTerms = dict([(t[0], t) for t in enriched if t[0] in ontology])
1124
1125    def getParents(term):
1126        parents = ontology.ExtractSuperGraph([term])
1127        parents = [id for id in parents if id in GOTerms and id != term]
1128        c = reduce(set.union, [set(ontology.ExtractSuperGraph([id])) -
1129                               set([id]) for id in parents], set())
1130        parents = [t for t in parents if t not in c]
1131        return parents
1132    parents = dict([(term, getParents(term)) for term in GOTerms])
1133    # print "Parentes", parents
1134
1135    def getChildren(term):
1136        return [id for id in GOTerms if term in parents[id]]
1137    topLevelTerms = [id for id in parents if not parents[id]]
1138    # print "Top level terms", topLevelTerms
1139    termsList = []
1140    fmt = "%" + ".%if" % precison
1141
1142    def collect(term, parent):
1143        termsList.append(GOTerms[term][1:4] + \
1144                         (fmt % GOTerms[term][4],
1145                          fmt % GOTerms[term][5],
1146                          ontology[term].name,
1147                          ", ".join(GOTerms[term][6])) + (parent,))
1148        parent = len(termsList) - 1
1149        for c in getChildren(term):
1150            collect(c, parent)
1151
1152    for topTerm in topLevelTerms:
1153        collect(topTerm, None)
1154    for entry in enriched:
1155        if entry[0] not in ontology:
1156            termsList.append(entry[1:4] + \
1157                             (fmt % entry[4],
1158                              fmt % entry[5],
1159                              entry[0],
1160                              ", ".join(entry[6])) + (None,))
1161
1162    drawEnrichmentGraphPIL_tostream(termsList, header, fh, width, height)
1163
1164
1165def drawEnrichmentGraphPIL_tostream(termsList, headers, fh, width=None, height=None):
1166    from PIL import Image, ImageDraw, ImageFont
1167    backgroundColor = (255, 255, 255)
1168    textColor = (0, 0, 0)
1169    graphColor = (0, 0, 255)
1170    fontSize = height == None and 12 or (height - 60) / len(termsList)
1171    font = ImageFont.load_default()
1172    try:
1173        font = ImageFont.truetype("arial.ttf", fontSize)
1174    except:
1175        pass
1176    getMaxTextHeightHint = lambda l: max([font.getsize(t)[1] for t in l])
1177    getMaxTextWidthHint = lambda l: max([font.getsize(t)[0] for t in l])
1178    maxFoldWidth = width != None and min(150, width / 6) or 150
1179    maxFoldEnrichment = max([t[0] for t in termsList])
1180    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
1181    foldWidths = [int(foldNormalizationFactor * term[0]) for term in termsList]
1182    treeStep = 10
1183    treeWidth = {}
1184    for i, term in enumerate(termsList):
1185        treeWidth[i] = (term[-1] == None and 1 or treeWidth[term[-1]] + 1)
1186    treeStep = width != None and min(treeStep, width / (6 * max(treeWidth.values())) or 2) or treeStep
1187    treeWidth = [w * treeStep + foldWidths[i] for i, w in treeWidth.items()]
1188    treeWidth = max(treeWidth) - maxFoldWidth
1189    verticalMargin = 10
1190    horizontalMargin = 10
1191##    print verticalMargin, maxFoldWidth, treeWidth
1192##    treeWidth = 100
1193    firstColumnStart = verticalMargin + maxFoldWidth + treeWidth + 10
1194    secondColumnStart = firstColumnStart + getMaxTextWidthHint([str(t[1]) for t in termsList] + [headers[0]]) + 2
1195    thirdColumnStart = secondColumnStart + getMaxTextWidthHint([str(t[2]) for t in termsList] + [headers[1]]) + 2
1196    fourthColumnStart = thirdColumnStart + getMaxTextWidthHint([str(t[3]) for t in termsList] + [headers[2]]) + 2
1197    fifthColumnStart = fourthColumnStart + getMaxTextWidthHint([str(t[4]) for t in termsList] + [headers[3]]) + 4
1198##    maxAnnotationTextWidth = width==None and getMaxTextWidthHint([str(t[4]) for t in termsList]+["Annotation"]) or (width - fourthColumnStart - verticalMargin) * 2 / 3
1199    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]]))
1200    sixthColumnStart = fifthColumnStart + maxAnnotationTextWidth + 4
1201    maxGenesTextWidth = width == None and getMaxTextWidthHint([str(t[6]) for t in termsList] + [headers[5]]) or (width - fifthColumnStart - verticalMargin) / 3
1202
1203    legendHeight = font.getsize("1234567890")[1] * 2
1204    termHeight = font.getsize("A")[1]
1205##    print fourthColumnStart, maxAnnotationTextWidth, verticalMargin
1206    width = sixthColumnStart + maxGenesTextWidth + verticalMargin
1207    height = len(termsList) * termHeight + 2 * (legendHeight + horizontalMargin)
1208
1209    image = Image.new("RGB", (width, height), backgroundColor)
1210    draw = ImageDraw.Draw(image)
1211
1212    def truncText(text, maxWidth, append=""):
1213        # print getMaxTextWidthHint([text]), maxAnnotationTextWidth
1214        if getMaxTextWidthHint([text]) > maxWidth:
1215            while getMaxTextWidthHint([text + "..." + append]) > maxWidth and text:
1216                text = text[:-1]
1217            if text:
1218                text = text + "..." + append
1219            else:
1220                text = append
1221        return text
1222    currentY = horizontalMargin + legendHeight
1223    connectAtX = {}
1224    for i, term in enumerate(termsList):
1225        draw.line([(verticalMargin, currentY + termHeight / 2), (verticalMargin + foldWidths[i], currentY + termHeight / 2)], width=termHeight - 2, fill=graphColor)
1226        draw.text((firstColumnStart, currentY), str(term[1]), font=font, fill=textColor)
1227        draw.text((secondColumnStart, currentY), str(term[2]), font=font, fill=textColor)
1228        draw.text((thirdColumnStart, currentY), str(term[3]), font=font, fill=textColor)
1229        draw.text((fourthColumnStart, currentY), str(term[4]), font=font, fill=textColor)
1230##        annotText = width!=None and truncText(str(term[5]), maxAnnotationTextWidth, str(term[5])) or str(term[4])
1231        annotText = width != None and truncText(str(term[5]), maxAnnotationTextWidth)
1232        draw.text((fifthColumnStart, currentY), annotText, font=font, fill=textColor)
1233        genesText = width != None and truncText(str(term[6]), maxGenesTextWidth) or str(term[6])
1234        draw.text((sixthColumnStart, currentY), genesText, font=font, fill=textColor)
1235        lineEnd = term[-1] == None and firstColumnStart - 10 or connectAtX[term[-1]]
1236        draw.line([(verticalMargin + foldWidths[i] + 1, currentY + termHeight / 2), (lineEnd, currentY + termHeight / 2)], width=1, fill=textColor)
1237        if term[-1] != None:
1238            draw.line([(lineEnd, currentY + termHeight / 2), (lineEnd, currentY + termHeight / 2 - termHeight * (i - term[-1]))], width=1, fill=textColor)
1239        connectAtX[i] = lineEnd - treeStep
1240        currentY += termHeight
1241
1242    currentY = horizontalMargin
1243    draw.text((firstColumnStart, currentY), headers[0], font=font, fill=textColor)
1244    draw.text((secondColumnStart, currentY), headers[1], font=font, fill=textColor)
1245    draw.text((thirdColumnStart, currentY), headers[2], font=font, fill=textColor)
1246    draw.text((fourthColumnStart, currentY), headers[3], font=font, fill=textColor)
1247    draw.text((fifthColumnStart, currentY), headers[4], font=font, fill=textColor)
1248    draw.text((sixthColumnStart, currentY), headers[5], font=font, fill=textColor)
1249
1250    horizontalMargin = 0
1251    # draw.line([(verticalMargin, height - horizontalMargin - legendHeight), (verticalMargin + maxFoldWidth, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1252    draw.line([(verticalMargin, horizontalMargin + legendHeight), (verticalMargin + maxFoldWidth, horizontalMargin + legendHeight)], width=1, fill=textColor)
1253    maxLabelWidth = getMaxTextWidthHint([" " + str(i) for i in range(int(maxFoldEnrichment + 1))])
1254    numOfLegendLabels = max(int(maxFoldWidth / maxLabelWidth), 2)
1255    for i in range(numOfLegendLabels + 1):
1256        # draw.line([(verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight/2), (verticalMargin + i*maxFoldWidth/10, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
1257        # draw.text((verticalMargin + i*maxFoldWidth/10 - font.getsize(str(i))[0]/2, height - horizontalMargin - legendHeight/2), str(i), font=font, fill=textColor)
1258
1259        label = str(int(i * maxFoldEnrichment / numOfLegendLabels))
1260        draw.line([(verticalMargin + i * maxFoldWidth / numOfLegendLabels, horizontalMargin + legendHeight / 2), (verticalMargin + i * maxFoldWidth / numOfLegendLabels, horizontalMargin + legendHeight)], width=1, fill=textColor)
1261        draw.text((verticalMargin + i * maxFoldWidth / numOfLegendLabels - font.getsize(label)[0] / 2, horizontalMargin), label, font=font, fill=textColor)
1262
1263    image.save(fh)
1264
1265
1266def drawEnrichmentGraphPylab_tostream(termsList, headers, fh, width=None, height=None, show=True):
1267    from matplotlib import pyplot as plt
1268    from matplotlib.patches import Rectangle
1269
1270    maxFoldWidth = width != None and min(150, width / 6) or 150
1271    maxFoldEnrichment = max([t[0] for t in termsList])
1272    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
1273##    foldWidths = [int(foldNormalizationFactor*term[0]) for term in termsList]
1274    foldWidths = [term[0] for term in termsList]
1275    treeStep = maxFoldEnrichment * 0.05
1276    treeWidth = {}
1277
1278    for i, term in enumerate(termsList):
1279        treeWidth[i] = (term[-1] == None and treeStep or treeWidth[term[-1]] + treeStep)
1280    maxTreeWidth = max(treeWidth)
1281
1282    connectAt = {}
1283    cellText = []
1284    axes1 = plt.axes([0.1, 0.1, 0.2, 0.8])
1285    for i, line in enumerate(termsList):
1286        enrichment, n, m, p_val, fdr_val, name, genes, parent = line
1287        r = Rectangle((0, len(termsList) - i - 0.4), enrichment, 0.8)
1288        plt.gca().add_patch(r)
1289        plt.plot([enrichment, connectAt.get(parent, maxFoldEnrichment + maxTreeWidth)], [len(termsList) - i, len(termsList) - i], color="black")
1290        connectAt[i] = connectAt.get(parent, maxFoldEnrichment + maxTreeWidth) - treeStep
1291        if parent != None:
1292            plt.plot([connectAt.get(parent)] * 2, [len(termsList) - i, len(termsList) - parent], color="black")
1293        cellText.append((str(n), str(m), p_val, fdr_val, name, genes))
1294
1295##    from Orange.orng.orngClustering import TableTextLayout
1296##    text = TableTextLayout((maxFoldEnrichment*1.1, len(termsList)), cellText)
1297    from Orange.orng.orngClustering import TablePlot
1298    if True:
1299        axes2 = plt.axes([0.3, 0.1, 0.6, 0.8], sharey=axes1)
1300        axes2.set_axis_off()
1301        table = TablePlot((0, len(termsList)), axes=plt.gca())
1302        for i, line in enumerate(cellText):
1303            for j, text in enumerate(line):
1304                table.add_cell(i, j, width=len(text), height=1, text=text, loc="left", edgecolor="w", facecolor="w")
1305
1306        table.set_figure(plt.gcf())
1307        plt.gca().add_artist(table)
1308        plt.gca()._set_artist_props(table)
1309##    plt.text(3, 3, "\n".join(["\t".join(text) for text in cellText]))
1310
1311##    table = plt.table(cellText=cellText, colLabels=headers, loc="right")
1312##    table.set_transform(plt.gca().transData)
1313##
1314##    table.set_xy(20,20)
1315    plt.show()
1316
1317
1318class Taxonomy(object):
1319    """Maps NCBI taxonomy ids to corresponding GO organism codes
1320    """
1321    common_org_map = {"297284": "9913", "30523": "9913",  # Bos taurus
1322                      "5782": "352472", "44689": "352472", "366501": "352472",  # Dictyostelium discoideum
1323                      "83333": "562",  # Escherichia coli
1324                      "52545": "4530", "4532": "4530", "65489": "4530", "4533": "4530", "77588": "4530", "29689": "4530",
1325                      "4538": "4530", "40148": "4530", "29690": "4530", "110450": "4530", "4534": "4530", "83309": "4530",
1326                      "4528": "4530", "127571": "4530", "40149": "4530", "83307": "4530", "63629": "4530", "4536": "4530",
1327                      "4535": "4530", "4537": "4530", "65491": "4530", "83308": "4530", "4529": "4530", "4530": "4530",
1328                      "39946": "4530", "39947": "4530", "110451": "4530", "364100": "4530", "364099": "4530", "4539": "4530",
1329                      }
1330    code_map = {"3702": "tair",  # Arabidopsis thaliana
1331                "9913": "goa_cow",  # Bos taurus
1332                "6239": "wb",  # Caenorhabditis elegans
1333                "3055": None,  # Chlamydomonas reinhardtii
1334                "7955": "zfin",  # Danio rerio (zebrafish)
1335                "352472": "dictyBase",  # Dictyostelium discoideum
1336                "7227": "fb",  # Drosophila melanogaster
1337                "562": "ecocyc",  # Escherichia coli
1338                "11103": None,  # Hepatitis C virus
1339                "9606": "goa_human",  # Homo sapiens
1340                "10090": "mgi",  # Mus musculus
1341                "2104": None,  # Mycoplasma pneumoniae
1342                "4530": "gramene_oryza",  # Oryza sativa
1343                "5833": "GeneDB_Pfalciparum",  # Plasmodium falciparum
1344                "4754": None,  # Pneumocystis carinii
1345                "10116": "rgd",  # Rattus norvegicus
1346                "4932": "sgd",  # Saccharomyces cerevisiae
1347                "4896": "GeneDB_Spombe",  # Schizosaccharomyces pombe
1348                "31033": None,  # Takifugu rubripes
1349                "8355": None,  # Xenopus laevis
1350                "4577": None  # Zea mays
1351                }
1352    version = 1
1353    __shared_state = {"tax": None}
1354
1355    def __init__(self):
1356        self.__dict__ = self.__shared_state
1357        if not self.tax:
1358            path = orngServerFiles.localpath_download("GO", "taxonomy.pickle")
1359            if os.path.isfile(path):
1360                self.tax = cPickle.load(open(path, "rb"))
1361            else:
1362                orngServerFiles.download("GO", "taxonomy.pickle")
1363                self.tax = cPickle.load(open(path, "rb"))
1364
1365    def __getitem__(self, key):
1366        key = self.common_org_map.get(key, key)
1367        return self.code_map[key]
1368
1369    def keys(self):
1370        return list(set(self.common_org_map.keys() + self.code_map.keys()))
1371
1372
1373def from_taxid(id):
1374    """ Return a set of GO organism codes that correspond to NCBI taxonomy id
1375    """
1376    return Taxonomy()[id]
1377
1378
1379def to_taxid(db_code):
1380    """ Return a set of NCBI taxonomy ids from db_code GO organism annotations
1381    """
1382    r = [key for key, val in Taxonomy().code_map.items() if db_code == val]
1383    return set(r)
1384
1385
1386class __progressCallbackWrapper:
1387    def __init__(self, callback):
1388        self.callback = callback
1389
1390    def __call__(self, bCount, bSize, fSize):
1391        fSize = 10000000 if fSize == -1 else fSize
1392        self.callback(100 * bCount * bSize / fSize)
1393
1394from .obiGenomicsUpdate import Update as UpdateBase
1395
1396
1397class Update(UpdateBase):
1398    def __init__(self, local_database_path=None, progressCallback=None):
1399        UpdateBase.__init__(self, local_database_path or getDataDir(), progressCallback)
1400
1401    def CheckModified(self, addr, date=None):
1402        return date < self.GetLastModified(addr) if date else True
1403
1404    def CheckModifiedOrg(self, org):
1405        return self.CheckModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz", self.LastModifiedOrg(org))
1406
1407    def LastModifiedOrg(self, org):
1408        return self.shelve.get((Update.UpdateAnnotation, (org,)), None)
1409
1410    def GetLastModified(self, addr):
1411        stream = urllib2.urlopen(addr)
1412        return datetime.strptime(stream.headers.get("Last-Modified"), "%a, %d %b %Y %H:%M:%S %Z")
1413##        return stream.headers.get("Last-Modified")
1414
1415    def GetAvailableOrganisms(self):
1416        source = urllib2.urlopen("http://www.geneontology.org/gene-associations/").read()
1417        return [s.split(".")[1] for s in sorted(set(re.findall("gene_association\.[a-zA-z0-9_]+?\.gz", source)))]
1418
1419    def GetDownloadedOrganisms(self):
1420        return [name.split(".")[1] for name in os.listdir(self.local_database_path) if name.startswith("gene_association")]
1421
1422    def IsUpdatable(self, func, args):
1423        if func == Update.UpdateOntology:
1424            return self.CheckModified("http://www.geneontology.org/ontology/gene_ontology.obo", self.shelve.get((Update.UpdateOntology, ()), None))
1425        elif func == Update.UpdateAnnotation:
1426            return self.CheckModifiedOrg(args[0])
1427
1428    def GetDownloadable(self):
1429        orgs = set(self.GetAvailableOrganisms()) - set(self.GetDownloadedOrganisms())
1430        ret = []
1431        if (Update.UpdateOntology, ()) not in self.shelve:
1432            ret.append((Update.UpdateOntology, ()))
1433        if orgs:
1434            ret.extend([(Update.UpdateAnnotation, (org,)) for org in orgs])
1435        return ret
1436
1437    def UpdateOntology(self):
1438        Ontology.DownloadOntology(os.path.join(self.local_database_path, "gene_ontology_edit.obo.tar.gz"), self.progressCallback)
1439        self._update(Update.UpdateOntology, (), self.GetLastModified("http://www.geneontology.org/ontology/gene_ontology.obo"))
1440
1441    def UpdateAnnotation(self, org):
1442        Annotations.DownloadAnnotations(org, os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"), self.progressCallback)
1443        self._update(Update.UpdateAnnotation, (org,), self.GetLastModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz"))
1444
1445    def UpdateTaxonomy(self, org):
1446        exclude = ["goa_uniprot", "goa_pdb", "GeneDB_tsetse", "reactome", "goa_zebrafish", "goa_rat", "goa_mouse"]
1447
1448        orgs = self.GetAvailableOrganisms()
1449        tax = defaultdict(set)
1450
1451        for org in orgs:
1452            if org in exclude:
1453                continue
1454            try:
1455                a = obiGO.Annotations(os.path.join(self.local_database_path, "gene_association." + org + ".tar.gz"))
1456                taxons = set(ann.taxon for ann in a.annotations)
1457                for taxId in [t.split(":")[-1] for t in taxons if "|" not in t]:  # exclude taxons with cardinality 2
1458                    tax[taxId].add(org)
1459            except Exception, ex:
1460                print ex
1461
1462        cPickle.dump(dict(tax), open(os.path.join(path, "taxonomy.pickle"), "wb"))
1463
1464
1465def _test1():
1466##    Ontology.DownloadOntology("ontology_arch.tar.gz")
1467##    Annotations.DownloadAnnotations("sgd", "annotations_arch.tar.gz")
1468    def _print(f):
1469        print f
1470    o = Ontology("ontology_arch.tar.gz")
1471    a = Annotations("annotations_arch.tar.gz", ontology=o)
1472
1473    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
1474##    profile.runctx("a.GetEnrichedTerms(sorted(a.geneNames)[:100])", {"a":a}, {})
1475    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
1476    d1 = a.GetEnrichedTerms(sorted(a.geneNames)[:1000])  # , progressCallback=_print)
1477
1478##    print a.GetEnrichedTerms(sorted(a.geneNames)[:100])#, progressCallback=_print)
1479
1480
1481def _test2():
1482    o = Ontology()
1483    a = Annotations("human", ontology=o)
1484    clusterGenes = sorted(a.geneNames)[:100]
1485    for i in range(10):
1486        terms = a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["P"])
1487        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["C"])
1488        a.GetEnrichedTerms(sorted(a.geneNames)[:200], aspect=["F"])
1489        print i
1490#    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1491
1492##    drawEnrichmentGraph([("bal", 1.0, 5, 6, 0.1, 0.4, ["vv"]),
1493##                        ("GO:0019079", 0.5, 5, 6, 0.1, 0.4, ["cc", "bb"]),
1494##                        ("GO:0022415", 0.4, 5, 7, 0.11, 0.4, ["cc1", "bb"])], open("graph.png", "wb"), None, None)
1495
1496
1497def _test3():
1498    o = Ontology()
1499    a = Annotations("sgd", ontology=o)
1500##    a = Annotations(list(a)[3:len(a)/3], ontology=o)
1501    clusterGenes = sorted(a.geneNames)[:1] + sorted(a.geneNames)[-1:]
1502##    clusterGenes = [g + "_" + str(i%5) for g in sorted(a.geneNames)[:2]]
1503    exonMap = dict([(gene, [gene + "_E%i" % i for i in range(10)]) for gene in a.geneNames])
1504    a.RemapGenes(exonMap)
1505##    o.reverseAliasMapper = o.aliasMapper = {}
1506    terms = a.GetEnrichedTerms(exonMap.values()[0][:2] + exonMap.values()[-1][2:])
1507##    terms = a.GetEnrichedTerms(clusterGenes)
1508    print terms
1509##    a.DrawEnrichmentGraph(filterByPValue(terms), len(clusterGenes), len(a.geneNames))
1510    a.DrawEnrichmentGraph(filterByPValue(terms, maxPValue=0.1), len(clusterGenes), len(a.geneNames))
1511
1512if __name__ == "__main__":
1513    _test2()
Note: See TracBrowser for help on using the repository browser.