source: orange-bioinformatics/_bioinformatics/obiGO.py @ 1716:4c4266f7c5a5

Revision 1716:4c4266f7c5a5, 62.3 KB checked in by markotoplak, 20 months ago (diff)

Removed the old version of obiKEGG. Renamed obiKEGG2 -> obiKEGG.

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