source: orange-bioinformatics/_bioinformatics/obiDicty.py @ 1828:857938f13079

Revision 1828:857938f13079, 57.7 KB checked in by Flashpoint <vid.flashpoint@…>, 10 months ago (diff)

Added an organisms shortname() function to the obiTaxonomy.py to add short names to geneset tags

RevLine 
[420]1import sys, pprint, time, re
2from itertools import *
[1029]3import urllib
[557]4import urllib2
[420]5import socket
[792]6import os
[889]7from collections import defaultdict
[1021]8import pickle
[420]9
[1632]10import orange
11from Orange.orng import orngServerFiles
12
[1148]13defaddress = "http://bcm.fri.uni-lj.si/microarray/api/index.php?"
[1292]14#defaddresspipa = "https://pipa.fri.uni-lj.si/pipa/script/api/orange.py?action="
15defaddresspipa = "https://pipa.fri.uni-lj.si/PIPA/PIPAapi/PIPAorange.py"
[1510]16defaddresspipa_coverage = "https://pipa.biolab.si/tbrowse/api/pipa.py"
[1292]17
18pipaparuser = "pipa_username"
19pipaparpass = "pipa_password"
[1148]20
[420]21#utility functions - from Marko's mMisc.py
22
23def splitN(origsize, maxchunk):
24    """
25    Splits an integer into chunks of given size. Each created chunk
26    except possibly the last one is of maximum allowed size.
27    Chunks are returned in list.
28    """
29    l = [maxchunk]*(origsize/maxchunk)
30    a = origsize % maxchunk
31    if a > 0: 
32        l.append(a)
33    return l
34
35def split(l, maxchunk):
36    """
37    Splits list l into chunks of size maxchunk. Each created chunk
38    except possibly the last one is of maximum allowed size.
39    """
40    sizes = splitN(len(l), maxchunk)
41    out = []
42    tillNow = 0
43    for s in sizes:
44        out.append(l[tillNow:tillNow+s])
45        tillNow += s
46    return out       
47
48def lloc(l,n):
49    """
50    List location in list of list structure.
51    Enable the use of negative locations:
52    -1 is the last element, -2 second last...
53    """
54    if n < 0:
55        return len(l[0])+n
56    else:
57        return n
58
59def loc(l,n):
60    """
61    List location.
62    Enable the use of negative locations:
63    -1 is the last element, -2 second last...
64    """
65    if n < 0:
66        return len(l)+n
67    else:
68        return n
69
70def nth(l,n):
71    """
72    Returns only nth elemnt in a list.
73    """
74    n = lloc(l,n)
75    return [ a[n] for a in l ]
76
77def imnth(l, ns):
78    """
79    Return only columns as specified in ns. Returns an generator.
80    """
81    ns = [ lloc(l,n) for n in ns ]
82    for a in l:
83        yield [ a[n] for n in ns ]
84
85def flatten(l,r=0):
86    """
87    Flatten a python structure into a list. Leave strings alone.
88    """
89    if type(l) == type("a"):
90        return [ l ]
91    try: #if enumerable then flatten it's elements
92        rec = [ flatten(a,r=r+1) for a in l ]
93        ret = []
94        for a in rec:
95            ret = ret + a
96        return ret
97    except:
98        return [ l ]
99
100def mxrange(lr):
101    """
102    Multiple xranges. Can be used to traverse matrices.
103    This function is very slow due to unknown number of
104    parameters.
105
106    >>> mxrange([3,5])
107    [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]
108    >>> mxrange([[3,5,1],[9,0,-3]])
109    [(3, 9), (3, 6), (3, 3), (4, 9), (4, 6), (4, 3)]
110    """
111    if len(lr) == 0:
112        yield ()
113    else:
114        #it can work with single numbers
115        index = lr[0]
116        if type(1) == type(index):
117            index = [ index ]
118        for a in range(*index):
119            for b in mxrange(lr[1:]):
120                yield tuple([a] + list(b))
121
[527]122
123def issequencens(x):
124    """
125    Is x a sequence and not string ? We say it is if it has a __getitem__
126    method and it is not an instance of basestring.
127    """
128    return hasattr(x, '__getitem__') and not isinstance(x, basestring)
129
[420]130#end utility functions
131
[1371]132import statc
133#median = statc.median
134
135def median(l):
136    if len(l) == 0:
137        return None
138    else:
139        return statc.median(l)
140
141
142
[1172]143socket.setdefaulttimeout(60)
[420]144
145verbose = 0
146
147class HttpGetException(Exception): pass
148
[527]149def replaceChars(address):
150    return address.replace(" ", "%20")
151
[420]152def httpGet(address, *args, **kwargs):
153    if verbose: 
[1172]154        print address, args, kwargs,  " "
[527]155    address = replaceChars(address)
[420]156    t1 = time.time()
[557]157    f = urllib2.urlopen(address, *args, **kwargs)
[420]158    read = f.read()
159    if verbose:
[1021]160        print "bytes", len(read),
[527]161    if verbose:
[420]162        print time.time() - t1
163    return read
164
165def txt2ll(s, separ=' ', lineSepar='\n'):
166    return [ a.split(separ) for a in s.split(lineSepar) ]
167
[1059]168class AuthenticationError(Exception):
169    pass
170
[420]171class DBInterface(object):
172 
[1029]173    def __init__(self, address):
[420]174        self.address = address
175
[1029]176    def raw(self, request, data=None, tryN=3):
[420]177        if verbose:
178            print "tryN", tryN
179
[1376]180        def try_down():
181            try:
182                if data == None:
183                    return httpGet(self.address + request)
184                else:
185                    return httpGet(self.address + request, data=urllib.urlencode(data))
186            except urllib2.HTTPError, e:
187                if e.code == 403:
188                    raise AuthenticationError()
189                else:
190                    raise e
191
192        if tryN > 1:
193            try:
[1377]194                return try_down()
[1376]195            except IOError:
196                return self.raw(request, data=data, tryN=tryN-1)
197        else:
[1377]198            return try_down()
[420]199
[1029]200    def get(self, request, data=None, tryN=3):
201        rawf = self.raw(request, data)
[420]202        if rawf == None:
[1377]203            print "rawf == None!"
[420]204            raise Exception("Connection error when contacting " + self.address + request)
[1059]205        if rawf.startswith("error: authentication failed"):
206            raise AuthenticationError()
207        elif rawf[:1] == "<" or rawf[:5] == "error" or rawf.startswith("MOD_PYTHON ERROR"): #an error occurred - starting some html input
[831]208            #TODO are there any other kinds of errors?
209            if tryN > 0:
[1017]210                if verbose:
211                    print "trying again"
[1029]212                return self.get(request, data=data, tryN=tryN-1)
[831]213            else:
[1214]214                if verbose:
215                    print rafw[:1000]
[831]216                raise Exception("Error with the database")
217
[420]218        a = txt2ll(rawf, separ='\t')
[831]219       
[420]220        if a[-1][0] == "": #remove empty line on end
221            a = a[:-1]
222        return a
223
224def _test():
225    import doctest
226    doctest.testmod()
227
228def splitTableOnColumn(ll,n):
229    omap = {}
230    n = lloc(ll,n)
231    for l in ll:
232        cell = omap.get(l[n], [])
233        cell.append(l)
234        omap[l[n]] = cell
235    return omap
236
237def neededColumns(legend, want):
238    return [ legend.index(a) for a in want ]
239
240def onlyColumns(ll, legend, want):
241    return list(imnth(ll, neededColumns(legend, want)))
242
243def ll2dic(ll, key=0, value=1):
244    """
245    Converts LL to map. Key is key position, value value position
246    """
247    names = nth(ll, lloc(ll, key))
248    if len(names) == len(set(names)):
249        return dict(list(imnth(ll, [ lloc(ll, key),  lloc(ll, value) ] )))
250    else:
251        raise Exception("all keys are not unique")
252
253def dic2ll(dic):
254
255    columns = sorted(dic.values()[0].keys())
256    ll = []
257    for key, d in dic.items():
258        ll.append([ key ] + [ d[a] for a in columns ])
259    return columns,ll
260
261def allUnique(els):
262    return len(els) == len(set(els))
263
264def reverseDic(d):
265    """
266    Create a reverse dictionary if a unique reverse is possible.
267    """
268    if allUnique(d.values()):
269        return dict([ (b,a) for a,b in d.items() ])
270
271def chainLookup(a, dics, force=[]):
272    """
273    Goes through a list of dictionaries. On each step try to
274    transform current key to value. The values becomes the key
275    for next search If unsuccessfull, query the next
276    dictionary with current key. Force last translation.
277    """
278    force = [ loc(dics, i) for i in force ]
279    for i,dic in enumerate(dics):
280        if i in force: #force specified translations
281            a = dic[a]
282        else:
283            if a in dic:
284                a = dic[a]
285    return a
286
[1021]287class DBCommon(object):
[420]288
[1021]289    def fromBuffer(self, addr):
290        return self.buffer.get(self.address + addr)
291
[1178]292    def toBuffer(self, addr, cont, version, autocommit=True):
[1021]293        if self.buffer:
[1178]294            return self.buffer.add(self.address + addr, cont, version=version, autocommit=autocommit)
295
296    def bufferCommit(self):
297        if self.buffer:
298            self.buffer.commit()
[1021]299
[1049]300    def bufferFun(self, bufkey, bufver, reload, fn, *args, **kwargs):
[1021]301        """
302        If bufkey is already present in buffer, return its contents.
303        If not, run function with arguments and save its result
304        into the buffer.
305        """
[1049]306        if self.inBuffer(bufkey) == bufver and reload == False:
[1021]307            res = self.fromBuffer(bufkey)
308        else:
309            res = fn(*args, **kwargs)
[1049]310            self.toBuffer(bufkey, res, bufver)
[1021]311        return res
312
[1510]313    def sq(self, s1, data=None, buffer=True, bufadd="", bufname=None, bufver="0", reload=False, bufferkey=None, db=None):
314        if db == None:
315            db = self.db
[1021]316        if buffer:
[1049]317            bufkey = bufadd + (bufname if bufname != None else s1)
[1155]318            if bufferkey != None:
[1148]319                bufkey = bufferkey(bufkey, data)
[1510]320            res = self.bufferFun(bufkey, bufver, reload, db.get, s1, data=data)
[1021]321        else:
[1510]322            res = db.get(s1, data=data)
[1021]323        return res[1:],res[0]
324
325    def inBuffer(self, addr):
326        if self.buffer:
327            return self.buffer.contains(self.address + addr)
328        else:
329            return False
330
331    def dictionarize(self, ids, fn, *args, **kwargs):
332        """
333        Creates a dictionary from id: function result.
334        Callback for each part done.
335        """
336        callback = kwargs.pop("callback", None)
337        odic = {}
338        for a,b in izip(ids, fn(*args, **kwargs)):
339            odic[a] = b
340            if callback: callback()
341        return odic
342
[1049]343    def downloadMulti_bufcommand_replace_multi(self, command, data=None, chunk=100, bufferkey=None, transformfn=None):
[1021]344        """
[1049]345        Get function which gives buffer address for an id and a function
346        which replaces $MULTI$.
[1021]347        """
348
[1029]349        def bufferkey1(command, data):
[1021]350            if transformfn:
351                return "TRANS " + command
352            else:
353                return command
354
[1029]355        def replace_multi(command, data, repl):
356            return command.replace("$MULTI$", repl),\
[1234]357                dict((a,b.replace("$MULTI$", repl)) for a,b in sorted(data.items())) if data != None else None
[1029]358
359        if bufferkey == None:
360            bufferkey=bufferkey1
361
362        bufcommand = lambda x, c=command, d=data: bufferkey(*replace_multi(c, d, x))
[1049]363        return bufcommand, replace_multi
364
365    def downloadMulti(self, command, ids, data=None, chunk=100, transformfn=None, bufferkey=None, separatefn=None, bufreload=False, bufver="0"):
366        """
367        Downloads multiple results at once.
368        Results in the same order as in ids.
369
370        Bufferkey transforms command and data into buffer key.
371        bufver is a function returning buffer version for a given id. if
372            a string is given, use it for all ids
373        """
374
375        sids = split(ids,chunk)
376   
377        bufverfn = None
378        if isinstance(bufver, basestring):
379            bufverfn = lambda x: bufver
380        else:
381            bufverfn = bufver
382
383        bufcommand, replace_multi = self.downloadMulti_bufcommand_replace_multi(command, data=data, chunk=chunk, bufferkey=bufferkey, transformfn=transformfn)
[1029]384
[1021]385        for i,sidp in enumerate(sids):
386
387            buffered = []
388            unbuffered = []
389       
390            for a in sidp:
[1049]391                if self.inBuffer(bufcommand(a)) == bufverfn(a) and bufreload == False:
[1021]392                    buffered.append(a)
393                else:
394                    unbuffered.append(a)
395
396            res = []
397            legend = []
398
399            if len(unbuffered) > 0:
[1029]400                com1, d1 = replace_multi(command, data, ",".join(unbuffered))
[1049]401                res, legend = self.sq(com1, data=d1, buffer=False) #get unbuffered part
[1021]402            else:
403                # get legend from buffer also
[1029]404                legend = self.fromBuffer(bufcommand(buffered[0]))[0]
[1021]405
406            #split on different values of the first column - first attribute
407
408            if not separatefn:
409                antss = splitTableOnColumn(res, 0)
410            else:
411                legend, antss = separatefn([legend]+res)
412
413            #if transform before saving is requested, do it
414            if transformfn:
415                nantss = {}
416                nlegend = None
417                for a,b in antss.items():
418                    nb, nlegend = transformfn(b, legend)
419                    nantss[a] = nb
420                legend = nlegend
421                antss = nantss
422 
423            #here save buffer
424            for a,b in antss.items():
[1178]425                self.toBuffer(bufcommand(a), [ legend ] + b, bufverfn(a), autocommit=False)
426            self.bufferCommit()
427           
[1021]428
[1049]429            #get buffered from the buffer
[1029]430            antssb = dict([ (b, self.fromBuffer(bufcommand(b))[1:]) for b in buffered ])
[1021]431            antss.update(antssb)
432
433            #put results in order
434            tl = []
435            for ci in sidp:
436                yield antss[ci], legend
437
[1154]438    def exampleTables(self, ids, chipsm=None, spotmap={}, callback=None, exclude_constant_labels=False, annots={}, chipfn=None, allowed_labels=None):
[1021]439        """
440        Create example tables from chip readings, spot mappings and
441        group specifications.
442
443        group is the output from "sortAnnotations" function.
444        spotmap is a dictionary of { spotid: gene }
445        chipsm is a dictionary of chip readings
446
447        Callback: number of chipids + 2
448        """
449
450        if verbose:
451            print "Creating example table"
452
453        if callback: callback()
454
455        amap = {}
456        amapnext = 0
457
458        togen = []
459
460        groupnames = []
461        groupvals = []
462        groupannots = []
463
464        if chipsm == None:
465            chipdl = chipfn(ids)
466
467        for chipid in ids:
468
469            if chipsm != None:
470                chipdata = chipsm[chipid]
471            else:
472                chipdata = chipdl.next()
473
474            if callback: callback()
475
476            #add to current position mapping
477            repeats = {}
478            for id,_ in chipdata:
479                rep = repeats.get(id, 0)
480                repeats[id] = rep+1
481                key = (id, rep)
482                if key not in amap:
483                    amap[key] = amapnext
484                    amapnext += 1
485
486            vals = [ None ] * len(amap)
487
488            repeats = {}
489            for id,v in chipdata:
490                rep = repeats.get(id, 0)
491                repeats[id] = rep+1
492                key = (id, rep)
493                putind = amap[key]
494                vals[putind] = v
495            groupvals.append(vals)
496
497            groupnames.append(chipid) 
498
499            newannots = [['id', str(chipid)]] #add chipid to annotations
500            if annots:
501                newannots += annots[chipid]
502            groupannots.append(newannots)
503
504        togen = (groupnames, groupvals, groupannots)
505       
506        if callback: callback()
507
508        ddb = [ None ]*len(amap)
509        for (a,rep),pos in amap.items():
510            if len(spotmap):
511                ddb[pos] = spotmap.get(a, "#"+a)
512            else:
513                ddb[pos] = a
514
515        #this is sorted position mapping: key -> sortedind
516        posMap = dict( (k,i) for i,k in enumerate(sorted(amap.keys())) )
517        revmap = dict( ( (i,k) for k,i in amap.items() ) )
518        #permutation[i] holds target of current [i]
519        permutation = [ posMap[revmap[i]] for i in range(len(amap)) ]
520
521        def enlength(a, tlen):
522            """ Adds Nones to the end of the list """
523            if len(a) < tlen:
524                return a + [ "None" ]*(tlen-len(a))
525            else:
526                return a
527
528        def enlengthl(l, tlen):
529            return [ enlength(a, tlen) for a in l ]
530   
531        groupnames, groupvals, groupannots = togen
532
533        et = createExampleTable(groupnames, 
534            enlengthl(groupvals, len(ddb)),
[1154]535            groupannots, ddb, exclude_constant_labels=exclude_constant_labels, permutation=permutation, allowed_labels=allowed_labels)
[1021]536
537        if callback: callback()
538
539        return et
540
[1029]541def bufferkeypipa(command, data):
542    """ Do not save password to the buffer! """
[1212]543    command = command + " v5" #add version
[1148]544    if data != None:
545        data = data.copy()
[1292]546        if pipaparpass in data:
547            data.pop(pipaparpass)
[1148]548        return command + " " +  urllib.urlencode(sorted(data.items()))
549    else:
550        return command
[1021]551
552class PIPA(DBCommon):
553
554    def __init__(self, address=defaddresspipa, buffer=None, username=None, password=None):
555        self.address = address
556        self.db=DBInterface(address)
[1510]557        self.db_coverage=DBInterface(defaddresspipa_coverage) #temporary
[1021]558        self.buffer = buffer
[1029]559        self.username = None
560        if username != None:
561            self.username = username
562            self.password = password
563
564    def add_auth(self, data=None):
565        if self.username == None:
566            return data
[1292]567        authdic = { pipaparuser: self.username, pipaparpass: self.password }
[1029]568        if data != None:
569            authdic.update(data)
570        return authdic
[1021]571
[1234]572    def annotations(self, reload=False, bufver="0"):
573        """
574        Returns a dictionary of (id: dictionary of annotations).
575        """
[1292]576        #res, legend = self.sq("mapping_list", data=self.add_auth(), reload=reload, bufferkey=bufferkeypipa, bufver=bufver)
577        res, legend = self.sq("", data=self.add_auth({"action": "expression_list"}), reload=reload, bufferkey=bufferkeypipa, bufver=bufver)
[1234]578        return dict( (sa[0], dict(zip(legend[1:], sa[1:]))) for sa in res )
[1021]579
[1234]580    def chips(self, ids, ctype, reload=False, bufver="0"):
[1021]581        def separatefn(res):
582            #each one is own rown
583            #genes are in the first row
[1061]584            #remove unknown
[1021]585            genes = res[0][1:]
586            cids = nth(res,0)[1:]
587
588            antss = {}
589            for i,cid in enumerate(cids):
590                row = i+1
591                vals = res[row][1:]
[1061]592                antss[cid] = [ list(a) for a in zip(genes, vals) if a[1] != "?" ]
[1021]593            return ['gene_id', 'value'], antss
594
[1234]595        download_command = "gene_expression"
[1292]596        datadict = {"action": download_command, "ids":"$MULTI$", 'transpose':'1'}
597        if ctype != None:
598            datadict["type"] = ctype
599        antss = self.downloadMulti("", ids, data=self.add_auth(datadict), chunk=10, separatefn=separatefn, bufferkey=bufferkeypipa, bufreload=reload, bufver=bufver)
[1021]600        for a,legend in antss:
601            yield a
602
[1565]603    def coverage(self, id_, genome, chromosome, reload=False, bufver="0", raw=False):
[1510]604        data = self.add_auth({"action": "coverage_strand", "chr_name": str(chromosome), 
[1565]605            "pixel_size": "1", "genome": genome, "sample_id": str(id_), "raw": str(raw).lower()})
[1510]606
607        def comp():
608
609            d = dict(data)
610
611            def get():
612                #sq returns [ "all lines after first", "firstline"] - take first column!
613                #do not use buffering on this stage!
614                a = self.sq("", data=d, reload=reload, buffer=False, bufferkey=bufferkeypipa, bufver=bufver, db=self.db_coverage)[1][0].strip("[]")
615                return map(int, a.split(", "))
616               
617            d["strand"] = "plus"
618            plus = get()
619            d["strand"] = "minus"
620            minus = get()
621
622            import numpy
623            return numpy.array([plus, minus])
624       
625        bufkey =  bufferkeypipa("", data)
626        return self.bufferFun(bufkey, bufver, reload, comp)
627
628
[1234]629    def gene_expression_types(self, reload=False, bufver="0"):
[1292]630        #res, legend = self.sq("gene_expression_type", data={self.add_auth(), reload=reload, bufferkey=bufferkeypipa, bufver=bufver)
631        res, legend = self.sq("", data=self.add_auth({"action":"gene_expression_type"}), reload=reload, bufferkey=bufferkeypipa, bufver=bufver)
[1234]632        return sorted(tuple(a) for a in res)
[1049]633
[1234]634    def chips_keynaming(self, ctype):
[1292]635        keynamingfn,_ = self.downloadMulti_bufcommand_replace_multi("", data=self.add_auth({"action": "gene_expression", "ids":"$MULTI$", 'transpose':'1', 'type':ctype}), chunk=100, bufferkey=bufferkeypipa, transformfn=None)
[1210]636        return keynamingfn
637
[1021]638    def get_data(self, exclude_constant_labels=False, average=median, 
[1292]639        ids=None, callback=None, bufver="0", transform=None, ctype=None, allowed_labels=None):
[1021]640        """
641        Get data in a single example table with labels of individual attributes
642        set to annotations for query and post-processing
643        instructions.
644
645        Parameters:
646            average: function used for combining multiple reading of the same spot on
647                a chip. If None, no averaging is done. Fuction should take a list
648                of floats and return an "averaged" float.
649            ids: a list of chip ids. If absent, make a search
650            exclude_constant_labels: if a label has the same value in whole
651                example table, remove it
652            format: if short, use short format for chip download
[1234]653            ctype: expression type, from gene_expression_types
[1021]654
655        Defaults: Median averaging.
656        """
657
658        def optcb():
659            if callback: callback()
660
661        cbc = CallBack(1, optcb, callbacks=10)
662
663        if not ids:
[1371]664            return None
[1021]665
666        cbc.end()
667
668        #downloads annotations
669        cbc = CallBack(len(ids), optcb, callbacks=10)
670
[1234]671        readall = self.annotations()
[1021]672
673        read = {}
674        for a,b in readall.items():
[1234]675            read[a] = b.items()
[1021]676
677        cbc.end()
678
679        #till now downloads were small
680
681        import time
682        tstart = time.time()
683
684        #here download actually happens
685        chipfn = None
686
[1234]687        chipfn = lambda x: self.chips(x, ctype, bufver=bufver)
[1021]688       
689        if verbose:
690            print "DOWNLOAD TIME", time.time() - tstart
691
692        cbc = CallBack(len(ids)*2+len(ids)+1, optcb, callbacks=999-30)
[1154]693        et = self.exampleTables(ids, spotmap={}, callback=cbc, annots=read, exclude_constant_labels=exclude_constant_labels, chipfn=chipfn, allowed_labels=allowed_labels)
[1021]694        cbc.end()
695
696        cbc = CallBack(1, optcb, callbacks=10)
697
[1148]698        #transformation is performed prior to averaging
699        if transform != None:
700            transformValues(et, fn=transform) #in place transform
701            cbc()
702
[1021]703        #if average function is given, use it to join same spotids
704        if average != None:
705            et = averageAttributes(et, fn=average)
706            cbc()
707
708        cbc.end()
709
710        return et
711
[1619]712def bufferkeypipax(command, data):
713    """ Do not save password to the buffer! """
[1686]714    command = command + " v8" #add version
[1619]715    if data != None:
716        data = data.copy()
717        if pipaparpass in data:
718            data.pop(pipaparpass)
719        return command + " " +  urllib.urlencode(sorted(data.items()))
720    else:
721        return command
[1615]722
[1620]723
[1615]724class PIPAx(PIPA):
725    """`PIPAx <http://pipa.biolab.si/?page_id=23>` api interface.
726    """
727
728    API_ADDRESS = "https://pipa.biolab.si/pipax/api.py"
729
730    def __init__(self, address=API_ADDRESS, buffer=None,
731                 username=None, password=None):
732        self.address = address
733        self.db = DBInterface(address)
734        self.buffer = buffer
735        self.username = username
736        self.password = password
737
738    def genomes(self, reload=False, bufver="0"):
739        """Return a list of available genomes as a list of
740        (genome_id, genome_name) tuples.
741
742        >>> pipax.genomes()
743        [('dd', 'Dictyostelium discoideum'), ('dp', 'Dictyostelium purpureum')]
744
745        """
746        data = {"action": "genomes"}
747        res, _ = self.sq("", data=data, reload=reload, bufver=bufver)
748        return [tuple(r) for r in res]
749
750    def mappings(self, reload=False, bufver="0"):
751        """Return a dictionary of {unique_id: dictionary_of_annotations}
752        where the keys for dictionary_of_annotations are
753        ["id", data_id", "data_name", "genomes_id"]
754
755        >>> mappings = pipax.mappings()
756
757        """
758        data = self.add_auth({"action": "mappings"})
759        res, legend = self.sq("", data=data, reload=reload,
[1619]760                              bufferkey=bufferkeypipax,
[1615]761                              bufver=bufver)
762
763        return dict((sa[0], dict(zip(legend[1:], sa[1:]))) for sa in res)
764
765    def result_types(self, reload=False, bufver="0"):
766        """Return a list of available result type templates.
767        """
768        data = {"action": "results_templates"}
769        res, _ = self.sq("", data=data, reload=reload, bufver=bufver)
770        return sorted(tuple(a) for a in res)
771
772    def results_list(self, rtype, reload=False, bufver="0"):
773        """Return a list of available results for a result type template
774        `rtype` (a key from the `result_types` return value).
775
776        """
777        data = {"action": "results",
778                "results_templates_id": rtype}
779        data = self.add_auth(data)
780        res, legend = self.sq("", data=data, reload=reload,
[1619]781                              bufferkey=bufferkeypipax,
[1615]782                              bufver=bufver)
783        # index by unique_id (last column)
784        res = splitTableOnColumn(res, -1)
785        return dict((id, dict(zip(legend, line[0]))) \
786                    for id, line in res.items())
787
788    def download_key_function(self):
789        data = self.add_auth({"action": "download",
790                              "ids": "$MULTI$"}
791                                )
792        keynamingfn, _ = self.downloadMulti_bufcommand_replace_multi("",
[1619]793                         data=data, chunk=100, bufferkey=bufferkeypipax,
[1615]794                         transformfn=None)
795        return keynamingfn
796
797    def get_data(self, ids=None, result_type=None,
798                 exclude_constant_labels=False, average=median,
799                 callback=None, bufver="0", transform=None,
[1680]800                 allowed_labels=None, reload=False):
[1615]801        """
802        Get data in a single example table with labels of individual
803        attributes set to annotations for query and post-processing
804        instructions.
805
806        :param ids: If `result_type` is specified then this must
807            be a list of  `unique_id`s or (id, data_id) tuples as
808            returned by `mappings`. Else the `ids` must be a list
809            of `unique_id`s as returned by `results_list`.
810        :type ids: list
811
812        :param result_type: Result template type id as returned by
813             `result_types`. Not specified by default (see `ids`.)
814        :type result_type: str
815
816        :param exclude_constant_labels: If a label has the same value
817            in whole example table, remove it.
818        :type exclude_constant_labels: bool
819
820        :param average: Function used for combining multiple reading of
821            the same spot on a chip. If None, no averaging is done.
822            Function should take a list of floats and return an "averaged"
823            float (default `median`).
824        :type average: function
825
826        """
827
828        def optcb():
829            if callback:
830                callback()
831
[1620]832        def argsort(a):
833            sort = sorted([(map(int, item), i) for i, item in enumerate(a)])
834            return [i for _, i in sort]
835
[1615]836        cbc = CallBack(len(ids), optcb, callbacks=10)
837
838        if result_type is not None:
[1620]839            ids_sort = argsort(ids)
[1680]840            res_list = self.results_list(result_type, reload=reload,
841                                         bufver=bufver)
[1615]842            # Map (data_id, mapping_id) to unique_id
843            res_types_to_unique_id = \
844                dict(((annot["data_id"], annot["mappings_id"]),
845                      annot["unique_id"]) \
846                     for annot in res_list.values())
847
[1680]848            mappings = self.mappings(reload=reload, bufver=bufver)
[1615]849
850            def id_map(mappings_unique_id):
851                cbc()
852                if isinstance(mappings_unique_id, tuple):
853                    data_id, mappings_id = mappings_unique_id
854                else:
855                    annot = mappings[mappings_unique_id]
856                    data_id = annot["data_id"]
857                    mappings_id = annot["id"]
858                return res_types_to_unique_id[data_id, mappings_id]
859
860            ids = map(id_map, ids)
861        else:
862            result_type_set = set(id.rsplit("_", 1)[-1] for id in ids)
863            if len(result_type_set) != 1:
864                raise ValueError("""\
865Can only retrieve a single result_template_type at a time"""
866)
867            result_type = result_type_set.pop()
[1680]868            res_list = self.results_list(result_type, reload=reload,
869                                         bufver=bufver)
[1620]870            sort_keys = [(int(res_list[id]["data_id"]),
871                          int(res_list[id]["mappings_id"]))\
872                         for id in ids]
873            ids_sort = argsort(sort_keys)
874
875        # Sort the ids for use by pipax api
876        ids_sorted = [ids[i] for i in ids_sort]
[1615]877
878        read = {}
879        for a, b in res_list.items():
880            read[a] = b.items()
881
882        cbc.end()
883
[1680]884        download_func = lambda x: self.download(x, reload=reload,
885                                                bufver=bufver)
[1615]886
[1620]887        cbc = CallBack(len(ids_sorted) + 3, optcb,
[1615]888                       callbacks=99 - 20)
[1620]889        et = self.exampleTables(ids_sorted, spotmap={}, callback=cbc,
890                            annots=read,
[1615]891                            exclude_constant_labels=exclude_constant_labels,
892                            chipfn=download_func,
893                            allowed_labels=allowed_labels)
894        cbc.end()
895
[1620]896        # Restore input ids order.
897        domain = orange.Domain([et.domain[id] for id in ids], None)
898        domain.addmetas(et.domain.getmetas())
899        et = orange.ExampleTable(domain, et)
900
[1615]901        cbc = CallBack(2, optcb, callbacks=10)
902
903        #transformation is performed prior to averaging
904        if transform != None:
905            transformValues(et, fn=transform)  # in place transform
906            cbc()
907
908        #if average function is given, use it to join same spotids
909        if average != None:
910            et = averageAttributes(et, fn=average)
911            cbc()
912
913        cbc.end()
914
915        return et
916
917    def download(self, result_ids, reload=False, bufver="0"):
[1620]918        """result_ids must be sorted (by `(data_id, mappings_id)`)
919        """
[1615]920        data = {"action": "download", "ids": "$MULTI$"}
921        data = self.add_auth(data)
922        antss = self.downloadMulti("", result_ids, data=data, chunk=10,
[1619]923                                   bufferkey=bufferkeypipax, bufreload=reload,
[1615]924                                   bufver=bufver)
925        for a, legend in antss:
926            yield a
927
928    def downloadMulti(self, command, ids, data=None, chunk=100,
929                      transformfn=None, bufferkey=None, separatefn=None,
930                      bufreload=False, bufver="0"):
931        """
932        Downloads multiple results at once.
933        Results in the same order as in ids.
934
935        Bufferkey transforms command and data into buffer key.
936        bufver is a function returning buffer version for a given id. if
937            a string is given, use it for all ids
938        """
939
940        sids = split(ids, chunk)
941
942        bufverfn = None
943        if isinstance(bufver, basestring):
944            bufverfn = lambda x: bufver
945        else:
946            bufverfn = bufver
947
948        bufcommand, replace_multi = \
949                self.downloadMulti_bufcommand_replace_multi(command,
950                                data=data, chunk=chunk, bufferkey=bufferkey,
951                                transformfn=transformfn
952                                )
953
954        for i, sidp in enumerate(sids):
955
956            buffered = []
957            unbuffered = []
958
959            for a in sidp:
960                if self.inBuffer(bufcommand(a)) == bufverfn(a) and \
961                        bufreload == False:
962                    buffered.append(a)
963                else:
964                    unbuffered.append(a)
965
966            res = []
967            legend = []
968
969            if len(unbuffered) > 0:
970                com1, d1 = replace_multi(command, data, ",".join(unbuffered))
971                # get unbuffered part
972                res, legend = self.sq(com1, data=d1, buffer=False)
973            else:
974                # get legend from buffer also
975                legend = self.fromBuffer(bufcommand(buffered[0]))[0]
976
977            legend = ["gene_id", "value"]
[1619]978            genes = nth(res, 0)
[1615]979
980            antss = {}
981            for i, cid in enumerate(unbuffered):
982                col = i + 1
983                vals = nth(res, col)
984                antss[cid] = [[a, b] for a, b in zip(genes, vals) if b != "?"]
985
986            #here save buffer
987            for a, b in antss.items():
988                self.toBuffer(bufcommand(a), [legend] + b, bufverfn(a),
989                              autocommit=False)
990            self.bufferCommit()
991
992            #get buffered from the buffer
993            antssb = dict([(b, self.fromBuffer(bufcommand(b))[1:]) \
994                           for b in buffered])
995            antss.update(antssb)
996
997            #put results in order
998            for ci in sidp:
999                yield antss[ci], legend
1000
1001
[1021]1002class DictyExpress(DBCommon):
[420]1003    """
1004    Type is object id
1005    """
1006   
1007    aoidPairs = txt2ll("""time extractions.developmental_time_point
1008sample biological_samples.sample
1009growthCond biological_samples.growth_condition
1010treatment biological_samples.treatment
1011replicate biological_sample_replicates.replicate
1012techReplicate chips.replicate
1013platform chips.chip_platform
1014isTimeSeries biological_sample_replicates.is_time_series""")
1015
1016    obidPairs = txt2ll("""norms normalizations
1017samples biological_samples
1018replicates biological_sample_replicates
1019analysis analysis
1020experiments experiments
1021extractions extractions
1022chips chips""")
1023
[927]1024    def __init__(self, address=defaddress, buffer=None):
[670]1025        self.address = address
[420]1026        self.db = DBInterface(address)
[527]1027        self.buffer = buffer
[420]1028        self.preload()
1029
1030    def preload(self):
1031
1032        # aoids are mappings from annotation name to annotation id
1033        self.aoids = ll2dic(self.__annotationTypes(), 1, 0)
1034        self.saoids = ll2dic(self.aoidPairs, 0, 1)
1035        self.aoidsr = reverseDic(self.aoids)
1036        self.saoidsr = reverseDic(self.saoids)
1037
1038        # obids are mappings from object id to annotation id
1039        self.obids = ll2dic(self.__objects(), 1, 0)
1040        self.sobids = ll2dic(self.obidPairs, 0, 1)
1041        self.obidsr = reverseDic(self.obids)
1042        self.sobidsr = reverseDic(self.sobids)
1043
1044    def aoidt(self, s):
1045        return chainLookup(s, [self.saoids, self.aoids], force=[-1])
1046
1047    def obidt(self, s):
1048        return chainLookup(s, [self.sobids, self.obids], force=[-1])
1049 
1050    def aoidtr(self, s, **kwargs):
1051        return chainLookup(s, [self.aoidsr, self.saoidsr], **kwargs)
1052
1053    def obidtr(self, s):
1054        return chainLookup(s, [self.obidsr, self.sobidsr])
1055
1056    def pq(self, q):
1057        """
1058        Prepare query.
1059        ||| separator between conditions,
1060        *** denotes equivalence
1061        """
1062        o =  "|||".join([ self.aoidt(a) + "***" + b for a,b in q.items()])
1063        return o
1064
1065    def geneInfo(self):
1066        res,legend = self.sq("action=gene_info")
1067        return res, legend
1068
[527]1069    def annotationOptions(self, ao=None, onlyDiff=False, **kwargs):
[420]1070        """
[527]1071        Returns annotation options for given query. Returns all possible
1072        annotations if the query is omitted.
1073
1074        If ao is choosen, only result
[420]1075        """
1076        params = ""
1077        if len(kwargs) > 0: params += "&query=" + self.pq(kwargs)
[527]1078        if ao: params += "&annotation_object_id=" +  self.aoidt(ao)
[670]1079        res,legend = self.sq("action=get_annotation_options%s" % (params), bufadd=self.address)
[420]1080        res = onlyColumns(res, legend, ['annotation_object_id', 'value' ])
1081
1082        #join columns with the annotation object id
1083        joined = {}
1084        for key,v in res:
1085            key = self.aoidtr(key)
1086            cur = joined.get(key, [])
1087            cur.append(v)
1088            joined[key] = cur
1089
1090        if onlyDiff:
1091            joined = dict([ (a,b) for a,b in joined.items() if len(b)>1 ])
1092
1093        return dict([ (a, sorted(b)) for a,b in joined.items() ])
1094
1095    def annotation(self, type, id):
1096        return list(self.annotations(type, [ id ]))[0]
1097
1098    def meaningfulAnnot(self, name):
1099        if name in self.saoids:
1100            return True
1101        else:
1102            return False
1103
1104    def keepOnlyMeaningful(self, annot):
1105        """
1106        Keep only meaningful annotations
1107        """
1108        if type(annot) == type({}):
1109            return dict( [ (a,b) for a,b in annot.items() \
1110                if self.meaningfulAnnot(a) ] )
1111        else:
1112            return [ [ a,b ] for a,b in annot \
1113                if self.meaningfulAnnot(a) ]
1114
1115
[527]1116    def annotations(self, type, ids=None, all=False):
[557]1117        """
[610]1118        Returns a generator returning annotations for specified type and ids.
1119        If ids are left blank, all annotations are outputed. Annotations are in the same order
1120        as input ids.
[705]1121        If all is True, all annotations are kept, else keep only "meaningful".
[557]1122        """
1123       
1124        inputids = False
1125        if ids != None:
1126            inputids = True
[558]1127            antss = self.downloadMulti(
1128                "action=get_annotations&ids=$MULTI$&object_id=%s" 
[557]1129                % (self.obidt(type)), ids)
1130        else:
1131            res,legend = self.sq(
1132                "action=get_annotations&object_id=%s"
1133                % (self.obidt(type)))
1134            antss = splitTableOnColumn(res, 0)
1135            ids = nth(antss.items(),0)
1136            antss = zip(nth(antss.items(),1), [ legend ]*len(antss))
[420]1137
[558]1138        for ants in izip(antss,ids):
1139            (res, legend), id = ants
[420]1140            res2 = onlyColumns(res, legend, ['name', 'value'])
1141            res2 = [ [ self.aoidtr(a),b ] for a,b in res2 ]
1142            if not all:
1143                res2 = self.keepOnlyMeaningful(res2)
[557]1144            if inputids:
1145                yield res2
1146            else:
1147                yield (id, res2)
[420]1148
1149    def search(self, type, **kwargs):
1150        """
1151        Break search for multiple values of one attribute to independant searches.
[527]1152        Search is case insensitive.
[420]1153       
1154        List of searchable annotation types: self.saoids.keys()
1155
1156        example usage:
1157        search("norms", platform='minichip', sample='abcC3-')
1158            finds all ids of normalized entries where platform is minchip and
1159            sample is abcC3-
1160        search("norms", platform='minichip', sample=[ 'abcC3-', 'abcG15-'])
1161            finds all ids of normalized entries where platform is minichip and
1162            sample is abcC3- or those where platform is minichip and sample
1163            is abcG15-
1164        """
1165       
1166        #search for all combinations of values - this is slow!
1167
1168        l = []
1169        for k,v in kwargs.items():
[527]1170            if not issequencens(v):
[420]1171                v = [ v ]
1172            l.append((k,v))
1173
1174        ares = []
1175
1176        for r in mxrange([ len(v) for k,v in l ]):
1177            dico = {}
1178            for i,a in enumerate(r):
1179                dico[l[i][0]] = l[i][1][a]
1180
1181            res,_ = self.sq("action=search&object_id=%s&query=%s" \
[670]1182                % (self.obidt(type), self.pq(dico)), bufadd=self.address)
[527]1183
[420]1184            ares += res
1185
1186        return sorted(set(nth(ares, 0)))
1187
[558]1188
[527]1189    def chipN(self, id):
[558]1190        return list(self.chipNs([id]))[0]
[420]1191
1192    def chipR(self, id):
[558]1193        return list(self.chipRs([id]))[0]
[420]1194 
[888]1195    def chipNs(self, ids, remove_bad=True):
1196        """
1197        If removebad = True removes those with weights 0.
1198        """
[558]1199         
1200        def sel(res, legend):
1201            #Drop unwanted columns - for efficiency
[888]1202            res = onlyColumns(res, legend, ["spot_id", 'M', 'weights'])
1203            legend = onlyColumns( [ legend ], legend, ["spot_id", 'M', 'weights'])[0]
1204
1205            def transf(a):
1206                if a[2] == 0:
1207                    return a[0], '', a[2]
1208                else:
1209                    return a[0], a[1], a[2]
1210
1211            res = [ transf(a) for a in res ]
[558]1212            res = onlyColumns(res, legend, ["spot_id", 'M'])
1213            legend = onlyColumns( [ legend ], legend, ["spot_id", 'M'])[0]
1214            return res, legend
1215
1216        antss = self.downloadMulti("action=get_normalized_data&ids=$MULTI$", ids, chunk=2, transformfn=sel)
1217        for a,legend in antss:
1218            yield a   
[420]1219
[605]1220    def chipNsN(self, ids, annots):
[793]1221        """
1222        Download chips using new shorter format.
1223        """
[605]1224        chip_map_ids = zip(ids,[ dict(a)['chips.chip_map_id'] for a in annots ])
1225
1226        def separateByChipsMaps(l):
1227            begin = 0
1228            cm = l[0][1]
1229            cp = 0
1230            for id,m in l[1:]:
1231                cp += 1
1232                if m != cm:
1233                    yield l[begin:cp]
1234                    cm = m
1235                    begin = cp
1236            yield l[begin:cp+1]
1237       
1238        sep = list(separateByChipsMaps(chip_map_ids))
[670]1239     
[605]1240        def sel(res, legend):
1241            #Drop unwanted columns - for efficiency
1242            res = onlyColumns(res, legend, ["spot_id", 'M'])
1243            legend = onlyColumns( [ legend ], legend, ["spot_id", 'M'])[0]
1244            return res, legend
1245
1246        def separatefn(res):
[670]1247            #each one is own rown
1248            #genes are in the first row
1249            genes = res[0][1:]
1250            cids = nth(res,0)[1:]
1251
[605]1252            antss = {}
1253            for i,cid in enumerate(cids):
[670]1254                row = i+1
1255                vals = res[row][1:]
[605]1256                antss[cid] = [ list(a) for a in zip(genes, vals) ]
1257            return ['spot_id', 'M'], antss
1258
1259        for part in sep:
[670]1260            pids = nth(part,0)
[605]1261            antss = self.downloadMulti("action=get_normalized_data&mergeexperiments=1&ids=$MULTI$", pids, chunk=10, transformfn=sel, separatefn=separatefn)
1262            for a, legend in antss:
1263                yield a
1264
[420]1265    def chipRs(self, id):
[558]1266        antss = self.downloadMulti("action=get_raw_data&ids=$MULTI$", ids, chunk=2)
1267        for a,legend in antss:
1268            yield a
1269 
[420]1270    def spotId(self):
1271        res,legend = self.sq("action=spot_id_mapping")
[831]1272        res2 = onlyColumns(res, legend, ["spot_id", 'ddb_g', 'genename'])
[420]1273        return res2
1274
1275    def annotationTypes(self):
1276        """
1277        Returns all annotation types.
1278        """
1279        return self.aoids.keys()
1280
1281    def objects(self):
1282        """
1283        Returns all objects.
1284        """
1285        return self.obids.keys()
1286
1287    def __annotationTypes(self):
1288        """
1289        Returns list of [ annotation_object_id, name ]
1290        """
1291        res, legend = self.sq("action=get_annotation_types")
1292        res2 = onlyColumns(res, legend, ["annotation_object_id", 'name'])
1293        return res2
1294
1295    def __objects(self):
1296        res, legend = self.sq("action=get_objects")
1297        res2 = onlyColumns(res, legend, ["object_id", 'object_name'])
1298        return res2
1299
1300    def spotMap(self):
1301        spotids = self.spotId()
1302        spotmap = [ (a[0],a[1]) for a in spotids ]
1303
1304        spotmapd = {}
1305
1306        for a,b in spotmap:
1307            if a in spotmapd:
1308                spotmapd[a] = spotmapd[a] + "-" + b
1309            else:
1310                spotmapd[a] = b
1311
1312        return spotmapd
1313
[1017]1314    def getData(self, *args, **kwargs):
1315        deprecatedError("Use get_single_data instead")
1316
[1021]1317    def get_data(self, type="norms", exclude_constant_labels=False, average=median, 
[1154]1318        ids=None, callback=None, format="short", transform=None, allowed_labels=None, **kwargs):
[420]1319        """
[1017]1320        Get data in a single example table with labels of individual attributes
1321        set to annotations for query and post-processing
[420]1322        instructions.
[610]1323
[420]1324        Parameters:
1325            average: function used for combining multiple reading of the same spot on
[610]1326                a chip. If None, no averaging is done. Fuction should take a list
1327                of floats and return an "averaged" float.
[1017]1328            ids: a list of chip ids. If absent, make a search
[838]1329            exclude_constant_labels: if a label has the same value in whole
[888]1330                example table, remove it
1331            format: if short, use short format for chip download
[420]1332
[1017]1333        Defaults: Median averaging.
[420]1334        """
[557]1335
1336        def optcb():
1337            if callback: callback()
1338
[794]1339        cbc = CallBack(1, optcb, callbacks=10)
1340
[527]1341        if not ids:
1342            #returns ids of elements that match the search function
1343            ids = self.search(type, **kwargs)
[420]1344
[557]1345        cbc.end()
1346
[420]1347        #downloads annotations
[557]1348        cbc = CallBack(len(ids), optcb, callbacks=10)
[605]1349
[794]1350        readall = self.dictionarize(ids, self.annotations, type, ids, all=True, callback=cbc)
[605]1351
1352        read = {}
1353        for a,b in readall.items():
1354            read[a] = self.keepOnlyMeaningful(b)
1355
[669]1356        annotsinlist = [] #annotations in the same order
[605]1357        for id in ids:
1358            annotsinlist.append(readall[id])
1359
[680]1360        if verbose:
1361            print zip(ids,[ dict(a)['chips.chip_map_id'] for a in annotsinlist ])
[670]1362
[557]1363        cbc.end()
[527]1364
[1017]1365        #till now downloads were small
[420]1366
[527]1367        import time
[1017]1368        tstart = time.time()
[527]1369
[557]1370        #here download actually happens
[894]1371        chipfn = None
[527]1372
[605]1373        if type == "norms":
[888]1374            if format == "short":
[1017]1375                chipfn = lambda x, al=annotsinlist: self.chipNsN(x, al)
[888]1376            else:
[894]1377                chipfn = self.chipNs
[605]1378        else:
[894]1379            chipfn = self.chipRs
[557]1380       
[670]1381        if verbose:
1382            print "DOWNLOAD TIME", time.time() - tstart
[420]1383
[1020]1384        cbc = CallBack(len(ids)*2+len(ids)+1, optcb, callbacks=999-30)
[1154]1385        et = self.exampleTables(ids, spotmap=self.spotMap(), callback=cbc, annots=read, exclude_constant_labels=exclude_constant_labels, chipfn=chipfn, allowed_labels=allowed_labels)
[557]1386        cbc.end()
[527]1387
[1017]1388        cbc = CallBack(1, optcb, callbacks=10)
[557]1389
[1148]1390        #transformation is performed prior to averaging
1391        if transform != None:
1392            transformValues(et, fn=transform) #in place transform
1393            cbc()
1394
[420]1395        #if average function is given, use it to join same spotids
1396        if average != None:
[1017]1397            et = averageAttributes(et, fn=average)
1398            cbc()
[420]1399
[557]1400        cbc.end()
1401
[1017]1402        return et
[1021]1403   
1404    def get_single_data(self, *args, **kwargs):
1405        return self.get_data(*args, **kwargs)
[420]1406
[1021]1407class DatabaseConnection(DictyExpress):
1408    pass
[794]1409
[838]1410def allAnnotationVals(annots):
1411    """
1412    All annotation valuess for given annotations
1413    in a dict of { name: set of possible values } pairs.
1414    """
[1021]1415    av = defaultdict(set)
[838]1416    for a in annots:
1417        for name,val in a:
[1021]1418            av[name].add(val)
[838]1419    return av
1420
1421def createExampleTable(names, vals, annots, ddb, cname="DDB", \
[1154]1422        exclude_constant_labels=False, permutation=None, always_include=["id"], allowed_labels=None):
[420]1423    """
1424    Create an ExampleTable for this group. Attributes are those in
1425    names.
1426    """
1427    attributes = [ orange.FloatVariable(n, numberOfDecimals=3) \
1428        for n in names ]
[794]1429
[838]1430    #exclusion of names with constant values
1431    annotsvals = allAnnotationVals(annots)
1432    oknames = set(annotsvals.keys())
1433    if exclude_constant_labels:
[1074]1434        oknames = set(nth(filter(lambda x: len(x[1]) > 1 or x[0] in always_include, 
[838]1435            annotsvals.items()), 0))
1436
[1154]1437    if allowed_labels != None:
1438        oknames = set(filter(lambda x: x in allowed_labels, oknames))
1439
[869]1440    #print oknames
[838]1441
[794]1442    for a,an in zip(attributes, annots):
[1371]1443        a.setattr("attributes", dict([(name,str(val)) for name,val in an if name in oknames]))
[794]1444
[420]1445    domain = orange.Domain(attributes, False)
1446    ddbv = orange.StringVariable(cname)
1447    id = orange.newmetaid()
1448    domain.addmeta(id, ddbv)
1449
[894]1450    examples = [ None ]*len(ddb)
1451    for i,(v,d) in enumerate(izip(izip(*vals), ddb)):
[420]1452        ex = orange.Example(domain, [ floatOrUnknown(a) for a in v ])
1453        ex[cname] = d
[894]1454
1455        if permutation:
1456            examples[permutation[i]] = ex
1457        else:
1458            examples[i] = ex
[420]1459
1460    return orange.ExampleTable(domain,examples)
1461
[1148]1462def transformValues(data, fn):
1463    """
1464    In place transformation.
1465    """
1466    for ex in data:
1467        for at in data.domain.attributes:
1468            if not ex[at].isSpecial():
1469                ex[at] = fn(ex[at])
[420]1470
1471def averageAttributes(data, joinc="DDB", fn=median):
1472    """
1473    Averages attributes with the same "join" parameter using
1474    specified function. Builds a now dataset. Order of
1475    "join" parameter stays the same with regard to first
1476    appearance.
1477    """
1478
[527]1479    if verbose:
1480        print "Averaging attributes"
1481
[420]1482    valueso = []
1483    valuess = set(valueso)
1484
1485    attributes = [ a for a in data.domain.attributes ]
1486    domain = data.domain
[1359]1487   
[420]1488    #accumulate values
1489
1490    for ex in data:
1491        j = str(ex[joinc])
1492
1493        if j not in valuess:
1494            valueso.append(j)
1495            valuess.add(j)
1496
[1359]1497    valuesddb = [ dict( (n,[]) for n in valueso) for at in attributes ]
1498   
1499    for ex in data:
1500        j = str(ex[joinc])
1501
1502        for a in range(len(attributes)):
[420]1503            val = ex[a]
1504            if not val.isSpecial():
[1359]1505                valuesddb[a][j].append(val.native())
[420]1506
[1359]1507    #create a new example table reusing the domain - apply fn
[420]1508    examples = []
1509    for v in valueso:
1510        example = orange.Example(domain, \
[1359]1511            [ floatOrUnknown(fn(valuesddb[a][v])) for a in range(len(attributes)) ] )
[420]1512        example[joinc] = v
1513        examples.append(example)
1514
1515    return orange.ExampleTable(domain, examples)
1516
1517def floatOrUnknown(a):
1518    """
1519    Converts an element to float if possible.
1520    If now, output "?".
1521    """
1522    try:
1523        return float(a)
1524    except:
1525        return "?"
1526
[794]1527class CallBack():
1528    """
1529    Converts "allparts" callbacks into by "callbacks"
1530    specified number of callbacks of function fn.
1531    """
1532
1533    def __init__(self, allparts, fn, callbacks=100):
1534        self.allparts = allparts
1535        self.lastreport = 0.00001
1536        self.getparts = 0
1537        self.increase = 1.0/callbacks
1538        self.callbacks = callbacks
1539        self.fn = fn
1540        self.cbs = 0
1541
1542    def __call__(self):
1543        self.getparts += 1
1544        done = float(self.getparts)/self.allparts
1545        while done > self.lastreport + self.increase:
1546            self.lastreport += self.increase
1547            self.fn()
1548            self.cbs += 1
1549
1550    def end(self):
1551        while self.cbs < self.callbacks:
1552            self.fn()
1553            self.cbs += 1
1554
[558]1555class BufferSQLite(object):
[527]1556
1557    def __init__(self, filename, compress=True):
[792]1558        self.compress = compress
1559        self.filename = filename
1560        self.conn = self.connect()
1561
1562    def clear(self):
1563        """
1564        Removes all entries in the buffer
1565        """
1566        self.conn.close()
1567        os.remove(self.filename)
1568        self.conn = self.connect()
1569
1570    def connect(self):
[527]1571        import sqlite3
[792]1572        conn = sqlite3.connect(self.filename)
1573        c = conn.cursor()
[527]1574        c.execute('''create table if not exists buf
1575        (address text primary key, time text, con blob)''')
1576        c.close()
[792]1577        conn.commit()
1578        return conn
[527]1579
1580    def contains(self, addr):
[1049]1581        """ Returns version or False, if it does not exists """
[527]1582        c = self.conn.cursor()
[1049]1583        c.execute('select time from buf where address=?', (addr,))
1584        lc = list(c)
[527]1585        c.close()
[1049]1586        if len(lc) == 0:
[527]1587            return False
1588        else:
[1049]1589            return lc[0][0]
[527]1590
[1049]1591    def list(self):
1592        c = self.conn.cursor()
1593        c.execute('select address from buf')
1594        return nth(list(c), 0)
1595
[1178]1596    def add(self, addr, con, version="0", autocommit=True):
[527]1597        import cPickle, zlib, sqlite3
[558]1598        if verbose:
1599            print "Adding", addr
[527]1600        c = self.conn.cursor()
1601        if self.compress:
1602            bin = sqlite3.Binary(zlib.compress(cPickle.dumps(con)))
1603        else:
1604            bin = sqlite3.Binary(cPickle.dumps(con))
[1049]1605        c.execute('insert or replace into buf values (?,?,?)', (addr, version, bin))
[527]1606        c.close()
[1178]1607        if autocommit:
1608            self.commit()
1609
1610    def commit(self):
[527]1611        self.conn.commit()
1612
1613    def get(self, addr):
1614        import cPickle, zlib
1615        if verbose:
1616            print "getting from buffer", addr
1617            t = time.time()
1618        c = self.conn.cursor()
1619        c.execute('select con from buf where address=?', (addr,))
1620        ls = list(c)
1621        first = ls[0][0]
1622        if verbose:
1623            print time.time() - t
1624        if self.compress:
1625            rc = cPickle.loads(zlib.decompress(first))
1626        else:
1627            rc = cPickle.loads(str(first))
1628        c.close()
1629
1630        if verbose:
1631            print time.time() - t
1632        return rc
1633
[777]1634def download_url(url, repeat=2):
1635    def do():
1636        return urllib2.urlopen(url)
1637
1638    if repeat <= 0:
1639        do()
1640    else:
1641        try:
1642            return do()
1643        except:
1644            return download_url(url, repeat=repeat-1)
1645
1646def empty_none(s):
1647    if s:
1648        return s
1649    else:
1650        return None
1651
[1036]1652def join_ats(atts):
1653    """ Joins attribute attributes together. If all values are the same,
1654    set the parameter to the common value, else return a list of the
1655    values in the same order as the attributes are imputed. """
1656    keys = reduce(lambda x,y: x | y, (set(at.keys()) for at in atts))
1657    od = {}
1658
1659    def man(x):
1660        if issequencens(x):
1661            return tuple(x)
1662        else:
1663            return x
1664
1665    for k in keys:
1666        s = set(man(at[k]) for at in atts)
1667        if len(s) == 1:
1668            od[k] = list(s)[0]
1669        else:
[1371]1670            od[k] = str([ at[k] for at in atts ])
[1036]1671    return od
1672
[1480]1673def join_replicates(data, ignorenames=["replicate", "id", "name", "map_stop1"], namefn=None, avg=median):
[1036]1674    """ Join replicates by median.
[1480]1675    Default parameters work for PIPA data.
1676    Sort elements in the same order as ignorenames!
1677    """
1678
[1036]1679    d = defaultdict(list)
1680
1681    if namefn == None:
[1049]1682        namefn = lambda att: ",".join(att["id"]) if issequencens(att["id"]) else att["id"]
[1036]1683
1684    #key function
1685    def key_g(att):
1686        dk = att.copy()
1687        for iname in ignorenames:
1688            dk.pop(iname, None)
1689       
1690        def man(x):
1691            if issequencens(x):
1692                return tuple(x)
1693            else:
1694                return x
1695
1696        return tuple(nth(sorted(((a, man(b)) for a,b in dk.items())), 1))
1697
1698    #prepare groups
1699    for i,a in enumerate(data.domain.attributes):
1700        att = a.attributes
1701        k = key_g(att)
1702        d[k].append(i)
1703
1704    d = dict(d) #want errors with wrong keys
1705
1706    natts = []
1707
1708    def nativeOrNone(val):
1709        if val.isSpecial(): 
1710            return None
1711        else: 
1712            return val.native()
1713
1714    def avgnone(l):
1715        """ Removes None and run avg function"""
1716        l = filter(lambda x: x != None, l)
1717        if len(l):
1718            return avg(l)
1719        else:
1720            return None
1721
[1480]1722    def data_type(vals): #data type converter
1723        try:
1724            _ = [ int(a) for a in vals ]
1725            return int
1726        except:
1727            try:
1728                _ = [ float(a) for a in vals ]
1729                return float
1730            except:
1731                return lambda x: x
1732
1733    all_values = defaultdict(set)
1734    for a in [ at.attributes for at in data.domain.attributes ]:
1735        for k,v in a.iteritems():
1736            all_values[k].add(v)
1737
1738    types = {}
1739    for k,vals in all_values.iteritems():
1740        types[k] = data_type(vals)
[1482]1741
[1036]1742    for group, elements in d.items():
1743        a = orange.FloatVariable()
[1480]1744        #here sort elements
1745   
1746        def sk(x):
[1482]1747            return [ types[n](x[n]) for n in ignorenames if n in all_values ]
[1480]1748
1749        elements = sorted(elements, key=lambda x: sk(data.domain.attributes[x].attributes))
1750
[1037]1751        a.attributes.update(join_ats([data.domain.attributes[i].attributes for i in elements]))
[1036]1752        a.name = namefn(a.attributes)
1753
1754        def avgel(ex, el):
1755            return orange.Value(avgnone([ nativeOrNone(ex[ind]) for ind in el ]))
1756
1757        a.getValueFrom = lambda ex,rw,el=elements: avgel(ex,el)
1758        natts.append(a)
1759
1760    ndom = orange.Domain(natts, data.domain.classVar)
1761    ndom.addmetas(data.domain.getmetas())
1762    return orange.ExampleTable(ndom, data)
1763
1764
[777]1765class DictyBase(object):
1766
1767    domain = "dictybase"
1768    filename = "information_mappings.pck"
[1828]1769    tags = [ "Dictyostelium discoideum", "gene", "essential", "dictyBase"] 
[777]1770 
1771    @classmethod
1772    def version(cls):
[832]1773        orngServerFiles.localpath_download(cls.domain, cls.filename)
[779]1774        return orngServerFiles.info(cls.domain, cls.filename)["datetime"]
[777]1775   
1776    @classmethod
1777    def download_information(cls):
1778        """
1779        Downloads gene information and parses it.
1780        Returns a dictionary {ID: (name, synonyms, products)}
1781        """
1782        s = download_url("http://www.dictybase.org/db/cgi-bin/dictyBase/download/download.pl?area=general&ID=gene_information.txt").read()
1783        out = []
1784        for l in txt2ll(s, separ='\t', lineSepar='\n')[1:]:
1785            if len(l) == 4:
1786                id = l[0]
1787                name = l[1]
1788                synonyms = filter(None, l[2].split(", "))
1789                products = l[3]
1790                out.append((id, name, synonyms, products))
1791        return dict((a,(b,c,d)) for a,b,c,d in out)
1792
1793    @classmethod
1794    def download_mappings(cls):
1795        """
1796        Downloads DDB-GeneID-UniProt mappings and parses them.
1797        Returns a list of (ddb, ddb_g, uniprot) triplets.
1798       
1799        2009/04/07: ddb's appear unique
1800        """
1801        s = download_url("http://www.dictybase.org/db/cgi-bin/dictyBase/download/download.pl?area=general&ID=DDB-GeneID-UniProt.txt").read()
1802        out = []
1803        for l in txt2ll(s, separ='\t', lineSepar='\n')[1:]:
[1335]1804            if len(l) == 4:
[777]1805                ddb = empty_none(l[0])
1806                ddb_g = empty_none(l[1])
[1335]1807                name = empty_none(l[2])
1808                uniprot = empty_none(l[3])
1809                out.append((ddb, ddb_g, name, uniprot))
[777]1810        return out
1811
1812    @classmethod
1813    def pickle_data(cls):
1814        info = cls.download_information()
1815        mappings = cls.download_mappings()
1816        return pickle.dumps((info,mappings), -1)
1817
1818    def __init__(self):
[781]1819        fn = orngServerFiles.localpath_download(self.domain, self.filename)
[779]1820        self.info, self.mappings = pickle.load(open(fn, 'rb'))
[777]1821
[557]1822if __name__=="__main__":
[1049]1823    verbose = 1
[777]1824
[1021]1825    def printet(et):
1826        et.save("ett.tab")
1827        print open("ett.tab").read()
1828
[1148]1829    """
[779]1830    a = DictyBase()
1831    print len(a.info)
1832
[1111]1833    dbc = DictyExpress(buffer=BufferSQLite("../tmpbufnew"))
[558]1834
[670]1835    print dbc.annotationOptions()
[527]1836
[557]1837    count = 0
1838    def cb():
1839        global count
1840        count += 1
[558]1841        #print "CBBB", count
[527]1842
[1155]1843    et = dbc.get_single_data(sample=[ "tagA-", "pkaC-"], callback=cb, exclude_constant_labels=True, allowed_labels=["sample"])
[838]1844    print et.domain
1845    print et.domain[0].attributes
[1021]1846    printet(et)
[1155]1847    """
[838]1848
[1292]1849    d = PIPA(buffer=BufferSQLite("../tmpbufnewpipa"))
[1029]1850
[1234]1851    print d.gene_expression_types()
1852
[1510]1853    cov = d.coverage("777", chromosome=1, genome="dd")
1854    print cov
1855
[1234]1856    allids = d.annotations().keys()
[1510]1857    allids = allids[:1]
[1292]1858    print ("list", allids)
[1234]1859    print d.annotations().items()[0]
1860    print ("annots", d.annotations().items()[:2])
[1049]1861
[1371]1862    allids = map(str, allids)
[1172]1863
[1148]1864    import math
[1021]1865
[1371]1866    data = d.get_data(ids=allids, ctype="3")
1867
1868    #data = d.get_data(ids=allids, transform=lambda x: math.log(x+1, 2), allowed_labels=["strain"], ctype="3")
[1234]1869
[1148]1870    printet(data)
Note: See TracBrowser for help on using the repository browser.