Revision 1636:10d234fdadb9, 47.3 KB checked in by mitar, 2 years ago (diff)

Restructuring because we will not be using namespaces.

Line
1import sys, pprint, time, re
2from itertools import *
3import urllib
4import urllib2
5import socket
6import os
7from collections import defaultdict
8import pickle
9
10import orange
11from Orange.orng import orngServerFiles
12
17
20
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
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
130#end utility functions
131
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
143socket.setdefaulttimeout(60)
144
145verbose = 0
146
147class HttpGetException(Exception): pass
148
151
153    if verbose:
154        print address, args, kwargs,  " "
156    t1 = time.time()
157    f = urllib2.urlopen(address, *args, **kwargs)
159    if verbose:
161    if verbose:
162        print time.time() - t1
164
165def txt2ll(s, separ=' ', lineSepar='\n'):
166    return [ a.split(separ) for a in s.split(lineSepar) ]
167
168class AuthenticationError(Exception):
169    pass
170
171class DBInterface(object):
172
175
176    def raw(self, request, data=None, tryN=3):
177        if verbose:
178            print "tryN", tryN
179
180        def try_down():
181            try:
182                if data == None:
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:
194                return try_down()
195            except IOError:
196                return self.raw(request, data=data, tryN=tryN-1)
197        else:
198            return try_down()
199
200    def get(self, request, data=None, tryN=3):
201        rawf = self.raw(request, data)
202        if rawf == None:
203            print "rawf == None!"
204            raise Exception("Connection error when contacting " + self.address + request)
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
208            #TODO are there any other kinds of errors?
209            if tryN > 0:
210                if verbose:
211                    print "trying again"
212                return self.get(request, data=data, tryN=tryN-1)
213            else:
214                if verbose:
215                    print rafw[:1000]
216                raise Exception("Error with the database")
217
218        a = txt2ll(rawf, separ='\t')
219
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
287class DBCommon(object):
288
291
292    def toBuffer(self, addr, cont, version, autocommit=True):
293        if self.buffer:
295
296    def bufferCommit(self):
297        if self.buffer:
298            self.buffer.commit()
299
300    def bufferFun(self, bufkey, bufver, reload, fn, *args, **kwargs):
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        """
306        if self.inBuffer(bufkey) == bufver and reload == False:
307            res = self.fromBuffer(bufkey)
308        else:
309            res = fn(*args, **kwargs)
310            self.toBuffer(bufkey, res, bufver)
311        return res
312
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
316        if buffer:
317            bufkey = bufadd + (bufname if bufname != None else s1)
318            if bufferkey != None:
319                bufkey = bufferkey(bufkey, data)
320            res = self.bufferFun(bufkey, bufver, reload, db.get, s1, data=data)
321        else:
322            res = db.get(s1, data=data)
323        return res[1:],res[0]
324
326        if self.buffer:
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
344        """
345        Get function which gives buffer address for an id and a function
346        which replaces \$MULTI\$.
347        """
348
349        def bufferkey1(command, data):
350            if transformfn:
351                return "TRANS " + command
352            else:
353                return command
354
355        def replace_multi(command, data, repl):
356            return command.replace("\$MULTI\$", repl),\
357                dict((a,b.replace("\$MULTI\$", repl)) for a,b in sorted(data.items())) if data != None else None
358
359        if bufferkey == None:
360            bufferkey=bufferkey1
361
362        bufcommand = lambda x, c=command, d=data: bufferkey(*replace_multi(c, d, x))
363        return bufcommand, replace_multi
364
366        """
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
384
385        for i,sidp in enumerate(sids):
386
387            buffered = []
388            unbuffered = []
389
390            for a in sidp:
391                if self.inBuffer(bufcommand(a)) == bufverfn(a) and bufreload == False:
392                    buffered.append(a)
393                else:
394                    unbuffered.append(a)
395
396            res = []
397            legend = []
398
399            if len(unbuffered) > 0:
400                com1, d1 = replace_multi(command, data, ",".join(unbuffered))
401                res, legend = self.sq(com1, data=d1, buffer=False) #get unbuffered part
402            else:
403                # get legend from buffer also
404                legend = self.fromBuffer(bufcommand(buffered[0]))[0]
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():
425                self.toBuffer(bufcommand(a), [ legend ] + b, bufverfn(a), autocommit=False)
426            self.bufferCommit()
427
428
429            #get buffered from the buffer
430            antssb = dict([ (b, self.fromBuffer(bufcommand(b))[1:]) for b in buffered ])
431            antss.update(antssb)
432
433            #put results in order
434            tl = []
435            for ci in sidp:
436                yield antss[ci], legend
437
438    def exampleTables(self, ids, chipsm=None, spotmap={}, callback=None, exclude_constant_labels=False, annots={}, chipfn=None, allowed_labels=None):
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)),
535            groupannots, ddb, exclude_constant_labels=exclude_constant_labels, permutation=permutation, allowed_labels=allowed_labels)
536
537        if callback: callback()
538
539        return et
540
541def bufferkeypipa(command, data):
542    """ Do not save password to the buffer! """
543    command = command + " v5" #add version
544    if data != None:
545        data = data.copy()
546        if pipaparpass in data:
547            data.pop(pipaparpass)
548        return command + " " +  urllib.urlencode(sorted(data.items()))
549    else:
550        return command
551
552class PIPA(DBCommon):
553
558        self.buffer = buffer
563
566            return data
568        if data != None:
569            authdic.update(data)
570        return authdic
571
573        """
574        Returns a dictionary of (id: dictionary of annotations).
575        """
578        return dict( (sa[0], dict(zip(legend[1:], sa[1:]))) for sa in res )
579
580    def chips(self, ids, ctype, reload=False, bufver="0"):
581        def separatefn(res):
582            #each one is own rown
583            #genes are in the first row
584            #remove unknown
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:]
592                antss[cid] = [ list(a) for a in zip(genes, vals) if a[1] != "?" ]
593            return ['gene_id', 'value'], antss
594
597        if ctype != None:
600        for a,legend in antss:
601            yield a
602
603    def coverage(self, id_, genome, chromosome, reload=False, bufver="0", raw=False):
604        data = self.add_auth({"action": "coverage_strand", "chr_name": str(chromosome),
605            "pixel_size": "1", "genome": genome, "sample_id": str(id_), "raw": str(raw).lower()})
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!
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
632        return sorted(tuple(a) for a in res)
633
634    def chips_keynaming(self, ctype):
636        return keynamingfn
637
638    def get_data(self, exclude_constant_labels=False, average=median,
639        ids=None, callback=None, bufver="0", transform=None, ctype=None, allowed_labels=None):
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
653            ctype: expression type, from gene_expression_types
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:
664            return None
665
666        cbc.end()
667
669        cbc = CallBack(len(ids), optcb, callbacks=10)
670
672
676
677        cbc.end()
678
680
681        import time
682        tstart = time.time()
683
685        chipfn = None
686
687        chipfn = lambda x: self.chips(x, ctype, bufver=bufver)
688
689        if verbose:
691
692        cbc = CallBack(len(ids)*2+len(ids)+1, optcb, callbacks=999-30)
693        et = self.exampleTables(ids, spotmap={}, callback=cbc, annots=read, exclude_constant_labels=exclude_constant_labels, chipfn=chipfn, allowed_labels=allowed_labels)
694        cbc.end()
695
696        cbc = CallBack(1, optcb, callbacks=10)
697
698        #transformation is performed prior to averaging
699        if transform != None:
700            transformValues(et, fn=transform) #in place transform
701            cbc()
702
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
712class DictyExpress(DBCommon):
713    """
714    Type is object id
715    """
716
717    aoidPairs = txt2ll("""time extractions.developmental_time_point
718sample biological_samples.sample
719growthCond biological_samples.growth_condition
720treatment biological_samples.treatment
721replicate biological_sample_replicates.replicate
722techReplicate chips.replicate
723platform chips.chip_platform
724isTimeSeries biological_sample_replicates.is_time_series""")
725
726    obidPairs = txt2ll("""norms normalizations
727samples biological_samples
728replicates biological_sample_replicates
729analysis analysis
730experiments experiments
731extractions extractions
732chips chips""")
733
737        self.buffer = buffer
739
741
742        # aoids are mappings from annotation name to annotation id
743        self.aoids = ll2dic(self.__annotationTypes(), 1, 0)
744        self.saoids = ll2dic(self.aoidPairs, 0, 1)
745        self.aoidsr = reverseDic(self.aoids)
746        self.saoidsr = reverseDic(self.saoids)
747
748        # obids are mappings from object id to annotation id
749        self.obids = ll2dic(self.__objects(), 1, 0)
750        self.sobids = ll2dic(self.obidPairs, 0, 1)
751        self.obidsr = reverseDic(self.obids)
752        self.sobidsr = reverseDic(self.sobids)
753
754    def aoidt(self, s):
755        return chainLookup(s, [self.saoids, self.aoids], force=[-1])
756
757    def obidt(self, s):
758        return chainLookup(s, [self.sobids, self.obids], force=[-1])
759
760    def aoidtr(self, s, **kwargs):
761        return chainLookup(s, [self.aoidsr, self.saoidsr], **kwargs)
762
763    def obidtr(self, s):
764        return chainLookup(s, [self.obidsr, self.sobidsr])
765
766    def pq(self, q):
767        """
768        Prepare query.
769        ||| separator between conditions,
770        *** denotes equivalence
771        """
772        o =  "|||".join([ self.aoidt(a) + "***" + b for a,b in q.items()])
773        return o
774
775    def geneInfo(self):
776        res,legend = self.sq("action=gene_info")
777        return res, legend
778
779    def annotationOptions(self, ao=None, onlyDiff=False, **kwargs):
780        """
781        Returns annotation options for given query. Returns all possible
782        annotations if the query is omitted.
783
784        If ao is choosen, only result
785        """
786        params = ""
787        if len(kwargs) > 0: params += "&query=" + self.pq(kwargs)
788        if ao: params += "&annotation_object_id=" +  self.aoidt(ao)
790        res = onlyColumns(res, legend, ['annotation_object_id', 'value' ])
791
792        #join columns with the annotation object id
793        joined = {}
794        for key,v in res:
795            key = self.aoidtr(key)
796            cur = joined.get(key, [])
797            cur.append(v)
798            joined[key] = cur
799
800        if onlyDiff:
801            joined = dict([ (a,b) for a,b in joined.items() if len(b)>1 ])
802
803        return dict([ (a, sorted(b)) for a,b in joined.items() ])
804
805    def annotation(self, type, id):
806        return list(self.annotations(type, [ id ]))[0]
807
808    def meaningfulAnnot(self, name):
809        if name in self.saoids:
810            return True
811        else:
812            return False
813
814    def keepOnlyMeaningful(self, annot):
815        """
816        Keep only meaningful annotations
817        """
818        if type(annot) == type({}):
819            return dict( [ (a,b) for a,b in annot.items() \
820                if self.meaningfulAnnot(a) ] )
821        else:
822            return [ [ a,b ] for a,b in annot \
823                if self.meaningfulAnnot(a) ]
824
825
826    def annotations(self, type, ids=None, all=False):
827        """
828        Returns a generator returning annotations for specified type and ids.
829        If ids are left blank, all annotations are outputed. Annotations are in the same order
830        as input ids.
831        If all is True, all annotations are kept, else keep only "meaningful".
832        """
833
834        inputids = False
835        if ids != None:
836            inputids = True
838                "action=get_annotations&ids=\$MULTI\$&object_id=%s"
839                % (self.obidt(type)), ids)
840        else:
841            res,legend = self.sq(
842                "action=get_annotations&object_id=%s"
843                % (self.obidt(type)))
844            antss = splitTableOnColumn(res, 0)
845            ids = nth(antss.items(),0)
846            antss = zip(nth(antss.items(),1), [ legend ]*len(antss))
847
848        for ants in izip(antss,ids):
849            (res, legend), id = ants
850            res2 = onlyColumns(res, legend, ['name', 'value'])
851            res2 = [ [ self.aoidtr(a),b ] for a,b in res2 ]
852            if not all:
853                res2 = self.keepOnlyMeaningful(res2)
854            if inputids:
855                yield res2
856            else:
857                yield (id, res2)
858
859    def search(self, type, **kwargs):
860        """
861        Break search for multiple values of one attribute to independant searches.
862        Search is case insensitive.
863
864        List of searchable annotation types: self.saoids.keys()
865
866        example usage:
867        search("norms", platform='minichip', sample='abcC3-')
868            finds all ids of normalized entries where platform is minchip and
869            sample is abcC3-
870        search("norms", platform='minichip', sample=[ 'abcC3-', 'abcG15-'])
871            finds all ids of normalized entries where platform is minichip and
872            sample is abcC3- or those where platform is minichip and sample
873            is abcG15-
874        """
875
876        #search for all combinations of values - this is slow!
877
878        l = []
879        for k,v in kwargs.items():
880            if not issequencens(v):
881                v = [ v ]
882            l.append((k,v))
883
884        ares = []
885
886        for r in mxrange([ len(v) for k,v in l ]):
887            dico = {}
888            for i,a in enumerate(r):
889                dico[l[i][0]] = l[i][1][a]
890
891            res,_ = self.sq("action=search&object_id=%s&query=%s" \
893
894            ares += res
895
896        return sorted(set(nth(ares, 0)))
897
898
899    def chipN(self, id):
900        return list(self.chipNs([id]))[0]
901
902    def chipR(self, id):
903        return list(self.chipRs([id]))[0]
904
906        """
907        If removebad = True removes those with weights 0.
908        """
909
910        def sel(res, legend):
911            #Drop unwanted columns - for efficiency
912            res = onlyColumns(res, legend, ["spot_id", 'M', 'weights'])
913            legend = onlyColumns( [ legend ], legend, ["spot_id", 'M', 'weights'])[0]
914
915            def transf(a):
916                if a[2] == 0:
917                    return a[0], '', a[2]
918                else:
919                    return a[0], a[1], a[2]
920
921            res = [ transf(a) for a in res ]
922            res = onlyColumns(res, legend, ["spot_id", 'M'])
923            legend = onlyColumns( [ legend ], legend, ["spot_id", 'M'])[0]
924            return res, legend
925
927        for a,legend in antss:
928            yield a
929
930    def chipNsN(self, ids, annots):
931        """
933        """
934        chip_map_ids = zip(ids,[ dict(a)['chips.chip_map_id'] for a in annots ])
935
936        def separateByChipsMaps(l):
937            begin = 0
938            cm = l[0][1]
939            cp = 0
940            for id,m in l[1:]:
941                cp += 1
942                if m != cm:
943                    yield l[begin:cp]
944                    cm = m
945                    begin = cp
946            yield l[begin:cp+1]
947
948        sep = list(separateByChipsMaps(chip_map_ids))
949
950        def sel(res, legend):
951            #Drop unwanted columns - for efficiency
952            res = onlyColumns(res, legend, ["spot_id", 'M'])
953            legend = onlyColumns( [ legend ], legend, ["spot_id", 'M'])[0]
954            return res, legend
955
956        def separatefn(res):
957            #each one is own rown
958            #genes are in the first row
959            genes = res[0][1:]
960            cids = nth(res,0)[1:]
961
962            antss = {}
963            for i,cid in enumerate(cids):
964                row = i+1
965                vals = res[row][1:]
966                antss[cid] = [ list(a) for a in zip(genes, vals) ]
967            return ['spot_id', 'M'], antss
968
969        for part in sep:
970            pids = nth(part,0)
972            for a, legend in antss:
973                yield a
974
975    def chipRs(self, id):
977        for a,legend in antss:
978            yield a
979
980    def spotId(self):
981        res,legend = self.sq("action=spot_id_mapping")
982        res2 = onlyColumns(res, legend, ["spot_id", 'ddb_g', 'genename'])
983        return res2
984
985    def annotationTypes(self):
986        """
987        Returns all annotation types.
988        """
989        return self.aoids.keys()
990
991    def objects(self):
992        """
993        Returns all objects.
994        """
995        return self.obids.keys()
996
997    def __annotationTypes(self):
998        """
999        Returns list of [ annotation_object_id, name ]
1000        """
1001        res, legend = self.sq("action=get_annotation_types")
1002        res2 = onlyColumns(res, legend, ["annotation_object_id", 'name'])
1003        return res2
1004
1005    def __objects(self):
1006        res, legend = self.sq("action=get_objects")
1007        res2 = onlyColumns(res, legend, ["object_id", 'object_name'])
1008        return res2
1009
1010    def spotMap(self):
1011        spotids = self.spotId()
1012        spotmap = [ (a[0],a[1]) for a in spotids ]
1013
1014        spotmapd = {}
1015
1016        for a,b in spotmap:
1017            if a in spotmapd:
1018                spotmapd[a] = spotmapd[a] + "-" + b
1019            else:
1020                spotmapd[a] = b
1021
1022        return spotmapd
1023
1024    def getData(self, *args, **kwargs):
1026
1027    def get_data(self, type="norms", exclude_constant_labels=False, average=median,
1028        ids=None, callback=None, format="short", transform=None, allowed_labels=None, **kwargs):
1029        """
1030        Get data in a single example table with labels of individual attributes
1031        set to annotations for query and post-processing
1032        instructions.
1033
1034        Parameters:
1035            average: function used for combining multiple reading of the same spot on
1036                a chip. If None, no averaging is done. Fuction should take a list
1037                of floats and return an "averaged" float.
1038            ids: a list of chip ids. If absent, make a search
1039            exclude_constant_labels: if a label has the same value in whole
1040                example table, remove it
1042
1043        Defaults: Median averaging.
1044        """
1045
1046        def optcb():
1047            if callback: callback()
1048
1049        cbc = CallBack(1, optcb, callbacks=10)
1050
1051        if not ids:
1052            #returns ids of elements that match the search function
1053            ids = self.search(type, **kwargs)
1054
1055        cbc.end()
1056
1058        cbc = CallBack(len(ids), optcb, callbacks=10)
1059
1060        readall = self.dictionarize(ids, self.annotations, type, ids, all=True, callback=cbc)
1061
1065
1066        annotsinlist = [] #annotations in the same order
1067        for id in ids:
1069
1070        if verbose:
1071            print zip(ids,[ dict(a)['chips.chip_map_id'] for a in annotsinlist ])
1072
1073        cbc.end()
1074
1076
1077        import time
1078        tstart = time.time()
1079
1081        chipfn = None
1082
1083        if type == "norms":
1084            if format == "short":
1085                chipfn = lambda x, al=annotsinlist: self.chipNsN(x, al)
1086            else:
1087                chipfn = self.chipNs
1088        else:
1089            chipfn = self.chipRs
1090
1091        if verbose:
1093
1094        cbc = CallBack(len(ids)*2+len(ids)+1, optcb, callbacks=999-30)
1095        et = self.exampleTables(ids, spotmap=self.spotMap(), callback=cbc, annots=read, exclude_constant_labels=exclude_constant_labels, chipfn=chipfn, allowed_labels=allowed_labels)
1096        cbc.end()
1097
1098        cbc = CallBack(1, optcb, callbacks=10)
1099
1100        #transformation is performed prior to averaging
1101        if transform != None:
1102            transformValues(et, fn=transform) #in place transform
1103            cbc()
1104
1105        #if average function is given, use it to join same spotids
1106        if average != None:
1107            et = averageAttributes(et, fn=average)
1108            cbc()
1109
1110        cbc.end()
1111
1112        return et
1113
1114    def get_single_data(self, *args, **kwargs):
1115        return self.get_data(*args, **kwargs)
1116
1117class DatabaseConnection(DictyExpress):
1118    pass
1119
1120def allAnnotationVals(annots):
1121    """
1122    All annotation valuess for given annotations
1123    in a dict of { name: set of possible values } pairs.
1124    """
1125    av = defaultdict(set)
1126    for a in annots:
1127        for name,val in a:
1129    return av
1130
1131def createExampleTable(names, vals, annots, ddb, cname="DDB", \
1132        exclude_constant_labels=False, permutation=None, always_include=["id"], allowed_labels=None):
1133    """
1134    Create an ExampleTable for this group. Attributes are those in
1135    names.
1136    """
1137    attributes = [ orange.FloatVariable(n, numberOfDecimals=3) \
1138        for n in names ]
1139
1140    #exclusion of names with constant values
1141    annotsvals = allAnnotationVals(annots)
1142    oknames = set(annotsvals.keys())
1143    if exclude_constant_labels:
1144        oknames = set(nth(filter(lambda x: len(x[1]) > 1 or x[0] in always_include,
1145            annotsvals.items()), 0))
1146
1147    if allowed_labels != None:
1148        oknames = set(filter(lambda x: x in allowed_labels, oknames))
1149
1150    #print oknames
1151
1152    for a,an in zip(attributes, annots):
1153        a.setattr("attributes", dict([(name,str(val)) for name,val in an if name in oknames]))
1154
1155    domain = orange.Domain(attributes, False)
1156    ddbv = orange.StringVariable(cname)
1157    id = orange.newmetaid()
1159
1160    examples = [ None ]*len(ddb)
1161    for i,(v,d) in enumerate(izip(izip(*vals), ddb)):
1162        ex = orange.Example(domain, [ floatOrUnknown(a) for a in v ])
1163        ex[cname] = d
1164
1165        if permutation:
1166            examples[permutation[i]] = ex
1167        else:
1168            examples[i] = ex
1169
1170    return orange.ExampleTable(domain,examples)
1171
1172def transformValues(data, fn):
1173    """
1174    In place transformation.
1175    """
1176    for ex in data:
1177        for at in data.domain.attributes:
1178            if not ex[at].isSpecial():
1179                ex[at] = fn(ex[at])
1180
1181def averageAttributes(data, joinc="DDB", fn=median):
1182    """
1183    Averages attributes with the same "join" parameter using
1184    specified function. Builds a now dataset. Order of
1185    "join" parameter stays the same with regard to first
1186    appearance.
1187    """
1188
1189    if verbose:
1190        print "Averaging attributes"
1191
1192    valueso = []
1193    valuess = set(valueso)
1194
1195    attributes = [ a for a in data.domain.attributes ]
1196    domain = data.domain
1197
1198    #accumulate values
1199
1200    for ex in data:
1201        j = str(ex[joinc])
1202
1203        if j not in valuess:
1204            valueso.append(j)
1206
1207    valuesddb = [ dict( (n,[]) for n in valueso) for at in attributes ]
1208
1209    for ex in data:
1210        j = str(ex[joinc])
1211
1212        for a in range(len(attributes)):
1213            val = ex[a]
1214            if not val.isSpecial():
1215                valuesddb[a][j].append(val.native())
1216
1217    #create a new example table reusing the domain - apply fn
1218    examples = []
1219    for v in valueso:
1220        example = orange.Example(domain, \
1221            [ floatOrUnknown(fn(valuesddb[a][v])) for a in range(len(attributes)) ] )
1222        example[joinc] = v
1223        examples.append(example)
1224
1225    return orange.ExampleTable(domain, examples)
1226
1227def floatOrUnknown(a):
1228    """
1229    Converts an element to float if possible.
1230    If now, output "?".
1231    """
1232    try:
1233        return float(a)
1234    except:
1235        return "?"
1236
1237class CallBack():
1238    """
1239    Converts "allparts" callbacks into by "callbacks"
1240    specified number of callbacks of function fn.
1241    """
1242
1243    def __init__(self, allparts, fn, callbacks=100):
1244        self.allparts = allparts
1245        self.lastreport = 0.00001
1246        self.getparts = 0
1247        self.increase = 1.0/callbacks
1248        self.callbacks = callbacks
1249        self.fn = fn
1250        self.cbs = 0
1251
1252    def __call__(self):
1253        self.getparts += 1
1254        done = float(self.getparts)/self.allparts
1255        while done > self.lastreport + self.increase:
1256            self.lastreport += self.increase
1257            self.fn()
1258            self.cbs += 1
1259
1260    def end(self):
1261        while self.cbs < self.callbacks:
1262            self.fn()
1263            self.cbs += 1
1264
1265class BufferSQLite(object):
1266
1267    def __init__(self, filename, compress=True):
1268        self.compress = compress
1269        self.filename = filename
1270        self.conn = self.connect()
1271
1272    def clear(self):
1273        """
1274        Removes all entries in the buffer
1275        """
1276        self.conn.close()
1277        os.remove(self.filename)
1278        self.conn = self.connect()
1279
1280    def connect(self):
1281        import sqlite3
1282        conn = sqlite3.connect(self.filename)
1283        c = conn.cursor()
1284        c.execute('''create table if not exists buf
1285        (address text primary key, time text, con blob)''')
1286        c.close()
1287        conn.commit()
1288        return conn
1289
1291        """ Returns version or False, if it does not exists """
1292        c = self.conn.cursor()
1294        lc = list(c)
1295        c.close()
1296        if len(lc) == 0:
1297            return False
1298        else:
1299            return lc[0][0]
1300
1301    def list(self):
1302        c = self.conn.cursor()
1304        return nth(list(c), 0)
1305
1307        import cPickle, zlib, sqlite3
1308        if verbose:
1310        c = self.conn.cursor()
1311        if self.compress:
1312            bin = sqlite3.Binary(zlib.compress(cPickle.dumps(con)))
1313        else:
1314            bin = sqlite3.Binary(cPickle.dumps(con))
1315        c.execute('insert or replace into buf values (?,?,?)', (addr, version, bin))
1316        c.close()
1317        if autocommit:
1318            self.commit()
1319
1320    def commit(self):
1321        self.conn.commit()
1322
1324        import cPickle, zlib
1325        if verbose:
1326            print "getting from buffer", addr
1327            t = time.time()
1328        c = self.conn.cursor()
1330        ls = list(c)
1331        first = ls[0][0]
1332        if verbose:
1333            print time.time() - t
1334        if self.compress:
1336        else:
1338        c.close()
1339
1340        if verbose:
1341            print time.time() - t
1342        return rc
1343
1345    def do():
1346        return urllib2.urlopen(url)
1347
1348    if repeat <= 0:
1349        do()
1350    else:
1351        try:
1352            return do()
1353        except:
1355
1356def empty_none(s):
1357    if s:
1358        return s
1359    else:
1360        return None
1361
1362def join_ats(atts):
1363    """ Joins attribute attributes together. If all values are the same,
1364    set the parameter to the common value, else return a list of the
1365    values in the same order as the attributes are imputed. """
1366    keys = reduce(lambda x,y: x | y, (set(at.keys()) for at in atts))
1367    od = {}
1368
1369    def man(x):
1370        if issequencens(x):
1371            return tuple(x)
1372        else:
1373            return x
1374
1375    for k in keys:
1376        s = set(man(at[k]) for at in atts)
1377        if len(s) == 1:
1378            od[k] = list(s)[0]
1379        else:
1380            od[k] = str([ at[k] for at in atts ])
1381    return od
1382
1383def join_replicates(data, ignorenames=["replicate", "id", "name", "map_stop1"], namefn=None, avg=median):
1384    """ Join replicates by median.
1385    Default parameters work for PIPA data.
1386    Sort elements in the same order as ignorenames!
1387    """
1388
1389    d = defaultdict(list)
1390
1391    if namefn == None:
1392        namefn = lambda att: ",".join(att["id"]) if issequencens(att["id"]) else att["id"]
1393
1394    #key function
1395    def key_g(att):
1396        dk = att.copy()
1397        for iname in ignorenames:
1398            dk.pop(iname, None)
1399
1400        def man(x):
1401            if issequencens(x):
1402                return tuple(x)
1403            else:
1404                return x
1405
1406        return tuple(nth(sorted(((a, man(b)) for a,b in dk.items())), 1))
1407
1408    #prepare groups
1409    for i,a in enumerate(data.domain.attributes):
1410        att = a.attributes
1411        k = key_g(att)
1412        d[k].append(i)
1413
1414    d = dict(d) #want errors with wrong keys
1415
1416    natts = []
1417
1418    def nativeOrNone(val):
1419        if val.isSpecial():
1420            return None
1421        else:
1422            return val.native()
1423
1424    def avgnone(l):
1425        """ Removes None and run avg function"""
1426        l = filter(lambda x: x != None, l)
1427        if len(l):
1428            return avg(l)
1429        else:
1430            return None
1431
1432    def data_type(vals): #data type converter
1433        try:
1434            _ = [ int(a) for a in vals ]
1435            return int
1436        except:
1437            try:
1438                _ = [ float(a) for a in vals ]
1439                return float
1440            except:
1441                return lambda x: x
1442
1443    all_values = defaultdict(set)
1444    for a in [ at.attributes for at in data.domain.attributes ]:
1445        for k,v in a.iteritems():
1447
1448    types = {}
1449    for k,vals in all_values.iteritems():
1450        types[k] = data_type(vals)
1451
1452    for group, elements in d.items():
1453        a = orange.FloatVariable()
1454        #here sort elements
1455
1456        def sk(x):
1457            return [ types[n](x[n]) for n in ignorenames if n in all_values ]
1458
1459        elements = sorted(elements, key=lambda x: sk(data.domain.attributes[x].attributes))
1460
1461        a.attributes.update(join_ats([data.domain.attributes[i].attributes for i in elements]))
1462        a.name = namefn(a.attributes)
1463
1464        def avgel(ex, el):
1465            return orange.Value(avgnone([ nativeOrNone(ex[ind]) for ind in el ]))
1466
1467        a.getValueFrom = lambda ex,rw,el=elements: avgel(ex,el)
1468        natts.append(a)
1469
1470    ndom = orange.Domain(natts, data.domain.classVar)
1472    return orange.ExampleTable(ndom, data)
1473
1474
1475class DictyBase(object):
1476
1477    domain = "dictybase"
1478    filename = "information_mappings.pck"
1479    tags = [ "Dictyostelium discoideum", "gene", "essential", "dictyBase" ]
1480
1481    @classmethod
1482    def version(cls):
1484        return orngServerFiles.info(cls.domain, cls.filename)["datetime"]
1485
1486    @classmethod
1488        """
1490        Returns a dictionary {ID: (name, synonyms, products)}
1491        """
1493        out = []
1494        for l in txt2ll(s, separ='\t', lineSepar='\n')[1:]:
1495            if len(l) == 4:
1496                id = l[0]
1497                name = l[1]
1498                synonyms = filter(None, l[2].split(", "))
1499                products = l[3]
1500                out.append((id, name, synonyms, products))
1501        return dict((a,(b,c,d)) for a,b,c,d in out)
1502
1503    @classmethod
1505        """
1507        Returns a list of (ddb, ddb_g, uniprot) triplets.
1508
1509        2009/04/07: ddb's appear unique
1510        """
1512        out = []
1513        for l in txt2ll(s, separ='\t', lineSepar='\n')[1:]:
1514            if len(l) == 4:
1515                ddb = empty_none(l[0])
1516                ddb_g = empty_none(l[1])
1517                name = empty_none(l[2])
1518                uniprot = empty_none(l[3])
1519                out.append((ddb, ddb_g, name, uniprot))
1520        return out
1521
1522    @classmethod
1523    def pickle_data(cls):
1526        return pickle.dumps((info,mappings), -1)
1527
1528    def __init__(self):
1530        self.info, self.mappings = pickle.load(open(fn, 'rb'))
1531
1532if __name__=="__main__":
1533    verbose = 1
1534
1535    def printet(et):
1536        et.save("ett.tab")
1538
1539    """
1540    a = DictyBase()
1541    print len(a.info)
1542
1543    dbc = DictyExpress(buffer=BufferSQLite("../tmpbufnew"))
1544
1545    print dbc.annotationOptions()
1546
1547    count = 0
1548    def cb():
1549        global count
1550        count += 1
1551        #print "CBBB", count
1552
1553    et = dbc.get_single_data(sample=[ "tagA-", "pkaC-"], callback=cb, exclude_constant_labels=True, allowed_labels=["sample"])
1554    print et.domain
1555    print et.domain[0].attributes
1556    printet(et)
1557    """
1558
1559    d = PIPA(buffer=BufferSQLite("../tmpbufnewpipa"))
1560
1561    print d.gene_expression_types()
1562
1563    cov = d.coverage("777", chromosome=1, genome="dd")
1564    print cov
1565
1566    allids = d.annotations().keys()
1567    allids = allids[:1]
1568    print ("list", allids)
1569    print d.annotations().items()[0]
1570    print ("annots", d.annotations().items()[:2])
1571
1572    allids = map(str, allids)
1573
1574    import math
1575
1576    data = d.get_data(ids=allids, ctype="3")
1577
1578    #data = d.get_data(ids=allids, transform=lambda x: math.log(x+1, 2), allowed_labels=["strain"], ctype="3")
1579
1580    printet(data)
Note: See TracBrowser for help on using the repository browser.