source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1768:9cb38945684f

Revision 1768:9cb38945684f, 59.3 KB checked in by Ales Erjavec <ales.erjavec@…>, 12 months ago (diff)

Removed deprecated and nonfunctioning 'drawEnrichmentGraph_tostream' function.

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