source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1771:bf833fa5d2a3

Revision 1771:bf833fa5d2a3, 60.6 KB checked in by Ales Erjavec <ales.erjavec@…>, 12 months ago (diff)

Fixed 'Ontology' and 'Annotations' stream parsing.

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