source: orange-bioinformatics/obiDicty.py @ 1292:e8f1724210d7

Revision 1292:e8f1724210d7, 45.4 KB checked in by markotoplak, 3 years ago (diff)

obiDicty: fixed for new PIPA interface.

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