source: orange-bioinformatics/obiAssess.py @ 1597:574fe4f8d1d3

Revision 1597:574fe4f8d1d3, 21.7 KB checked in by markotoplak, 2 years ago (diff)

Added GSA to obiGeneSetSig.

Line 
1"""
2Construction of gene set scores for each sample.
3
4Learners in this module build models needed for construction
5of features from individual genes.
6
7The other classes just take example and return a
8dictionary of { name: score } for that example.
9"""
10
11import obiGsea
12import obiGeneSets
13import orange
14import Orange
15import stats
16import statc
17import numpy
18import math
19import obiExpression
20import obiGene
21from collections import defaultdict
22
23def normcdf(x, mi, st):
24    return 0.5*(2. - stats.erfcc((x - mi)/(st*math.sqrt(2))))
25
26class AT_edelmanParametric(object):
27
28    def __init__(self, **kwargs):
29        for a,b in kwargs.items():
30            setattr(self, a, b)
31
32    def __call__(self, nval):
33
34        if self.mi1 == None or self.mi2 == None or self.st1 == None or self.st2 == None:
35            return 0 
36
37        try:
38            val = nval.value
39            if nval.isSpecial():
40                return 0 
41        except:
42            val = nval
43
44        try:
45            if val >= self.mi1:
46                p1 = 1 - normcdf(val, self.mi1, self.st1)
47            else:
48                p1 = normcdf(val, self.mi1, self.st1)
49
50            if val >= self.mi2:
51                p2 = 1 - normcdf(val, self.mi2, self.st2)
52            else:
53                p2 = normcdf(val, self.mi2, self.st2)
54
55            #print p1, p2
56            return math.log(p1/p2)
57        except:
58            #print p1, p2, "exception"
59            return 0
60
61class AT_edelmanParametricLearner(object):
62    """
63    Returns attribute transfromer for Edelman parametric measure for a given attribute in the
64    dataset.
65    Edelman et al, 06.
66
67    Modified a bit.
68    """
69
70    def __init__(self, a=None, b=None):
71        """
72        a and b are choosen class values.
73        """
74        self.a = a
75        self.b = b
76
77    def __call__(self, i, data):
78        cv = data.domain.classVar
79        #print data.domain
80
81        if self.a == None: self.a = cv.values[0]
82        if self.b == None: self.b = cv.values[1]
83
84        def avWCVal(value):
85            return [ex[i].value for ex in data if ex[-1].value == value and not ex[i].isSpecial() ]
86
87        list1 = avWCVal(self.a)
88        list2 = avWCVal(self.b)
89
90        mi1 = mi2 = st1 = st2 = None
91
92        try:
93            mi1 = statc.mean(list1)
94            st1 = statc.std(list1)
95        except:
96            pass
97   
98        try:
99            mi2 = statc.mean(list2)
100            st2 = statc.std(list2)
101        except:
102            pass
103
104        return AT_edelmanParametric(mi1=mi1, mi2=mi2, st1=st1, st2=st2)
105
106class AT_loess(object):
107
108    def __init__(self, **kwargs):
109        for a,b in kwargs.items():
110            setattr(self, a, b)
111
112    def __call__(self, nval):
113
114        val = nval.value
115        if nval.isSpecial():
116            return 0.0 #middle value
117        #return first class probablity
118
119        import math
120
121        def saveplog(a,b):
122            try:
123                return math.log(a/b)
124            except:
125                if a < b:
126                    return -10
127                else:
128                    return +10
129
130        try:
131            ocene = self.condprob(val)
132            if sum(ocene) < 0.01:
133                return 0.0
134            return saveplog(ocene[0], ocene[1])
135
136        except:
137            return 0.0
138
139class AT_loessLearner(object):
140
141    def __call__(self, i, data):
142        cv = data.domain.classVar
143        #print data.domain
144        try:
145            ca = orange.ContingencyAttrClass(data.domain.attributes[i], data)
146            a = orange.ConditionalProbabilityEstimatorConstructor_loess(ca, nPoints=5)
147            return AT_loess(condprob=a)
148        except:
149            return AT_loess(condprob=None)
150
151def nth(l, n):
152    return [a[n] for a in l]
153
154class Assess(object):
155
156    def __init__(self, **kwargs):
157        for a,b in kwargs.items():
158            setattr(self, a, b)
159
160    def __call__(self, example):
161        enrichmentScores = [] 
162
163        lcor = [ self.attrans[at](example[at]) for at in range(len(self.attrans)) ]
164
165        ordered = obiGsea.orderedPointersCorr(lcor)
166
167        def rev(l):
168           return numpy.argsort(l)
169
170        rev2 = rev(ordered)
171
172        gsetsnumit = self.gsetsnum.items()
173
174        enrichmentScores = dict( 
175            (name, obiGsea.enrichmentScoreRanked(subset, lcor, ordered, rev2=rev2)[0]) \
176            for name,subset in gsetsnumit)
177   
178        return enrichmentScores
179
180
181class AssessLearner(object):
182    """
183    Uses the underlying GSEA code to select genes.
184    Takes data and creates attribute transformations.
185    """
186   
187    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, rankingf=None):
188        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
189            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
190       
191        if rankingf == None:
192            rankingf = AT_edelmanParametricLearner()
193
194        #rank individual attributes on the training set
195        attrans = [ rankingf(iat, data) for iat, at in enumerate(data.domain.attributes) ]
196
197        return Assess(attrans=attrans, gsetsnum=gsetsnum)
198
199def selectGenesetsData(data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None):
200    """
201    Returns gene sets and data which falling under upper criteria.
202    """
203    gso = obiGsea.GSEA(data, matcher=matcher, classValues=classValues, atLeast=0)
204    gso.addGenesets(geneSets)
205    okgenesets = gso.selectGenesets(minSize=minSize, maxSize=maxSize, minPart=minPart).keys()
206    gsetsnum = gso.to_gsetsnum(okgenesets)
207    return gso.data, okgenesets, gsetsnum
208
209def corgs_activity_score(ex, corg):
210    """ activity score for a sample for pathway given by corgs """
211    #print [ ex[i].value for i in corg ] #FIXME what to do with unknown values?
212    return sum(ex[i].value if ex[i].value != '?' else 0.0 for i in corg)/len(corg)**0.5
213
214class CORGs(object):
215
216    def __init__(self, **kwargs):
217        for a,b in kwargs.items():
218            setattr(self, a, b)
219
220    def __call__(self, example):
221        return dict( (name,corgs_activity_score(example, corg)) \
222            for name, corg in self.corgs.items() )
223
224class CORGsLearner(object):
225   
226    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None):
227        """
228        WARNING: input has to be z_ij table! each gene needs to be normalized
229        (mean=0, stdev=1) for all samples.
230        """
231        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
232            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
233   
234        tscorecache = {}
235
236        def tscorec(data, at, cache=None):
237            """ Cached attribute  tscore calculation """
238            if cache != None and at in cache: return cache[at]
239            ma = obiExpression.MA_t_test()(at,data)
240            if cache != None: cache[at] = ma
241            return ma
242
243        def compute_corg(data, inds):
244            """
245            Compute CORG for this geneset specified with gene inds
246            in the example table. Output is the list of gene inds
247            in CORG.
248
249            """
250            #order member genes by their t-scores: decreasing, if av(t-score) >= 0,
251            #else increasing
252            tscores = [ tscorec(data, at, tscorecache) for at in inds ]
253            sortedinds = nth(sorted(zip(inds,tscores), key=lambda x: x[1], \
254                reverse=statc.mean(tscores) >= 0), 0)
255
256            def S(corg):
257                """ Activity score separation - S(G) in
258                the article """
259                asv = orange.FloatVariable(name='AS')
260                asv.getValueFrom = lambda ex,rw: orange.Value(asv, corgs_activity_score(ex, corg))
261                data2 = orange.ExampleTable(orange.Domain([asv], data.domain.classVar), data)
262                return abs(tscorec(data2, 0)) #FIXME absolute - nothing in the article about it
263                   
264            #greedily find CORGS procing the best separation
265            g = S(sortedinds[:1])
266            bg = 1
267            for a in range(2, len(sortedinds)+1):
268                tg = S(sortedinds[:a])
269                if tg > g:
270                    g = tg
271                    bg = a
272                else:
273                    break
274               
275            return sortedinds[:a]
276
277        #now, on the learning set produce the CORGS for each geneset and
278        #save it for use in further prediction
279
280        corgs = {}
281
282        for name, inds in gsetsnum.items():
283            inds = sorted(set(inds)) # take each gene only once!
284            corgs[name] = compute_corg(data, inds)
285
286        return CORGs(corgs=corgs)
287
288class GSA(object):
289
290    def __init__(self, **kwargs):
291        for a,b in kwargs.items():
292            setattr(self, a, b)
293
294    def __call__(self, example):
295        return dict( (name, statc.mean([example[i].value for i in inds if example[i].value != "?"]) ) \
296            for name, inds in self.subsets.items() )
297
298class GSALearner(object):
299   
300    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None):
301        """
302        WARNING: input has to be z_ij table! each gene needs to be normalized
303        (mean=0, stdev=1) for all samples.
304        """
305        import scipy.stats
306
307        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
308            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
309   
310        def tscorec(data, at, cache=None):
311            ma = obiExpression.MA_t_test()(at,data)
312            return ma
313
314        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
315
316        def to_z_score(t):
317            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
318
319        zscores = map(to_z_score, tscores)
320
321        subsets = {}
322
323        for name, inds in gsetsnum.items():
324            inds = sorted(set(inds)) # take each gene only once!
325
326            D = statc.mean([max(zscores[i],0) for i in inds]) \
327                + statc.mean([min(zscores[i],0) for i in inds])
328
329            if D >= 0:
330                subsets[name] = [ i for i in inds if zscores[i] > 0.0 ]
331            else:
332                subsets[name] = [ i for i in inds if zscores[i] < 0.0 ]
333
334        return GSA(subsets=subsets)
335
336def pls_transform(example, constt):
337    """
338    Uses calculated PLS weights to transform the given example.
339    """
340
341    inds, xmean, W, P = constt
342    dom = orange.Domain([example.domain.attributes[i1] for i1 in inds ], False)
343    newex = orange.ExampleTable(dom, [example])
344    ex = newex.toNumpy()[0]
345    ex = ex - xmean # same input transformation
346
347    nc = W.shape[1]
348
349    TR = numpy.empty((1, nc))
350    XR = ex
351
352    dot = numpy.dot
353
354    for i in range(nc):
355       t = dot(XR, W[:,i].T)
356       XR = XR - t*numpy.array([P[:,i]])
357       TR[:,i] = t
358
359    return TR
360
361def PLSCall(data, y=None, x=None, nc=None, weight=None, save_partial=False):
362
363    def normalize(vector):
364        return vector / numpy.linalg.norm(vector)
365
366    if y == None:
367        y = [ data.domain.classVar ]
368    if x == None:
369        x = [v for v in data.domain.variables if v not in y]
370
371    Ncomp = nc if nc is not None else len(x)
372       
373    dataX = orange.ExampleTable(orange.Domain(x, False), data)
374    dataY = orange.ExampleTable(orange.Domain(y, False), data)
375
376    # transformation to numpy arrays
377    X = dataX.toNumpy()[0]
378    Y = dataY.toNumpy()[0]
379
380    # data dimensions
381    n, mx = numpy.shape(X)
382    my = numpy.shape(Y)[1]
383
384    # Z-scores of original matrices
385    YMean = numpy.mean(Y, axis = 0)
386    XMean = numpy.mean(X, axis = 0)
387   
388    X = (X-XMean)
389    Y = (Y-YMean)
390
391    P = numpy.empty((mx,Ncomp))
392    T = numpy.empty((n,Ncomp))
393    W = numpy.empty((mx,Ncomp))
394    E,F = X,Y
395
396    dot = numpy.dot
397    norm = numpy.linalg.norm
398
399    #PLS1 - from Gutkin, shamir, Dror: SlimPLS
400
401    for i in range(Ncomp):
402        w = dot(E.T,F)
403        w = w/norm(w) #normalize w in Gutkin et al the do w*c, where c is 1/norm(w)
404        t = dot(E, w) #t_i -> a row vector
405        p = dot(E.T, t)/dot(t.T, t) #p_i t.T is a row vector - this is inner(t.T, t.T)
406        q = dot(F.T, t)/dot(t.T, t) #q_i
407           
408        E = E - dot(t, p.T)
409        F = F - dot(t, q.T)
410
411        T[:,i] = t.T
412        W[:,i] = w.T
413        P[:,i] = p.T
414
415    return XMean, W, P, T
416
417class PLS(object):
418
419    def __init__(self, **kwargs):
420        for a,b in kwargs.items():
421            setattr(self, a, b)
422
423    def __call__(self, example):
424
425        od = {}
426
427        for name, constt in self.constructt.items():
428            ts = pls_transform(example, constt)[0]
429            if len(ts) == 1:
430                od[name] = ts[0]
431            else:
432                for i,s in enumerate(ts):
433                   od[name + "_LC_" + str(i+1)] = s
434 
435        return od
436
437class PLSLearner(object):
438    """ Transforms gene sets using Principal Leasts Squares. """
439   
440    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, components=1):
441        """
442        If more that 1 components are used, _LC_componetsNumber is appended to
443        the name of the gene set.
444        """
445
446        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
447            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
448   
449        constructt = {}
450
451        #build weight matrices for every gene set
452        for name, inds in gsetsnum.items():
453            dom2 = orange.Domain([ data.domain.attributes[i1] for i1 in inds ], data.domain.classVar)
454            data_gs = orange.ExampleTable(dom2, data)
455            xmean, W, P, T = PLSCall(data_gs, nc=components, y=[data_gs.domain.classVar], save_partial=True)
456            constructt[name] = inds, xmean, W, P
457
458        return PLS(constructt=constructt)
459
460class PCA(object):
461
462    def __init__(self, **kwargs):
463        for a,b in kwargs.items():
464            setattr(self, a, b)
465
466    def __call__(self, example):
467        od = {}
468
469        for name, constt in self.constructt.items():
470            ts = pca_transform(example, constt)[0]
471            od[name] = ts
472
473        return od
474
475def pca_transform(example, constt):
476    inds, evals, evect, xmean = constt
477    dom = orange.Domain([example.domain.attributes[i1] for i1 in inds ], False)
478    newex = orange.ExampleTable(dom, [example])
479    arr = newex.toNumpy()[0]
480    arr = arr - xmean # same input transformation - a row in a matrix
481
482    ev0 = evect[0] #this is a row in a matrix - do a dot product
483    a = numpy.dot(arr, ev0)
484    return a
485
486def pca(data, snapshot=0):
487    "Perform PCA on M, return eigenvectors and eigenvalues, sorted."
488    M = data.toNumpy("a")[0]
489    XMean = numpy.mean(M, axis = 0)
490    M = M - XMean
491
492    T, N = numpy.shape(M)
493    # if there are less rows T than columns N, use snapshot method
494    if (T < N) or snapshot:
495        C = numpy.dot(M, numpy.transpose(M))
496        evals, evecsC = numpy.linalg.eigh(C) #columns of evecsC are eigenvectors
497        evecs = numpy.dot(M.T, evecsC)/numpy.sqrt(numpy.abs(evals))
498    else:
499        K = numpy.dot(numpy.transpose(M), M)
500        evals, evecs = numpy.linalg.eigh(K)
501   
502    evecs = numpy.transpose(evecs)
503
504    # sort the eigenvalues and eigenvectors, decending order
505    order = (numpy.argsort(numpy.abs(evals))[::-1])
506    evecs = numpy.take(evecs, order, 0)
507    evals = numpy.take(evals, order)
508    return evals, evecs, XMean
509
510class PCALearner(object):
511    """ Transforms gene sets using Principal Leasts Squares. """
512   
513    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None):
514
515        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
516            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
517   
518        constructt = {}
519
520        #build weight matrices for every gene set
521        for name, inds in gsetsnum.items():
522            dom2 = orange.Domain([ data.domain.attributes[i1] for i1 in inds ], data.domain.classVar)
523
524            data_gs = orange.ExampleTable(dom2, data)
525            evals, evect, xmean = pca(data_gs)
526            constructt[name] = inds, evals, evect, xmean
527
528        return PCA(constructt=constructt)
529
530
531class SimpleFun(object):
532
533    def __init__(self, **kwargs):
534        for a,b in kwargs.items():
535            setattr(self, a, b)
536
537    def __call__(self, example):
538        #for  name,ids in self.gsets.items():
539        #    print name, [example[i].value for i in ids], self.fn([example[i].value for i in ids])
540        return dict( (name, self.fn([example[i].value for i in ids])) \
541            for name,ids in self.gsets.items() )
542
543class SimpleFunLearner(object):
544    """
545    Just applies a function taking attribute values of an example and
546    produces to all gene sets.   
547    """
548    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, fn=None):
549        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
550            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
551        return SimpleFun(gsets=gsetsnum, fn=fn)
552
553class MedianLearner(object):
554   
555    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, fn=None):
556       sfl =  SimpleFunLearner()
557       return sfl(data, matcher, geneSets, \
558            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues, fn=statc.median)
559
560class MeanLearner(object):
561   
562    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, fn=None):
563       sfl =  SimpleFunLearner()
564       return sfl(data, matcher, geneSets, \
565            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues, fn=statc.mean)
566
567def impute_missing(data):
568    #remove attributes with only unknown values
569    newatts = []
570    for at in data.domain.attributes:
571        svalues = [ 1 for a in data if a[at].isSpecial() ]
572        real = len(data) - len(svalues)
573        if real > 0:
574            newatts.append(at)
575
576    dom2 = orange.Domain(newatts, data.domain.classVar)
577    data = orange.ExampleTable(dom2, data)
578
579    #impute
580    import orngTree 
581    imputer = orange.ImputerConstructor_model() 
582    imputer.learnerContinuous = imputer.learnerDiscrete = orange.MajorityLearner()
583    imputer = imputer(data)
584
585    data = imputer(data)
586    return data
587
588def setSig_example(ldata, ex, genesets):
589    """
590    Create an dictionary with (geneset_name, geneset_expression) for every
591    geneset on input.
592
593    ldata is class-annotated data
594    """
595    #seznam ocen genesetov za posamezen primer v ucni mzozici
596    pathways = {}
597
598    def setSig_example_geneset(ex, data):
599        """ ex contains only selected genes """
600
601        distances = [ [], [] ]   
602
603        def pearsonr(v1, v2):
604            try:
605                return statc.pearsonr(v1, v2)[0]
606            except:
607                return numpy.corrcoef([v1, v2])[0,1]
608
609        def pearson(ex1, ex2):
610            attrs = range(len(ex1.domain.attributes))
611            vals1 = [ ex1[i].value for i in attrs ]
612            vals2 = [ ex2[i].value for i in attrs ]
613            return pearsonr(vals1, vals2)
614
615        def ttest(ex1, ex2):
616            try:
617                return stats.lttest_ind(ex1, ex2)[0]
618            except:
619                return 0.0
620       
621        #maps class value to its index
622        classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.classVar.values) ])
623     
624        #create distances to all learning data - save or other class
625        for c in data:
626            distances[classValueMap[c[-1].value]].append(pearson(c, ex))
627
628        return ttest(distances[0], distances[1])
629           
630    for name, gs in genesets.items(): #for each geneset
631        #for each gene set: take the attribute subset and work on the attribute subset only
632        #only select the subset of genes from the learning data
633        domain = orange.Domain([ldata.domain.attributes[ai] for ai in gs], ldata.domain.classVar)
634        datao = orange.ExampleTable(domain, ldata)
635        example = orange.Example(domain, ex) #domains need to be the same
636     
637        ocena = setSig_example_geneset(example, datao)
638        pathways[name] = ocena
639       
640    return pathways
641
642class SetSig(object):
643
644    def __init__(self, **kwargs):
645        for a,b in kwargs.items():
646            setattr(self, a, b)
647
648    def __call__(self, example):
649        return setSig_example(self.learndata, example, self.genesets)
650
651class SetSigLearner(object):
652
653    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None):
654        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
655            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
656        return SetSig(learndata=data, genesets=gsetsnum)
657
658if __name__ == "__main__":
659
660    data = Orange.data.Table("iris")
661    gsets = obiGeneSets.collections({
662        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
663        "f3": ['sepal length', 'sepal width', 'petal length'],
664        "l3": ['sepal width', 'petal length', 'petal width'],
665        })
666
667    fp = 120
668    ldata = orange.ExampleTable(data.domain, data[:fp])
669    tdata = orange.ExampleTable(data.domain, data[fp:])
670
671    matcher = obiGene.matcher([])
672
673    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
674    #ass = AssessLearner()(data, matcher, gsets, rankingf=AT_loessLearner())
675    #ass = AssessLearner()(data, matcher, gsets, minPart=0.0)
676    #ass = MeanLearner()(data, matcher, gsets, default=False)
677    #ass = MedianLearner()(data, matcher, gsets)
678    #ass = PLSLearner()(data, matcher, gsets, classValues=choosen_cv, minPart=0.0)
679    #ass = SetSigLearner()(ldata, matcher, gsets, classValues=choosen_cv, minPart=0.0)
680    #ass = PCALearner()(ldata, matcher, gsets, classValues=choosen_cv, minPart=0.0)
681    ass = GSALearner()(ldata, matcher, gsets, classValues=choosen_cv, minPart=0.0)
682
683    ar = defaultdict(list)
684    for d in (list(ldata) + list(tdata))[:5]:
685        for a,b in ass(d).items():
686            ar[a].append(b)
687
688    def pp1(ar):
689        ol =  sorted(ar.items())
690        print '\n'.join([ a.id + ": " +str(b) for a,b in ol])
691
692    pp1(ar)
693
Note: See TracBrowser for help on using the repository browser.