source: orange-bioinformatics/_bioinformatics/obiGeneSets.py @ 1643:2cfa80dac3d3

Revision 1643:2cfa80dac3d3, 16.8 KB checked in by mitar, 2 years ago (diff)

Fixing some imports. Marking widgets as prototypes.

Line 
1"""
2Getting genesets from KEGG and GO.
3
4Maintainer: Marko Toplak
5"""
6
7from __future__ import absolute_import, with_statement
8
9import cPickle as pickle, os, tempfile, sys
10from collections import defaultdict
11
12import orange
13from Orange.orng import orngServerFiles
14
15from . import obiGO, obiKEGG, obiTaxonomy
16
17sfdomain = "gene_sets"
18
19def nth(l,n):
20    return [ a[n] for a in l]
21
22class GeneSet(object):
23
24    def __init__(self, genes=None, name=None, id=None, \
25        description=None, link=None, organism=None, hierarchy=None, pair=None):
26        """
27        pair can be (id, listofgenes) - it is used before anything else.
28        """
29        if genes == None:
30            genes = []
31
32        self.hierarchy = hierarchy
33        self.genes = set(genes)
34        self.name = name
35        self.id = id
36        self.description = description
37        self.link = link
38        self.organism = organism
39
40        if pair:
41            self.id, self.genes = pair[0], set(pair[1])
42
43    """
44    the following functions are needed for sets of gene sets to be able
45    to assess equality
46    """
47
48    def __hash__(self):
49        return self.id.__hash__() + self.name.__hash__()
50
51    def __eq__(self, other):
52        if isinstance(other, self.__class__):
53            return self.__dict__ == other.__dict__
54        else:
55            return False
56
57    def __ne__(self, other):
58        return not self.__eq__(other)
59
60    def size(self):
61        return len(self.genes)
62
63    def cname(self, source=True, name=True):
64        """ Constructs a gene set name with the hierarchy. """
65        oname = self.id
66        if source and self.hierarchy:
67            oname = "[ " + ", ".join(self.hierarchy) + " ] " + oname
68        if name and self.name:
69            oname = oname + " " + self.name
70        return oname
71
72    def to_odict(self, source=True, name=True):
73        """
74        Returns a pair (id, listofgenes), like in old format.
75        """
76        return self.cname(source=source, name=name), self.genes
77
78    def __repr__(self):
79        return "GeneSet(" + ", ".join( [ 
80            "id=" + str(self.id),
81            "genes=" + str(self.genes),
82            "name=" + str(self.name),
83            "link=" + str(self.link),
84            "hierarchy=" + str(self.hierarchy)
85        ]) + ")"
86
87class GeneSetIDException(Exception):
88    pass
89
90class GeneSets(set):
91   
92    def __init__(self, input=None):
93        """
94        odict are genesets in old dict format.
95        gs are genesets in new format
96        """
97        if input != None and len(input) > 0:
98            self.update(input)
99
100    def update(self, input):
101        from . import obiGeneSets
102        if isinstance(input, obiGeneSets.GeneSets):
103            super(GeneSets, self).update(input)
104        elif hasattr(input, "items"):
105            for i, g in input.items():
106                self.add(obiGeneSets.GeneSet(pair=(i, g)))
107        else:
108            for i in input:
109                if isinstance(i, obiGeneSets.GeneSet):
110                    self.add(i)
111                else:
112                    i, g = i
113                    self.add(obiGeneSets.GeneSet(pair=(i, g)))
114
115    def to_odict(self):
116        """ Return gene sets in old dictionary format. """
117        return dict(gs.to_odict() for gs in self)
118
119    def set_hierarchy(self, hierarchy):
120        """ Sets hierarchy for all gene sets """
121        for gs in self:
122            gs.hierarchy = hierarchy
123
124    def __repr__(self):
125        return "GeneSets(" + set.__repr__(self) + ")"
126
127    def common_org(self):
128        """ Returns the common organism. """
129        if len(self) == 0:
130            raise GenesetRegException("Empty gene sets.")
131
132        organisms = set(a.organism for a in self)
133
134        try:
135            return only_option(organisms)
136        except:
137            raise GenesetRegException("multiple organisms: " + str(organisms))
138
139    def hierarchies(self):
140        """ Returns all hierachies """
141        if len(self) == 0:
142            raise GenesetRegException("Empty gene sets.")
143        return set(a.hierarchy for a in self)
144
145    def common_hierarchy(self):
146        hierarchies = self.hierarchies()
147
148        def common_hierarchy1(hierarchies):
149            def hier(l): return set(map(lambda x: x[:currentl], hierarchies))
150            currentl = max(map(len, hierarchies))
151            while len(hier(currentl)) > 1:
152                currentl -= 1
153            return only_option(hier(currentl))
154
155        return common_hierarchy1(hierarchies)
156
157    def split_by_hierarchy(self):
158        """ Splits gene sets by hierarchies. """
159        from . import obiGeneSets
160        hd = dict((h,obiGeneSets.GeneSets()) for h in  self.hierarchies())
161        for gs in self:
162            hd[gs.hierarchy].add(gs)
163        return hd.values()
164
165def goGeneSets(org):
166    """Returns gene sets from GO."""
167    from . import obiGeneSets
168
169    ontology = obiGO.Ontology()
170    annotations = obiGO.Annotations(org, ontology=ontology)
171
172    genesets = []
173    link_fmt = "http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=%s"
174    for termn, term in ontology.terms.items():
175        genes = annotations.GetAllGenes(termn)
176        hier = ("GO", term.namespace)
177        if len(genes) > 0:
178            gs = obiGeneSets.GeneSet(id=termn, name=term.name, genes=genes, hierarchy=hier, organism=org, link=link_fmt % termn) 
179            genesets.append(gs)
180
181    return obiGeneSets.GeneSets(genesets)
182
183def keggGeneSets(org):
184    """
185    Returns gene sets from KEGG pathways.
186    """
187    from . import obiKEGG2 as obiKEGG, obiGeneSets
188   
189    kegg = obiKEGG.KEGGOrganism(org)
190
191    genesets = []
192    for id in kegg.pathways():
193        pway = obiKEGG.KEGGPathway(id)
194        hier = ("KEGG","pathways")
195        gs = obiGeneSets.GeneSet(id=id,
196                                 name=pway.title,
197                                 genes=kegg.get_genes_by_pathway(id),
198                                 hierarchy=hier,
199                                 organism=org,
200                                 link=pway.link)
201        genesets.append(gs)
202
203    return obiGeneSets.GeneSets(genesets)
204
205def omimGeneSets():
206    """
207    Return gene sets from OMIM (Online Mendelian Inheritance in Man) diseses
208    """
209    from . import obiOMIM
210    genesets = [GeneSet(id=disease.id, name=disease.name, genes=obiOMIM.disease_genes(disease), hierarchy=("OMIM",), organism="9606",
211                    link=("http://www.ncbi.nlm.nih.gov/entrez/dispomim.cgi?id=" % disease.id if disease.id else None)) \
212                    for disease in obiOMIM.diseases()]
213    return GeneSets(genesets)
214
215def miRNAGeneSets(org):
216    """
217    Return gene sets from miRNA targets
218    """
219    from . import obimiRNA, obiKEGG
220    org_code = obiKEGG.from_taxid(org)
221    link_fmt = "http://www.mirbase.org/cgi-bin/mirna_entry.pl?acc=%s"
222    mirnas = [(id, obimiRNA.get_info(id)) for id in obimiRNA.ids(org_code)]
223    genesets = [GeneSet(id=mirna.matACC, name=mirna.matID, genes=mirna.targets.split(","), hierarchy=("miRNA", "Targets"),
224                        organism=org, link=link_fmt % mirna.matID) for id, mirna in mirnas]
225    return GeneSets(genesets)
226
227def go_miRNASets(org, ontology=None, enrichment=True, pval=0.05, treshold=0.04):
228    from . import obimiRNA, obiGO
229    mirnas = obimiRNA.ids(int(org))
230    if ontology is None:
231        ontology = obiGO.Ontology()
232         
233    annotations = obiGO.Annotations(org, ontology=ontology)
234   
235    go_sets = obimiRNA.get_GO(mirnas, annotations, enrichment=enrichment, pval=pval, goSwitch=False)
236    print go_sets
237   
238    go_sets = obimiRNA.filter_GO(go_sets, annotations, treshold=treshold)
239   
240    from . import obiGeneSets as gs
241    link_fmt = "http://amigo.geneontology.org/cgi-bin/amigo/term-details.cgi?term=%s"
242    gsets = [gs.GeneSet(id=key, name=ontology[key].name, genes=value, hierarchy=("miRNA", "go_sets",),
243                        organism=org, link=link_fmt % key) for key, value in go_sets.items()]
244    gset = gs.GeneSets(gsets)
245    return gset
246
247
248def loadGMT(contents, name):
249    """
250    Eech line consists of tab separated elements. First is
251    the geneset name, next is it's description.
252   
253    For now the description is skipped.
254    """
255
256    from . import obiGeneSets
257
258    def hline(s):
259        tabs = [tab.strip() for tab in s.split("\t")]
260        return obiGeneSets.GeneSet(id=tabs[0], description=tabs[1],
261                                   hierarchy=(name,), genes=tabs[2:])
262
263    def handleNELines(s, fn):
264        """
265        Run function on nonempty lines of a string.
266        Return a list of results for each line.
267        """
268        lines = (l.strip() for l in s.splitlines())
269        return [fn(l) for l in lines if l]
270
271    return obiGeneSets.GeneSets(handleNELines(contents, hline))
272
273"""
274We have multiple paths for gene set data:
275buffer/bigfiles/gene_sets
276and
277buffer/gene_sets_local
278both have available.txt
279"""
280
281def omakedirs(dir):
282    try:
283        os.makedirs(dir)
284    except OSError:
285        pass
286
287def local_path():
288    """ Returns local path for gene sets. Creates it if it does not exists
289    yet. """
290    from Orange.orng import orngEnviron
291    pth = os.path.join(orngEnviron.directoryNames["bufferDir"], "gene_sets_local")
292    omakedirs(pth)
293    return pth
294
295def build_index(dir):
296    """ Returns gene set availability index for some folder. """
297    pass
298
299class GenesetRegException(Exception): pass
300
301def only_option(a):
302    if len(a) == 1:
303        return list(a)[0]
304    else:
305        raise Exception()
306
307def filename(hierarchy, organism):
308    """ Obtain a filename for given hierarchy and organism. """
309    return "gs_" + "_._".join(hierarchy + \
310        (organism if organism != None else "",)) + ".pck"
311
312def filename_parse(fn):
313    """ Returns a hierarchy and the organism from the filename."""
314    fn = fn[3:-4]
315    parts = fn.split("_._")
316    hierarchy = tuple(parts[:-1])
317    org = parts[-1] if parts[-1] != "" else None
318    return hierarchy, org
319
320def is_genesets_file(fn):
321    return fn.startswith("gs_") and fn.endswith(".pck")
322
323def list_local():
324    """ Returns available gene sets from the local repository:
325    a list of (hierarchy, organism, on_local) """
326    pth = local_path()
327    gs_files = filter(is_genesets_file, os.listdir(pth))
328    return [ filename_parse(fn) + (True,) for fn in gs_files ]
329   
330def list_serverfiles_from_flist(flist):
331    gs_files = filter(is_genesets_file, flist)
332    localfiles = set(orngServerFiles.listfiles(sfdomain))
333    return [ filename_parse(fn) + \
334        ((True,) if fn in localfiles else (False,)) for fn in gs_files ]
335
336def list_serverfiles_conn(serverfiles=None):
337    """ Returns available gene sets from the server files
338    repository: a list of (hierarchy, organism, on_local) """
339    if serverfiles == None:
340        serverfiles = orngServerFiles.ServerFiles()
341    flist = serverfiles.listfiles(sfdomain)
342    return list_serverfiles_from_flist(flist)
343
344def list_serverfiles():
345    fname = orngServerFiles.localpath_download(sfdomain, "index.pck")
346    flist = pickle.load(open(fname, 'r'))
347    return list_serverfiles_from_flist(flist)
348
349def list_all():
350    """
351    return a list of (hier, org, avalable_locally)
352    If something for a specific (hier, org) is not downloaded
353    yet, show it as not-local. """
354    flist = list_local() + list_serverfiles()
355    d = {}
356    for h,o,local in flist:
357        d[h,o] = min(local, d.get((h,o),True))
358    return [ (h,o,local) for (h,o),local in d.items() ]
359
360def update_server_list(serverfiles_upload, serverfiles_list=None):
361    if serverfiles_list == None:
362        serverfiles_list = orngServerFiles.ServerFiles()
363
364    flist = map(lambda x: filename(*x[:2]), list_serverfiles_conn(serverfiles_list))
365
366    tfname = pickle_temp(flist)
367   
368    try:
369        fn = "index.pck"
370        title = "Gene sets: index"
371        tags = [ "gene sets", "index", "essential" ]
372        serverfiles_upload.upload(sfdomain, fn, tfname, title, tags)
373        serverfiles_upload.unprotect(sfdomain, fn)
374    except Exception,e:
375        raise e
376    finally:
377        os.remove(tfname)
378
379def register_local(genesets):
380    """ Registers using the common hierarchy and organism. """
381    pth = local_path()
382
383    org = genesets.common_org()
384    hierarchy = genesets.common_hierarchy()
385    fn = filename(hierarchy, org)
386
387    with open(os.path.join(pth, fn), "w") as f:
388        pickle.dump(genesets, f)
389
390    return fn
391
392def pickle_temp(obj):
393    """ Pickle a file to a temporary file returns its name """
394    fd,tfname = tempfile.mkstemp()
395    os.close(fd)
396    f = open(tfname, 'wb')
397    pickle.dump(obj, f)
398    f.close()
399    return tfname
400
401def register_serverfiles(genesets, serverFiles):
402    """ Registers using the common hierarchy and organism. """
403    org = genesets.common_org()
404    hierarchy = genesets.common_hierarchy()
405    fn = filename(hierarchy, org)
406
407    #save to temporary file
408    tfname = pickle_temp(genesets)
409   
410    try:
411        taxname = obiTaxonomy.name(org)
412        title = "Gene sets: " + ", ".join(hierarchy) + \
413            ((" (" + taxname + ")") if org != None else "")
414        tags = list(hierarchy) + [ "gene sets", taxname ] + \
415            ([ "essential" ] if org in obiTaxonomy.essential_taxids() else [] )
416        serverFiles.upload(sfdomain, fn, tfname, title, tags)
417        serverFiles.unprotect(sfdomain, fn)
418    finally:
419        os.remove(tfname)
420
421    update_server_list(serverFiles)
422
423def register(genesets, serverFiles=None):
424    """
425    Hierarchy is induced from the gene set names.
426    """
427    if serverFiles == None:
428        register_local(genesets)
429    else:
430        register_serverfiles(genesets, serverFiles)
431
432def build_hierarchy_dict(files):
433    hierd = defaultdict(list)
434    for ind,f in enumerate(files):
435        hier, org = f
436        for i in range(len(hier)+1):
437            hierd[(hier[:i], org)].append(ind)
438    return hierd
439
440def load_local(hierarchy, organism):
441    files = map(lambda x: x[:2], list_local())
442
443    hierd = build_hierarchy_dict(files)
444
445    out = GeneSets()
446    for (h, o) in [ files[i] for i in hierd[(hierarchy, organism)]]:
447        fname = os.path.join(local_path(), filename(h, o))
448        out.update(pickle.load(open(fname, 'r')))
449    return out
450
451def load_serverfiles(hierarchy, organism):
452    files = map(lambda x: x[:2], list_serverfiles())
453    hierd = build_hierarchy_dict(files)
454
455    out = GeneSets()
456    for (h, o) in [ files[i] for i in hierd[(hierarchy, organism)]]:
457        fname = orngServerFiles.localpath_download(sfdomain, 
458            filename(h, o))
459        out.update(pickle.load(open(fname, 'r')))
460    return out
461
462def load(hierarchy, organism):
463    """ First try to load from the local registred folder, then
464    from the server files. """
465    ret = load_local(hierarchy, organism)
466    ret.update(load_serverfiles(hierarchy, organism))
467    return ret
468
469def collections(*args):
470    """
471    Input is a list of collections.
472    Collection can either be a tuple (hierarchy, orgranism), where
473    hierarchy is a tuple also.
474    """
475    from . import obiGeneSets
476    result = obiGeneSets.GeneSets()
477
478    for collection in args:
479        try:
480            result.update(collection)
481        except ValueError:
482            if issequencens(collection): #have a hierarchy, organism specification
483                new = load(*collection)
484                result.update(new)
485            else:
486                if collection.lower()[-4:] == ".gmt": #format from webpage
487                    result.update(loadGMT(open(collection,"rt").read(), collection))
488                else:
489                    raise Exception("collection() accepts files in .gmt format only.")
490
491    return result
492
493def issequencens(x):
494    "Is x a sequence and not string ? We say it is if it has a __getitem__ method and it is not an instance of basestring."
495    return hasattr(x, '__getitem__') and not isinstance(x, basestring)
496
497class TException(Exception): pass
498
499def upload_genesets(rsf):
500    """
501    Builds the default gene sets and
502    """
503    orngServerFiles.update_local_files()
504
505    genesetsfn = [ keggGeneSets, goGeneSets, miRNAGeneSets]
506    organisms = obiTaxonomy.common_taxids()
507    for fn in genesetsfn:
508        for org in organisms:
509        #for org in [ "9606" ]:
510            print "Uploading ORG", org, fn
511            try:
512                genesets = fn(org).split_by_hierarchy()
513                for gs in genesets:
514                    print "registering", gs.common_hierarchy()
515                    register_serverfiles(gs, rsf)
516                    print "successful", gs.common_hierarchy()
517            except Exception, e:
518                print "Not successful"
519
520if __name__ == "__main__":
521
522    #gs = keggGeneSets("9606")
523    akegg = collections((("KEGG",),"9606"))
524    a = collections({"PATH": ["g1","g2","g3"]}, "steroltalk.gmt", akegg)
525    for e in a:
526        print e
527    #print len(collections(keggGeneSets("9606"),(("KEGG",),"9606"), "C5.BP.gmt"))
528    #print len(collections((("KEGG",),"9606"), "C5.BP.gmt"))
529    #print sorted(list_all())
530    #print len(collections((("KEGG",),"9606"), (("GO",), "9606"), "C5.BP.gmt"))
531    #register_local(keggGeneSets("9606"))
532    #register_local(goGeneSets("9606"))
533    #register_serverfiles(gs, rsf)
534    #print list_serverfiles_conn()
535    #print "Server list from index", list_serverfiles()
536    rsf = orngServerFiles.ServerFiles(username=sys.argv[1], password=sys.argv[2])
537    upload_genesets(rsf)
Note: See TracBrowser for help on using the repository browser.