source: orange-bioinformatics/_bioinformatics/obiGeneSets.py @ 1707:7d692228fea1

Revision 1707:7d692228fea1, 16.3 KB checked in by markotoplak, 20 months ago (diff)

Speed up of setsig (obiGeneSetSig).

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