source: orange-bioinformatics/orangecontrib/bio/obiGO.py @ 1889:2ff9da5aed53

Revision 1889:2ff9da5aed53, 60.7 KB checked in by markotoplak, 6 months ago (diff)

Attributes m and n in Hypergeometric test can now be swapped and results do not change. The implementation is more resistant to rounding errors.

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