source: orange-bioinformatics/orangecontrib/bio/obiGO.py @ 1873:0810c5708cc5

Revision 1873:0810c5708cc5, 60.7 KB checked in by Ales Erjavec <ales.erjavec@…>, 7 months ago (diff)

Moved '_bioinformatics' into orangecontrib namespace.

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
[1834]27from . import obiProb, obiGene, obiTaxonomy
[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()
[1771]287
288        if file is not None:
[538]289            self.ParseFile(file, progressCallback)
[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:
[1771]304            # Load the default ontology file.
[639]305            fool = self.Load(progressCallback)
[1767]306            # A fool and his attributes are soon parted
307            self.__dict__ = fool.__dict__
[538]308
309    @classmethod
310    def Load(cls, progressCallback=None):
[1334]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.
[1771]314
[538]315        """
[1767]316        filename = os.path.join(default_database_path,
317                                "gene_ontology_edit.obo.tar.gz")
[588]318        if not os.path.isfile(filename) and not os.path.isdir(filename):
[571]319            orngServerFiles.download("GO", "gene_ontology_edit.obo.tar.gz")
[1771]320
321        return cls(filename, progressCallback=progressCallback)
[1767]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        """
[1771]328        if isinstance(file, basestring):
[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)
[1771]333            elif os.path.isdir(file):
334                f = open(os.path.join(file, "gene_ontology_edit.obo"))
[588]335            else:
[1771]336                raise ValueError("Cannot open %r for parsing" % file)
[538]337        else:
338            f = file
[1767]339
[538]340        data = f.readlines()
341        data = "".join([line for line in data if not line.startswith("!")])
[571]342        self.header = data[: data.index("[Term]")]
[1767]343        c = re.compile("\[.+?\].*?\n\n", re.DOTALL)
344        data = c.findall(data)
[639]345
[1015]346        milestones = orngMisc.progressBarMilestones(len(data), 90)
[538]347        for i, block in enumerate(builtinOBOObjects + data):
348            if block.startswith("[Term]"):
[639]349                term = Term(block, self)
[538]350                self.terms[term.id] = term
351            elif block.startswith("[Typedef]"):
[639]352                typedef = Typedef(block, self)
[538]353                self.typedefs[typedef.id] = typedef
354            elif block.startswith("[Instance]"):
[639]355                instance = Instance(block, self)
[538]356                self.instances[instance.id] = instance
[542]357            if progressCallback and i in milestones:
[1767]358                progressCallback(90.0 * i / len(data))
359
[538]360        self.aliasMapper = {}
[964]361        self.reverseAliasMapper = defaultdict(set)
[1015]362        milestones = orngMisc.progressBarMilestones(len(self.terms), 10)
[964]363        for i, (id, term) in enumerate(self.terms.iteritems()):
[538]364            for typeId, parent in term.related:
365                self.terms[parent].relatedTo.add((typeId, id))
[676]366            try:
[1767]367                self.aliasMapper.update([(alt_id, id)
368                                         for alt_id in term.alt_id])
[964]369                self.reverseAliasMapper[id].union_update(term.alt_id)
370            except AttributeError:
[676]371                pass
[964]372            if progressCallback and i in milestones:
[1767]373                progressCallback(90.0 + 10.0 * i / len(self.terms))
[538]374
[571]375    def GetDefinedSlimsSubsets(self):
[1334]376        """ Return a list of defined subsets.
[571]377        """
[1767]378        return [line.split()[1] for line in self.header.splitlines()
379                if line.startswith("subsetdef:")]
[571]380
[567]381    def SetSlimsSubset(self, subset):
[1334]382        """ Set the slims term subset to subset. If subset is a string it
383        must equal one of the defined subsetdef.
[567]384        """
385        if type(subset) == str:
[1767]386            self.slimsSubset = [id for id, term in self.terms.items()
387                                if subset in getattr(term, "subset", set())]
[567]388        else:
389            self.slimsSubset = set(subset)
[578]390
[994]391    def GetSlimsSubset(self, subset):
[1767]392        return [id for id, term in self.terms.items()
393                if subset in getattr(term, "subset", set())]
[994]394
[578]395    def GetSlimTerms(self, termId):
[639]396        """ Return a list of slim terms for termId.
[578]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:
[1767]407                queue.update(set(id for typeId, id in self[term].related) -
408                             visited)
[578]409        return slims
[567]410
[538]411    def ExtractSuperGraph(self, terms):
[639]412        """ Return all super terms of terms up to the most general one.
[538]413        """
[647]414        terms = [terms] if type(terms) == str else terms
[538]415        visited = set()
416        queue = set(terms)
417        while queue:
418            term = queue.pop()
419            visited.add(term)
[1767]420            queue.update(set(id for typeId, id in self[term].related) -
421                         visited)
[538]422        return visited
423
424    def ExtractSubGraph(self, terms):
[639]425        """ Return all sub terms of terms.
[538]426        """
[647]427        terms = [terms] if type(terms) == str else terms
[538]428        visited = set()
429        queue = set(terms)
430        while queue:
431            term = queue.pop()
432            visited.add(term)
[1767]433            queue.update(set(id for typeId, id in self[term].relatedTo) -
434                         visited)
[538]435        return visited
436
437    def GetTermDepth(self, term, cache_={}):
[1334]438        """ Return the minimum depth of a term (length of the shortest
439        path to this term from the top level term).
[565]440        """
[1025]441        if term not in cache_:
[1767]442            cache_[term] = min([self.GetTermDepth(parent) + 1
443                                for typeId, parent in self[term].related] or
444                               [1])
[1025]445        return cache_[term]
[538]446
[639]447    def __getitem__(self, id):
448        """
[1770]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)
[538]457
[619]458    def __iter__(self):
[1770]459        """
460        Iterate over all ids in ontology.
[639]461        """
[625]462        return iter(self.terms)
463
464    def __len__(self):
[1770]465        """
466        Return number of objects in ontology.
[639]467        """
[625]468        return len(self.terms)
[619]469
[640]470    def __contains__(self, id):
[1770]471        """
472        Return `True` if a term with `id` is present in the ontology.
473        """
[791]474        return id in self.terms or id in self.aliasMapper
[640]475
[538]476    @staticmethod
477    def DownloadOntology(file, progressCallback=None):
[1767]478        tFile = tarfile.open(file, "w:gz") if isinstance(file, basestring) \
479                else file
[538]480        tmpDir = os.path.join(orngEnviron.bufferDir, "tmp_go/")
481        try:
482            os.mkdir(tmpDir)
483        except Exception:
484            pass
[1767]485        urlretrieve("http://www.geneontology.org/ontology/gene_ontology_edit.obo",
[1334]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")
[542]490        tFile.close()
491        os.remove(os.path.join(tmpDir, "gene_ontology_edit.obo"))
[538]492
[1609]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:
[1767]498            filename = os.path.join(default_database_path,
499                                    "gene_ontology_edit@rev%s.obo" % rev)
[1609]500        r = urllib2.urlopen(url)
[1767]501
[1610]502        with open(filename + ".part", "wb") as f:
503            shutil.copyfileobj(r, f)
[1767]504
[1610]505        os.rename(filename + ".part", filename)
[1767]506
507
[570]508_re_obj_name_ = re.compile("([a-zA-z0-9-_]+)")
509
[1767]510
[570]511class AnnotationRecord(object):
[1334]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)
[1767]521
[570]522    """
[1334]523    __slots__ = annotationFields + ["geneName", "GOId", "evidence",
524                                    "aspect", "alias", "additionalAliases"]
[1767]525
[570]526    def __init__(self, fullText):
[1351]527        """\
528        :param fulText: A single line from the annotation file.
[1767]529
[1351]530        """
[1256]531        for slot, val in zip(annotationFields, fullText.split("\t")):
[1351]532            setattr(self, slot, intern(val))
[1334]533        self.geneName = self.DB_Object_Symbol
534        self.GOId = self.GO_ID
535        self.evidence = self.Evidence_Code
536        self.aspect = self.Aspect
[1351]537        self.alias = list(map(intern, self.DB_Object_Synonym.split("|")))
[570]538
[882]539        self.additionalAliases = []
[570]540        if ":" in self.DB_Object_Name:
[1767]541            self.additionalAliases = []  # _re_obj_name_.findall(self.DB_Object_Name.split(":")[0])
[570]542
543    def __getattr__(self, name):
544        if name in annotationFieldsDict:
[1256]545            return getattr(self, self.__slots__[annotationFieldsDict[name]])
[570]546        else:
547            raise AttributeError(name)
548
549
[538]550class Annotations(object):
[565]551    """Annotations object holds the annotations.
552    """
[928]553    version = 2
[1767]554
555    def __init__(self, file=None, ontology=None, genematcher=None,
556                 progressCallback=None, rev=None):
[1334]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.
[1767]561
[565]562        """
[538]563        self.ontology = ontology
564        self.allAnnotations = defaultdict(list)
565        self.geneAnnotations = defaultdict(list)
566        self.termAnnotations = defaultdict(list)
[791]567        self._geneNames = None
568        self._geneNamesDict = None
569        self._aliasMapper = None
[538]570        self.additionalAliases = {}
571        self.annotations = []
572        self.header = ""
[791]573        self.genematcher = genematcher
574        self.taxid = None
[1771]575
[679]576        if type(file) in [list, set, dict, Annotations]:
577            for ann in file:
578                self.AddAnnotation(ann)
[791]579            if type(file, Annotations):
[1771]580                self.taxid = file.taxid
581
[1609]582        elif isinstance(file, basestring) and os.path.exists(file):
[538]583            self.ParseFile(file, progressCallback)
[791]584            try:
585                self.taxid = to_taxid(os.path.basename(file).split(".")[1]).pop()
586            except IOError:
587                pass
[1771]588
589        elif isinstance(file, basestring):
590            # Assuming organism code/name
[1609]591            if rev is not None:
592                if not _CVS_REVISION_RE.match(rev):
593                    raise ValueError("Invalid revision format")
[1767]594
[1609]595                if rev.startswith("rev"):
596                    rev = rev[3:]
597                code = self.organism_name_search(file)
598                filename = os.path.join(default_database_path,
[1767]599                                        "gene_association.%s@rev%s.tar.gz" %
600                                        (code, rev))
601
[1609]602                if not os.path.exists(filename):
[1771]603                    self.DownloadAnnotationsAtRev(
604                        code, rev, filename, progressCallback)
[1767]605
[1609]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()
[1771]612        elif file is not None:
613            self.ParseFile(file, progressCallback)
614
[791]615        if not self.genematcher and self.taxid:
[1771]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
[791]625        if self.genematcher:
626            self.genematcher.set_targets(self.geneNames)
[1767]627
[745]628    @classmethod
629    def organism_name_search(cls, org):
[1767]630        ids = to_taxid(org)
[1632]631        from . import obiTaxonomy as tax
[745]632        if not ids:
[1767]633            ids = [org] if org in Taxonomy().common_org_map.keys() + \
634                  Taxonomy().code_map.keys() else []
[928]635        if not ids:
636            ids = tax.to_taxid(org, mapTo=Taxonomy().keys())
[745]637        if not ids:
638            ids = tax.search(org, exact=True)
[928]639            ids = set(ids).intersection(Taxonomy().keys())
[745]640        if not ids:
641            ids = tax.search(org)
[928]642            ids = set(ids).intersection(Taxonomy().keys())
643        codes = set([from_taxid(id) for id in ids])
[745]644        if len(codes) > 1:
[1767]645            raise tax.MultipleSpeciesException(
646                ", ".join(["%s: %s" % (str(from_taxid(id)), tax.name(id))
647                           for id in ids]))
[745]648        elif len(codes) == 0:
[1767]649            raise tax.UnknownSpeciesIdentifier(org)
[745]650        return codes.pop()
651
[804]652    @classmethod
653    def organism_version(cls, name):
654        name = organism_name_search(name)
[1767]655        orngServerFiles.localpath_download(
656            "GO", "gene_association.%s.tar.gz" % name)
[1334]657        return ("v%i." % cls.version) + orngServerFiles.info("GO",
658                        "gene_association.%s.tar.gz" % name)["datetime"]
[804]659
[565]660    def SetOntology(self, ontology):
[1334]661        """ Set the ontology to use in the annotations mapping.
662        """
[565]663        self.allAnnotations = defaultdict(list)
664        self._ontology = ontology
665
666    def GetOntology(self):
667        return self._ontology
668
[1334]669    ontology = property(GetOntology, SetOntology,
670                        doc="Ontology object for annotations")
[1767]671
[538]672    @classmethod
[791]673    def Load(cls, org, ontology=None, genematcher=None, progressCallback=None):
[1334]674        """A class method that tries to load the association file for the
675        given organism from default_database_path.
[538]676        """
[791]677        code = organism_name_search(org)
[1767]678
[1834]679        print "CODE: %s" % code
680        print "ORG: %s" % org
681
[791]682        file = "gene_association.%s.tar.gz" % code
[612]683
[791]684        path = os.path.join(orngServerFiles.localpath("GO"), file)
[1712]685
[612]686        if not os.path.exists(path):
[1712]687            sf = orngServerFiles.ServerFiles()
688            available = sf.listfiles("GO")
689            if file not in available:
[1716]690                from . import obiKEGG
[1834]691                raise obiTaxonomy.UnknownSpeciesIdentifier(org + str(code))
[791]692            orngServerFiles.download("GO", file)
[1712]693
[1767]694        return cls(path, ontology=ontology, genematcher=genematcher,
695                   progressCallback=progressCallback)
696
[538]697    def ParseFile(self, file, progressCallback=None):
[1334]698        """ Parse and load the annotations from file. Report progress
699        with progressCallback.
[639]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
[565]705        """
[1771]706        if isinstance(file, basestring):
[588]707            if os.path.isfile(file) and tarfile.is_tarfile(file):
708                f = tarfile.open(file).extractfile("gene_association")
[1609]709            elif os.path.isfile(file) and file.endswith(".gz"):
[1767]710                f = gzip.open(file)
[588]711            elif os.path.isfile(file):
712                f = open(file)
[1771]713            elif os.path.isdir(file):
714                f = open(os.path.join(file, "gene_association"))
[588]715            else:
[1771]716                raise ValueError("Cannot open %r for parsing." % file)
[538]717        else:
718            f = file
[1682]719        lines = [line for line in f.read().splitlines() if line.strip()]
[1767]720
[1015]721        milestones = orngMisc.progressBarMilestones(len(lines), 100)
[1767]722        for i, line in enumerate(lines):
[538]723            if line.startswith("!"):
724                self.header = self.header + line + "\n"
725                continue
[1767]726
727            a = AnnotationRecord(line)
[679]728            self.AddAnnotation(a)
[964]729#            self.annotations.append(a)
[538]730            if progressCallback and i in milestones:
[1767]731                progressCallback(100.0 * i / len(lines))
[538]732
[679]733    def AddAnnotation(self, a):
[1334]734        """ Add a single `AnotationRecord` instance to this `Annotations`
735        object.
736        """
[679]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
[1767]741
[791]742        self.geneAnnotations[a.geneName].append(a)
[679]743        self.annotations.append(a)
744        self.termAnnotations[a.GOId].append(a)
745        self.allAnnotations = defaultdict(list)
[1767]746
[791]747        self._geneNames = None
748        self._geneNamesDict = None
749        self._aliasMapper = None
750
751    @property
752    def geneNamesDict(self):
[1334]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)
[791]757        return self._geneNamesDict
758
759    @property
760    def geneNames(self):
[1334]761        if getattr(self, "_geneNames", None) is None:
[791]762            self._geneNames = set([ann.geneName for ann in self.annotations])
763        return self._geneNames
764
765    @property
766    def aliasMapper(self):
[1334]767        if getattr(self, "_aliasMapper", None) is None:
768            self._aliasMapper = {}
769            for ann in self.annotations:
[1767]770                self._aliasMapper.update([(alias, ann.geneName)
771                                          for alias in ann.alias +
772                                           [ann.geneName, ann.DB_Object_ID]])
[791]773        return self._aliasMapper
[679]774
[539]775    def GetGeneNamesTranslator(self, genes):
[1334]776        """ Return a dictionary mapping canonical names (DB_Object_Symbol)
777        to `genes`.
[1767]778
[1334]779        """
[539]780        def alias(gene):
[791]781            if self.genematcher:
782                return self.genematcher.umatch(gene)
783            else:
[1334]784                return gene if gene in self.geneNames else \
785                        self.aliasMapper.get(gene,
786                             self.additionalAliases.get(gene, None))
[1767]787
[539]788        return dict([(alias(gene), gene) for gene in genes if alias(gene)])
789
[803]790    def _CollectAnnotations(self, id, visited):
[639]791        """ Recursive function collects and caches all annotations for id
792        """
[803]793        if id not in self.allAnnotations and id not in visited:
[677]794            if id in self.ontology.reverseAliasMapper:
[1767]795                annotations = [self.termAnnotations.get(alt_id, [])
796                               for alt_id in
797                               self.ontology.reverseAliasMapper[id]] + \
798                              [self.termAnnotations[id]]
[677]799            else:
[1334]800                ## annotations for this term alone
[1767]801                annotations = [self.termAnnotations[id]]
[803]802            visited.add(id)
[538]803            for typeId, child in self.ontology[id].relatedTo:
[803]804                aa = self._CollectAnnotations(child, visited)
[1334]805                if type(aa) == set:
806                    ## if it was already reduced in GetAllAnnotations
[538]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):
[1334]814        """ Return a set of all annotations (instances if `AnnotationRectord`)
815        for GO term `id` and all it's subterms.
[1767]816
[1334]817        :param id: GO term id
818        :type id: str
[1767]819
[538]820        """
[803]821        visited = set()
[677]822        id = self.ontology.aliasMapper.get(id, id)
[1767]823        if id not in self.allAnnotations or \
824                type(self.allAnnotations[id]) == list:
[538]825            annot_set = set()
[803]826            for annots in self._CollectAnnotations(id, set()):
[538]827                annot_set.update(annots)
828            self.allAnnotations[id] = annot_set
829        return self.allAnnotations[id]
830
[1767]831    def GetAllGenes(self, id, evidenceCodes=None):
[1334]832        """ Return a list of genes annotated by specified `evidenceCodes`
833        to GO term 'id' and all it's subterms."
[1767]834
[1334]835        :param id: GO term id
836        :type id: str
[1767]837
[1334]838        :param evidneceCodes: List of evidence codes to consider when
839                              matching annotations to terms.
840        :type evidenceCodes: list-of-strings
[538]841        """
842        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
843        annotations = self.GetAllAnnotations(id)
[1767]844        return list(set([ann.geneName for ann in annotations
845                         if ann.Evidence_Code in evidenceCodes]))
[538]846
[1767]847    def GetEnrichedTerms(self, genes, reference=None, evidenceCodes=None,
[1422]848                         slimsOnly=False, aspect=None, prob=obiProb.Binomial(),
[1334]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).
[1767]853
[1334]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
[1422]859        :param aspect: Which aspects to use. Use all by default. "P", "F", "C"
860            or a set containing these elements.
[538]861        """
[539]862        revGenesDict = self.GetGeneNamesTranslator(genes)
[538]863        genes = set(revGenesDict.keys())
[697]864        if reference:
865            refGenesDict = self.GetGeneNamesTranslator(reference)
866            reference = set(refGenesDict.keys())
867        else:
868            reference = self.geneNames
[1422]869
870        if aspect == None:
871            aspects_set = set(["P", "C", "F"])
[1767]872        elif isinstance(aspect, basestring):
873            aspects_set = set([aspect])
[1422]874        else:
[1767]875            aspects_set = aspect
[1422]876
[538]877        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
[1767]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
[538]890        annotationsDict = defaultdict(set)
891        for ann in annotations:
892            annotationsDict[ann.GO_ID].add(ann)
[1767]893
[1149]894        if slimsOnly and not self.ontology.slimsSubset:
[1767]895            warnings.warn("Unspecified slims subset in the ontology! "
896                          "Using 'goslim_generic' subset", UserWarning)
[1149]897            self.ontology.SetSlimsSubset("goslim_generic")
[1767]898
[1769]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)
[538]909        res = {}
[1767]910
[1015]911        milestones = orngMisc.progressBarMilestones(len(terms), 100)
[538]912        for i, term in enumerate(terms):
[567]913            if slimsOnly and term not in self.ontology.slimsSubset:
914                continue
[679]915            allAnnotations = self.GetAllAnnotations(term).intersection(refAnnotations)
916##            allAnnotations.intersection_update(refAnnotations)
[619]917            allAnnotatedGenes = set([ann.geneName for ann in allAnnotations])
[679]918            mappedGenes = genes.intersection(allAnnotatedGenes)
[726]919##            if not mappedGenes:
920##                print >> sys.stderr, term, sorted(genes)
921##                print >> sys.stderr, sorted(allAnnotatedGenes)
922##                return
[538]923            if len(reference) > len(allAnnotatedGenes):
924                mappedReferenceGenes = reference.intersection(allAnnotatedGenes)
925            else:
926                mappedReferenceGenes = allAnnotatedGenes.intersection(reference)
[1334]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))
[538]931            if progressCallback and i in milestones:
932                progressCallback(100.0 * i / len(terms))
[1098]933        if useFDR:
[1767]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]))])
[538]938        return res
939
[1767]940    def GetAnnotatedTerms(self, genes, directAnnotationOnly=False,
941                          evidenceCodes=None, progressCallback=None):
[639]942        """ Return all terms that are annotated by genes with evidenceCodes.
[538]943        """
[647]944        genes = [genes] if type(genes) == str else genes
[539]945        revGenesDict = self.GetGeneNamesTranslator(genes)
[578]946        genes = set(revGenesDict.keys())
[538]947        evidenceCodes = set(evidenceCodes or evidenceDict.keys())
[1767]948        annotations = [ann for gene in genes for ann in self.geneAnnotations[gene]
[1334]949                       if ann.Evidence_Code in evidenceCodes]
[538]950        dd = defaultdict(set)
951        for ann in annotations:
[860]952            dd[ann.GO_ID].add(revGenesDict.get(ann.geneName, ann.geneName))
[578]953        if not directAnnotationOnly:
[1769]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)
[538]964            for i, term in enumerate(terms):
[679]965                termAnnots = self.GetAllAnnotations(term).intersection(annotations)
966##                termAnnots.intersection_update(annotations)
[1767]967                dd[term].update([revGenesDict.get(ann.geneName, ann.geneName)
968                                 for ann in termAnnots])
[578]969        return dict(dd)
[619]970
[1767]971    def DrawEnrichmentGraph(self, terms, clusterSize, refSize=None,
972                            file="graph.png", width=None, height=None,
973                            precison=3):
[666]974        refSize = len(self.geneNames) if refSize == None else refSize
[1767]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)
[664]990
[679]991    def __add__(self, iterable):
992        """ Return a new Annotations object with combined annotations
993        """
[1767]994        return Annotations([a for a in self] + [a for a in iterable],
995                           ontology=self.ontology)
[679]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
[1767]1005
[619]1006    def __iter__(self):
[639]1007        """ Iterate over all AnnotationRecord objects in annotations
1008        """
[619]1009        return iter(self.annotations)
1010
1011    def __len__(self):
[639]1012        """ Return the number of annotations
1013        """
[619]1014        return len(self.annotations)
1015
1016    def __getitem__(self, index):
[639]1017        """ Return the i-th annotation record
1018        """
[619]1019        return self.annotations[index]
[679]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        """
[1767]1042        """
[679]1043        for gene in map:
1044            annotations = self.geneAnnotations[gene]
1045            for ann in annotations:
1046                for name in map[gene]:
[1767]1047                    ann1 = copy.copy(ann)
[679]1048                    ann1.geneName = name
1049                    self.add(ann1)
[970]1050        self.genematcher = obiGene.GMDirect()
[971]1051        try:
1052            del self._geneNames
1053        except Exception:
1054            pass
1055        self.genematcher.set_targets(self.geneNames)
[1767]1056
[538]1057    @staticmethod
1058    def DownloadAnnotations(org, file, progressCallback=None):
[1771]1059        if isinstance(file, basestring):
1060            tFile = tarfile.open(file, "w:gz")
1061        else:
1062            tFile = file
1063
[538]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"
[1767]1070        urlretrieve(("http://www.geneontology.org/gene-associations/" +
1071                     fileName),
[1334]1072                    os.path.join(tmpDir, fileName),
[1767]1073                    progressCallback and __progressCallbackWrapper(progressCallback))
[538]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()
[1767]1079
1080        tFile.add(os.path.join(tmpDir, "gene_association." + org),
1081                  "gene_association")
[1334]1082        annotation = Annotations(os.path.join(tmpDir, "gene_association." + org),
1083                    genematcher=obiGene.GMDirect(), progressCallback=progressCallback)
[633]1084        cPickle.dump(annotation.geneNames, open(os.path.join(tmpDir, "gene_names.pickle"), "wb"))
[538]1085        tFile.add(os.path.join(tmpDir, "gene_names.pickle"), "gene_names.pickle")
1086        tFile.close()
[542]1087        os.remove(os.path.join(tmpDir, "gene_association." + org))
1088        os.remove(os.path.join(tmpDir, "gene_names.pickle"))
[1767]1089
[1609]1090    @staticmethod
[1767]1091    def DownloadAnnotationsAtRev(org, rev, filename=None,
1092                                 progressCallback=None):
[1609]1093        if filename is None:
1094            filename = os.path.join(default_database_path,
[1767]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))
[1609]1099        url += ";content-type=application%2Fx-gzip"
1100        r = urllib2.urlopen(url)
[1767]1101
[1610]1102        with open(filename + ".part", "wb") as f:
1103            shutil.copyfileobj(r, f)
[1767]1104
[1610]1105        os.rename(filename + ".part", filename)
[538]1106
[1643]1107from .obiTaxonomy import pickled_cache
[753]1108
[1767]1109
1110@pickled_cache(None, [("GO", "taxonomy.pickle"),
1111                      ("Taxonomy", "ncbi_taxonomy.tar.gz")])
[753]1112def organism_name_search(name):
1113    return Annotations.organism_name_search(name)
1114
[1767]1115
[662]1116def filterByPValue(terms, maxPValue=0.1):
[1334]1117    """ Filters the terms by the p-value. Asumes terms is a dict with
1118    the same structure as returned from GetEnrichedTerms.
[1767]1119
[662]1120    """
[1767]1121    return dict(filter(lambda (k, e): e[1] <= maxPValue, terms.items()))
1122
[662]1123
1124def filterByFrequency(terms, minF=2):
[1334]1125    """ Filters the terms by the cluster frequency. Asumes terms is
1126    a dict with the same structure as returned from GetEnrichedTerms.
[1767]1127
[662]1128    """
[1767]1129    return dict(filter(lambda (k, e): len(e[0]) >= minF, terms.items()))
1130
[662]1131
1132def filterByRefFrequency(terms, minF=4):
[1334]1133    """ Filters the terms by the reference frequency. Asumes terms is
1134    a dict with the same structure as returned from GetEnrichedTerms.
[1767]1135
[662]1136    """
[1767]1137    return dict(filter(lambda (k, e): e[2] >= minF, terms.items()))
[662]1138
1139
[1768]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)
[1767]1145
1146
[1768]1147def drawEnrichmentGraph_tostreamMk2(enriched, fh, width, height, header=None,
1148                                    ontology=None, precison=4):
[664]1149    ontology = ontology if ontology else Ontology()
[1767]1150    header = header if header else ["List", "Total", "p-value", "FDR",
1151                                    "Names", "Genes"]
[664]1152    GOTerms = dict([(t[0], t) for t in enriched if t[0] in ontology])
[1767]1153
[664]1154    def getParents(term):
1155        parents = ontology.ExtractSuperGraph([term])
1156        parents = [id for id in parents if id in GOTerms and id != term]
[1767]1157        c = reduce(set.union, [set(ontology.ExtractSuperGraph([id])) -
1158                               set([id]) for id in parents], set())
[664]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])
[1767]1162    # print "Parentes", parents
1163
[664]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]]
[1767]1167    # print "Top level terms", topLevelTerms
1168    termsList = []
[664]1169    fmt = "%" + ".%if" % precison
[1767]1170
[664]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,))
[1767]1177        parent = len(termsList) - 1
[664]1178        for c in getChildren(term):
1179            collect(c, parent)
[1767]1180
[664]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
[797]1191    drawEnrichmentGraphPIL_tostream(termsList, header, fh, width, height)
[1767]1192
1193
[664]1194def drawEnrichmentGraphPIL_tostream(termsList, headers, fh, width=None, height=None):
[662]1195    from PIL import Image, ImageDraw, ImageFont
1196    backgroundColor = (255, 255, 255)
1197    textColor = (0, 0, 0)
1198    graphColor = (0, 0, 255)
[1767]1199    fontSize = height == None and 12 or (height - 60) / len(termsList)
[662]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])
[1767]1207    maxFoldWidth = width != None and min(150, width / 6) or 150
[662]1208    maxFoldEnrichment = max([t[0] for t in termsList])
[1767]1209    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
1210    foldWidths = [int(foldNormalizationFactor * term[0]) for term in termsList]
[662]1211    treeStep = 10
1212    treeWidth = {}
1213    for i, term in enumerate(termsList):
[1767]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()]
[662]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
[1767]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
[662]1227##    maxAnnotationTextWidth = width==None and getMaxTextWidthHint([str(t[4]) for t in termsList]+["Annotation"]) or (width - fourthColumnStart - verticalMargin) * 2 / 3
[1767]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
[662]1233    termHeight = font.getsize("A")[1]
1234##    print fourthColumnStart, maxAnnotationTextWidth, verticalMargin
[664]1235    width = sixthColumnStart + maxGenesTextWidth + verticalMargin
[1767]1236    height = len(termsList) * termHeight + 2 * (legendHeight + horizontalMargin)
[662]1237
1238    image = Image.new("RGB", (width, height), backgroundColor)
1239    draw = ImageDraw.Draw(image)
1240
1241    def truncText(text, maxWidth, append=""):
[1767]1242        # print getMaxTextWidthHint([text]), maxAnnotationTextWidth
1243        if getMaxTextWidthHint([text]) > maxWidth:
1244            while getMaxTextWidthHint([text + "..." + append]) > maxWidth and text:
[662]1245                text = text[:-1]
1246            if text:
[1767]1247                text = text + "..." + append
[662]1248            else:
1249                text = append
1250        return text
1251    currentY = horizontalMargin + legendHeight
1252    connectAtX = {}
1253    for i, term in enumerate(termsList):
[1767]1254        draw.line([(verticalMargin, currentY + termHeight / 2), (verticalMargin + foldWidths[i], currentY + termHeight / 2)], width=termHeight - 2, fill=graphColor)
[662]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)
[664]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])
[1767]1260        annotText = width != None and truncText(str(term[5]), maxAnnotationTextWidth)
[664]1261        draw.text((fifthColumnStart, currentY), annotText, font=font, fill=textColor)
[1767]1262        genesText = width != None and truncText(str(term[6]), maxGenesTextWidth) or str(term[6])
[664]1263        draw.text((sixthColumnStart, currentY), genesText, font=font, fill=textColor)
[1767]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)
[662]1268        connectAtX[i] = lineEnd - treeStep
[1767]1269        currentY += termHeight
[662]1270
1271    currentY = horizontalMargin
[664]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)
[662]1278
1279    horizontalMargin = 0
[1767]1280    # draw.line([(verticalMargin, height - horizontalMargin - legendHeight), (verticalMargin + maxFoldWidth, height - horizontalMargin - legendHeight)], width=1, fill=textColor)
[662]1281    draw.line([(verticalMargin, horizontalMargin + legendHeight), (verticalMargin + maxFoldWidth, horizontalMargin + legendHeight)], width=1, fill=textColor)
[1767]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)
[662]1287
[1767]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
[662]1292    image.save(fh)
1293
[1767]1294
[791]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
[1767]1298
1299    maxFoldWidth = width != None and min(150, width / 6) or 150
[791]1300    maxFoldEnrichment = max([t[0] for t in termsList])
[1767]1301    foldNormalizationFactor = float(maxFoldWidth) / maxFoldEnrichment
[791]1302##    foldWidths = [int(foldNormalizationFactor*term[0]) for term in termsList]
1303    foldWidths = [term[0] for term in termsList]
[1767]1304    treeStep = maxFoldEnrichment * 0.05
[791]1305    treeWidth = {}
[745]1306
[791]1307    for i, term in enumerate(termsList):
[1767]1308        treeWidth[i] = (term[-1] == None and treeStep or treeWidth[term[-1]] + treeStep)
[791]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:
[1767]1321            plt.plot([connectAt.get(parent)] * 2, [len(termsList) - i, len(termsList) - parent], color="black")
[791]1322        cellText.append((str(n), str(m), p_val, fdr_val, name, genes))
1323
[1632]1324##    from Orange.orng.orngClustering import TableTextLayout
[791]1325##    text = TableTextLayout((maxFoldEnrichment*1.1, len(termsList)), cellText)
[1632]1326    from Orange.orng.orngClustering import TablePlot
[791]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):
[1767]1333                table.add_cell(i, j, width=len(text), height=1, text=text, loc="left", edgecolor="w", facecolor="w")
[791]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)
[1767]1342##
[791]1343##    table.set_xy(20,20)
1344    plt.show()
[1767]1345
1346
[612]1347class Taxonomy(object):
[1767]1348    """Maps NCBI taxonomy ids to corresponding GO organism codes
[619]1349    """
[1767]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",
[924]1358                      }
[1767]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
[924]1368                "9606": "goa_human",  # Homo sapiens
[1767]1369                "10090": "mgi",  # Mus musculus
[924]1370                "2104": None,  # Mycoplasma pneumoniae
1371                "4530": "gramene_oryza",  # Oryza sativa
1372                "5833": "GeneDB_Pfalciparum",  # Plasmodium falciparum
1373                "4754": None,  # Pneumocystis carinii
[1767]1374                "10116": "rgd",  # Rattus norvegicus
[924]1375                "4932": "sgd",  # Saccharomyces cerevisiae
[1767]1376                "4896": "GeneDB_Spombe",  # Schizosaccharomyces pombe
1377                "31033": None,  # Takifugu rubripes
[924]1378                "8355": None,  # Xenopus laevis
[1767]1379                "4577": None  # Zea mays
[924]1380                }
[822]1381    version = 1
[612]1382    __shared_state = {"tax": None}
[1767]1383
[612]1384    def __init__(self):
1385        self.__dict__ = self.__shared_state
1386        if not self.tax:
[1256]1387            path = orngServerFiles.localpath_download("GO", "taxonomy.pickle")
[612]1388            if os.path.isfile(path):
[633]1389                self.tax = cPickle.load(open(path, "rb"))
[612]1390            else:
1391                orngServerFiles.download("GO", "taxonomy.pickle")
[633]1392                self.tax = cPickle.load(open(path, "rb"))
[1767]1393
[612]1394    def __getitem__(self, key):
[924]1395        key = self.common_org_map.get(key, key)
1396        return self.code_map[key]
[1767]1397
[928]1398    def keys(self):
1399        return list(set(self.common_org_map.keys() + self.code_map.keys()))
[1767]1400
1401
[612]1402def from_taxid(id):
[639]1403    """ Return a set of GO organism codes that correspond to NCBI taxonomy id
1404    """
[924]1405    return Taxonomy()[id]
[612]1406
[1767]1407
[612]1408def to_taxid(db_code):
[639]1409    """ Return a set of NCBI taxonomy ids from db_code GO organism annotations
1410    """
[924]1411    r = [key for key, val in Taxonomy().code_map.items() if db_code == val]
[637]1412    return set(r)
[1767]1413
[612]1414
[538]1415class __progressCallbackWrapper:
1416    def __init__(self, callback):
1417        self.callback = callback
[1767]1418
[538]1419    def __call__(self, bCount, bSize, fSize):
1420        fSize = 10000000 if fSize == -1 else fSize
[1767]1421        self.callback(100 * bCount * bSize / fSize)
1422
[1632]1423from .obiGenomicsUpdate import Update as UpdateBase
[450]1424
[448]1425
1426class Update(UpdateBase):
1427    def __init__(self, local_database_path=None, progressCallback=None):
1428        UpdateBase.__init__(self, local_database_path or getDataDir(), progressCallback)
[1767]1429
[448]1430    def CheckModified(self, addr, date=None):
[674]1431        return date < self.GetLastModified(addr) if date else True
[1767]1432
[448]1433    def CheckModifiedOrg(self, org):
1434        return self.CheckModified("http://www.geneontology.org/gene-associations/gene_association." + org + ".gz", self.LastModifiedOrg(org))
[1767]1435
[448]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)
[450]1441        return datetime.strptime(stream.headers.get("Last-Modified"), "%a, %d %b %Y %H:%M:%S %Z")
1442##        return stream.headers.get("Last-Modified")
[448]1443
[538]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
[448]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])
[1767]1456
[448]1457    def GetDownloadable(self):
[538]1458        orgs = set(self.GetAvailableOrganisms()) - set(self.GetDownloadedOrganisms())
[448]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):
[538]1467        Ontology.DownloadOntology(os.path.join(self.local_database_path, "gene_ontology_edit.obo.tar.gz"), self.progressCallback)
[448]1468        self._update(Update.UpdateOntology, (), self.GetLastModified("http://www.geneontology.org/ontology/gene_ontology.obo"))
[538]1469
[612]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"))
[1767]1473
[612]1474    def UpdateTaxonomy(self, org):
[619]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)
[1767]1486                for taxId in [t.split(":")[-1] for t in taxons if "|" not in t]:  # exclude taxons with cardinality 2
[619]1487                    tax[taxId].add(org)
1488            except Exception, ex:
1489                print ex
[1767]1490
[633]1491        cPickle.dump(dict(tax), open(os.path.join(path, "taxonomy.pickle"), "wb"))
[1767]1492
[612]1493
[538]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)
[1767]1501
1502    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
[538]1503##    profile.runctx("a.GetEnrichedTerms(sorted(a.geneNames)[:100])", {"a":a}, {})
[1767]1504    a.GetEnrichedTerms(sorted(a.geneNames)[:100])  # , progressCallback=_print)
1505    d1 = a.GetEnrichedTerms(sorted(a.geneNames)[:1000])  # , progressCallback=_print)
1506
[538]1507##    print a.GetEnrichedTerms(sorted(a.geneNames)[:100])#, progressCallback=_print)
[664]1508
[1767]1509
[664]1510def _test2():
[726]1511    o = Ontology()
[1256]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))
[1767]1520
[726]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)
[679]1524
[1767]1525
[679]1526def _test3():
1527    o = Ontology()
1528    a = Annotations("sgd", ontology=o)
[726]1529##    a = Annotations(list(a)[3:len(a)/3], ontology=o)
[679]1530    clusterGenes = sorted(a.geneNames)[:1] + sorted(a.geneNames)[-1:]
1531##    clusterGenes = [g + "_" + str(i%5) for g in sorted(a.geneNames)[:2]]
[1767]1532    exonMap = dict([(gene, [gene + "_E%i" % i for i in range(10)]) for gene in a.geneNames])
[679]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))
[1767]1540
[538]1541if __name__ == "__main__":
[726]1542    _test2()
Note: See TracBrowser for help on using the repository browser.