source: orange-bioinformatics/obiAssess.py @ 1592:c39f563cc473

Revision 1592:c39f563cc473, 21.6 KB checked in by markotoplak, 2 years ago (diff)

Added Mean and Median to the new gene set signatures module.

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
485    return a
486
487def pca(data, snapshot=0):
488    "Perform PCA on M, return eigenvectors and eigenvalues, sorted."
489    M = data.toNumpy("a")[0]
490    XMean = numpy.mean(M, axis = 0)
491    print XMean.shape, M.shape
492    M = M - XMean
493
494    T, N = numpy.shape(M)
495    # if there are less rows T than columns N, use snapshot method
496    if (T < N) or snapshot:
497        C = numpy.dot(M, numpy.transpose(M))
498        evals, evecsC = numpy.linalg.eigh(C) #columns of evecsC are eigenvectors
499        evecs = numpy.dot(M.T, evecsC)/numpy.sqrt(numpy.abs(evals))
500    else:
501        K = numpy.dot(numpy.transpose(M), M)
502        evals, evecs = numpy.linalg.eigh(K)
503   
504    evecs = numpy.transpose(evecs)
505
506    # sort the eigenvalues and eigenvectors, decending order
507    order = (numpy.argsort(numpy.abs(evals))[::-1])
508    evecs = numpy.take(evecs, order, 0)
509    evals = numpy.take(evals, order)
510    return evals, evecs, XMean
511
512class PCALearner(object):
513    """ Transforms gene sets using Principal Leasts Squares. """
514   
515    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None):
516
517        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
518            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
519   
520        constructt = {}
521
522        #build weight matrices for every gene set
523        for name, inds in gsetsnum.items():
524            dom2 = orange.Domain([ data.domain.attributes[i1] for i1 in inds ], data.domain.classVar)
525
526            data_gs = orange.ExampleTable(dom2, data)
527            evals, evect, xmean = pca(data_gs)
528            constructt[name] = inds, evals, evect, xmean
529
530        return PCA(constructt=constructt)
531
532
533class SimpleFun(object):
534
535    def __init__(self, **kwargs):
536        for a,b in kwargs.items():
537            setattr(self, a, b)
538
539    def __call__(self, example):
540        #for  name,ids in self.gsets.items():
541        #    print name, [example[i].value for i in ids], self.fn([example[i].value for i in ids])
542        return dict( (name, self.fn([example[i].value for i in ids])) \
543            for name,ids in self.gsets.items() )
544
545class SimpleFunLearner(object):
546    """
547    Just applies a function taking attribute values of an example and
548    produces to all gene sets.   
549    """
550    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, fn=None):
551        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
552            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
553        return SimpleFun(gsets=gsetsnum, fn=fn)
554
555class MedianLearner(object):
556   
557    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, fn=None):
558       sfl =  SimpleFunLearner()
559       return sfl(data, matcher, geneSets, \
560            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues, fn=statc.median)
561
562class MeanLearner(object):
563   
564    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None, fn=None):
565       sfl =  SimpleFunLearner()
566       return sfl(data, matcher, geneSets, \
567            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues, fn=statc.mean)
568
569def impute_missing(data):
570    #remove attributes with only unknown values
571    newatts = []
572    for at in data.domain.attributes:
573        svalues = [ 1 for a in data if a[at].isSpecial() ]
574        real = len(data) - len(svalues)
575        if real > 0:
576            newatts.append(at)
577
578    dom2 = orange.Domain(newatts, data.domain.classVar)
579    data = orange.ExampleTable(dom2, data)
580
581    #impute
582    import orngTree 
583    imputer = orange.ImputerConstructor_model() 
584    imputer.learnerContinuous = imputer.learnerDiscrete = orange.MajorityLearner()
585    imputer = imputer(data)
586
587    data = imputer(data)
588    return data
589
590def setSig_example(ldata, ex, genesets):
591    """
592    Create an dictionary with (geneset_name, geneset_expression) for every
593    geneset on input.
594
595    ldata is class-annotated data
596    """
597    #seznam ocen genesetov za posamezen primer v ucni mzozici
598    pathways = {}
599
600    def setSig_example_geneset(ex, data):
601        """ ex contains only selected genes """
602
603        distances = [ [], [] ]   
604
605        def pearsonr(v1, v2):
606            try:
607                return statc.pearsonr(v1, v2)[0]
608            except:
609                return numpy.corrcoef([v1, v2])[0,1]
610
611        def pearson(ex1, ex2):
612            attrs = range(len(ex1.domain.attributes))
613            vals1 = [ ex1[i].value for i in attrs ]
614            vals2 = [ ex2[i].value for i in attrs ]
615            return pearsonr(vals1, vals2)
616
617        def ttest(ex1, ex2):
618            try:
619                return stats.lttest_ind(ex1, ex2)[0]
620            except:
621                return 0.0
622       
623        #maps class value to its index
624        classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.classVar.values) ])
625     
626        #create distances to all learning data - save or other class
627        for c in data:
628            distances[classValueMap[c[-1].value]].append(pearson(c, ex))
629
630        return ttest(distances[0], distances[1])
631           
632    for name, gs in genesets.items(): #for each geneset
633        #for each gene set: take the attribute subset and work on the attribute subset only
634        #only select the subset of genes from the learning data
635        domain = orange.Domain([ldata.domain.attributes[ai] for ai in gs], ldata.domain.classVar)
636        datao = orange.ExampleTable(domain, ldata)
637        example = orange.Example(domain, ex) #domains need to be the same
638     
639        ocena = setSig_example_geneset(example, datao)
640        pathways[name] = ocena
641       
642    return pathways
643
644class SetSig(object):
645
646    def __init__(self, **kwargs):
647        for a,b in kwargs.items():
648            setattr(self, a, b)
649
650    def __call__(self, example):
651        return setSig_example(self.learndata, example, self.genesets)
652
653class SetSigLearner(object):
654
655    def __call__(self, data, matcher, geneSets, minSize=3, maxSize=1000, minPart=0.1, classValues=None):
656        data, oknames, gsetsnum = selectGenesetsData(data, matcher, geneSets, \
657            minSize=minSize, maxSize=maxSize, minPart=minPart, classValues=classValues)
658        return SetSig(learndata=data, genesets=gsetsnum)
659
660if __name__ == "__main__":
661
662    data = Orange.data.Table("iris")
663    gsets = obiGeneSets.collections({
664        "ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
665        "f3": ['sepal length', 'sepal width', 'petal length'],
666        "l3": ['sepal width', 'petal length', 'petal width'],
667        })
668
669    fp = 120
670    ldata = orange.ExampleTable(data.domain, data[:fp])
671    tdata = orange.ExampleTable(data.domain, data[fp:])
672
673    matcher = obiGene.matcher([])
674
675    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
676    #ass = AssessLearner()(data, matcher, gsets, rankingf=AT_loessLearner())
677    #ass = MeanLearner()(data, matcher, gsets, default=False)
678    ass = MedianLearner()(data, matcher, gsets)
679    #ass = PLSLearner()(data, matcher, gsets, classValues=choosen_cv)
680    #ass = SetSigLearner()(ldata, matcher, gsets, classValues=choosen_cv, minPart=0.0)
681    #ass = PCALearner()(ldata, matcher, gsets, classValues=choosen_cv, minPart=0.0)
682    #ass = GSALearner()(ldata, matcher, gsets, classValues=choosen_cv, minPart=0.0)
683
684    ar = defaultdict(list)
685    for d in (list(ldata) + list(tdata))[:5]:
686        for a,b in ass(d).items():
687            ar[a].append(b)
688
689    def pp1(ar):
690        ol =  sorted(ar.items())
691        print '\n'.join([ a.id + ": " +str(b) for a,b in ol])
692
693    pp1(ar)
694
Note: See TracBrowser for help on using the repository browser.