source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1770:48af55b8e628

Revision 1770:48af55b8e628, 60.3 KB checked in by Ales Erjavec <ales.erjavec@…>, 12 months ago (diff)

'Ontology.getitem' now raises a KeyError for an unknown term id.

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