source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1834:91f840200ae3

Revision 1834:91f840200ae3, 60.7 KB checked in by Flashpoint <vid.flashpoint@…>, 10 months ago (diff)

Changed the exception value from KEGG.OrganismNotFound to obiTaxonomy.UnknownSpeciesIdentifier

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