source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1707:7d692228fea1

Revision 1707:7d692228fea1, 19.6 KB checked in by markotoplak, 20 months ago (diff)

Speed up of setsig (obiGeneSetSig).

Line 
1from __future__ import absolute_import
2
3import math
4from collections import defaultdict
5
6import scipy.stats
7
8import numpy
9
10import Orange, Orange.utils, statc
11
12if __name__ == "__main__":
13    __package__ = "Orange.bio"
14
15from .obiGsea import takeClasses
16from .obiAssess import pca, PLSCall, corgs_activity_score
17from . import obiExpression, obiGene, obiGeneSets, obiGsea, stats
18
19class GeneSetTrans(object):
20
21    __new__ = Orange.utils._orange__new__(object)
22
23    def _mat_ni(self, data):
24        """ With cached gene matchers. """
25        if data.domain not in self._cache:
26            self._cache[data.domain] = mat_ni(data, self.matcher)
27        return self._cache[data.domain]
28
29    def _match_instance(self, instance, geneset, takegenes=None):
30        nm, name_ind = self._mat_ni(instance)
31        genes = [ nm.umatch(gene) for gene in geneset ]
32        if takegenes:
33            genes = [ genes[i] for i in takegenes ]
34        return nm, name_ind, genes
35
36    def _match_data(self, data, geneset, odic=False):
37        nm, name_ind = self._mat_ni(data)
38        genes = [ nm.umatch(gene) for gene in geneset ]
39        if odic:
40            to_geneset = dict(zip(genes, geneset))
41        takegenes = [ i for i,a in enumerate(genes) if a != None ]
42        genes = [ genes[i] for i in takegenes ]
43        if odic:
44            return nm, name_ind, genes, takegenes, to_geneset
45        else:
46            return nm, name_ind, genes, takegenes
47
48    def __init__(self, matcher=None, gene_sets=None, min_size=3, max_size=1000, min_part=0.1, class_values=None, cv=False):
49        self.matcher = matcher
50        self.gene_sets = gene_sets
51        self.min_size = min_size
52        self.max_size = max_size
53        self.min_part = min_part
54        self.class_values = class_values
55        self._cache = {}
56        self.cv = cv
57
58    def __call__(self, data, weight_id=None):
59
60        #selection of classes and gene sets
61        data = takeClasses(data, classValues=self.class_values)
62        nm,_ =  self._mat_ni(data)
63        gene_sets = select_genesets(nm, self.gene_sets, self.min_size, self.max_size, self.min_part)
64
65        #build a new domain
66        newfeatures = self.build_features(data, gene_sets)
67        newdomain = Orange.data.Domain(newfeatures, data.domain.class_var)
68
69        #build a data set with cross validation
70        if self.cv == False:
71            return Orange.data.Table(newdomain, data)
72        else:
73            # The domain has the transformer that is build on all samples,
74            # while the transformed data table uses cross-validation
75            # internally
76            folds = 5
77            cvi = Orange.data.sample.SubsetIndicesCV(data, folds)
78            data_cv = [ [] for _ in range(len(data)) ]
79            for f in range(folds):
80                learn = data.select(cvi, f, negate=True)
81                test = data.select(cvi, f)
82                lf = self.build_features(learn, gene_sets)
83                transd = Orange.data.Domain(lf, data.domain.class_var)
84                trans_test = Orange.data.Table(transd, test)
85                for ex, pos in \
86                    zip(trans_test, [ i for i,n in enumerate(cvi) if n == f ]):
87                    data_cv[pos] = ex.native(0)
88            return Orange.data.Table(newdomain, data_cv)
89
90    def build_features(self, data, gene_sets):
91        return [ self.build_feature(data, gs) for gs in gene_sets ]
92
93def normcdf(x, mi, st):
94    return 0.5*(2. - stats.erfcc((x - mi)/(st*math.sqrt(2))))
95
96class AT_edelmanParametric(object):
97
98    def __init__(self, **kwargs):
99        for a,b in kwargs.items():
100            setattr(self, a, b)
101
102    def __call__(self, nval):
103
104        if self.mi1 == None or self.mi2 == None or self.st1 == None or self.st2 == None:
105            return 0.0 
106
107        val = nval
108
109        try:
110            if val >= self.mi1:
111                p1 = 1 - normcdf(val, self.mi1, self.st1)
112            else:
113                p1 = normcdf(val, self.mi1, self.st1)
114
115            if val >= self.mi2:
116                p2 = 1 - normcdf(val, self.mi2, self.st2)
117            else:
118                p2 = normcdf(val, self.mi2, self.st2)
119
120            #print p1, p2
121            return math.log(p1/p2)
122        except:
123            #print p1, p2, "exception"
124            return 0
125
126class AT_edelmanParametricLearner(object):
127    """
128    Returns attribute transfromer for Edelman parametric measure for a
129    given attribute in the dataset.  Edelman et al, 06. Modified a bit.
130    """
131
132    def __init__(self, a=None, b=None):
133        """
134        a and b are choosen class values.
135        """
136        self.a = a
137        self.b = b
138
139    def __call__(self, i, data):
140        cv = data.domain.classVar
141        #print data.domain
142
143        if self.a == None: self.a = cv.values[0]
144        if self.b == None: self.b = cv.values[1]
145
146        def avWCVal(value):
147            return [ex[i].value for ex in data if ex[-1].value == value and not ex[i].isSpecial() ]
148
149        list1 = avWCVal(self.a)
150        list2 = avWCVal(self.b)
151
152        mi1 = mi2 = st1 = st2 = None
153
154        try:
155            mi1 = statc.mean(list1)
156            st1 = statc.std(list1)
157        except:
158            pass
159   
160        try:
161            mi2 = statc.mean(list2)
162            st2 = statc.std(list2)
163        except:
164            pass
165
166        return AT_edelmanParametric(mi1=mi1, mi2=mi2, st1=st1, st2=st2)
167
168class AT_loess(object):
169
170    def __init__(self, **kwargs):
171        for a,b in kwargs.items():
172            setattr(self, a, b)
173
174    def __call__(self, nval):
175
176        val = nval
177
178        def saveplog(a,b):
179            try:
180                return math.log(a/b)
181            except:
182                if a < b:
183                    return -10
184                else:
185                    return +10
186
187        try:
188            ocene = self.condprob(val)
189            if sum(ocene) < 0.01:
190                return 0.0
191            return saveplog(ocene[0], ocene[1])
192
193        except:
194            return 0.0
195
196class AT_loessLearner(object):
197
198    def __call__(self, i, data):
199        try:
200            ca = Orange.statistics.contingency.VarClass(data.domain.attributes[i], data)
201            a =  Orange.statistics.estimate.ConditionalLoess(ca, nPoints=5)
202            return AT_loess(condprob=a)
203        except:
204            return AT_loess(condprob=None)
205
206def nth(l, n):
207    return [a[n] for a in l]
208
209class Assess(GeneSetTrans):
210    """
211    Uses the underlying GSEA code to select genes.
212    Takes data and creates attribute transformations.
213    """
214
215    def __init__(self, rankingf=None, **kwargs):
216        self.rankingf = rankingf
217        if self.rankingf == None:
218            self.rankingf = AT_edelmanParametricLearner()
219        super(Assess, self).__init__(**kwargs)
220
221    def build_features(self, data, gene_sets):
222
223        attributes = []
224
225        #attrans: { i_orig: ranking_function }
226        attrans = [ self.rankingf(iat, data) for iat, at in enumerate(data.domain.attributes) ]
227
228        nm_all, _ =  self._mat_ni(data)
229
230        for gs in gene_sets:
231
232            at = Orange.feature.Continuous(name=str(gs))
233
234            geneset = list(gs.genes)
235            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
236            genes = set(genes)
237           
238            def t(ex, w, geneset=geneset, takegenes=takegenes, nm=nm, attrans=attrans):
239
240                nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
241
242                ex_atts = [ at.name for at in ex.domain.attributes ]
243                new_atts = [ name_ind[nm.umatch(an)] if nm.umatch(an) != None else None
244                    for an in ex_atts ]
245                #new_atts: indices of genes in original data for that sample
246                #POSSIBLE REVERSE IMPLEMENTATION (slightly different
247                #for data from different chips):
248                #save pairs together and sort (or equiv. dictionary transformation)
249
250                indexes = filter(lambda x: x[0] != None, zip(new_atts, range(len(ex_atts))))
251
252                lcor = [ attrans[index_in_data](ex[index_in_ex].value) 
253                    for index_in_data, index_in_ex in indexes if
254                    ex[index_in_ex].value != '?' ]
255                #indexes in original lcor, sorted from higher to lower values
256                ordered = obiGsea.orderedPointersCorr(lcor) 
257                #subset = list of indices, lcor = correlations, ordered = order
258                subset = [ name_ind2[g] for g in genes2 ]
259                return obiGsea.enrichmentScoreRanked(subset, lcor, ordered)[0] 
260
261            at.get_value_from = t
262            attributes.append(at)
263
264        return attributes
265   
266def setSig_example_geneset(ex, data, no_unknowns):
267    """ Gets learning data and example with the same domain, both
268    containing only genes from the gene set. """
269
270    distances = [ [], [] ]   
271
272    def pearson(ex1, ex2):
273        vals1 = ex1.native(0)[:-1]
274        vals2 = ex2.native(0)[:-1]
275
276        #leaves undefined elements out
277        if not no_unknowns:
278            common = [ True if v1 != "?" and v2 != "?" else False \
279                for v1,v2 in zip(vals1,vals2) ]
280            vals1 = [ v for v,c in zip(vals1, common) if c ]
281            vals2 = [ v for v,c in zip(vals2, common) if c ]
282
283        return numpy.corrcoef([vals1, vals2])[0,1]
284
285    def ttest(ex1, ex2):
286        try:
287            return stats.lttest_ind(ex1, ex2)[0]
288        except:
289            return 0.0
290   
291    #maps class value to its index
292    classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.class_var.values) ])
293 
294    #create distances to all learning data - save or other class
295    for c in data:
296        distances[classValueMap[c[-1].value]].append(pearson(c, ex))
297
298    return ttest(distances[0], distances[1])
299
300def mat_ni(data, matcher):
301    nm = matcher([at.name for at in data.domain.attributes])
302    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes))
303    return nm, name_ind
304
305def select_genesets(nm, gene_sets, min_size=3, max_size=1000, min_part=0.1):
306    """ Returns a list of gene sets that have sizes in limits """
307
308    def ok_sizes(gs):
309        """compares sizes of genesets to limitations"""
310        transl = filter(lambda x: x != None, [ nm.umatch(gene) for gene in gs.genes ])
311        if len(transl) >= min_size \
312            and len(transl) <= max_size \
313            and float(len(transl))/len(gs.genes) >= min_part:
314            return True
315        return False
316
317    return filter(ok_sizes, gene_sets) 
318
319def vou(ex, gn, indices):
320    """ returns the value or "?" for the given gene name gn"""
321    if gn not in indices:
322        return "?"
323    else:
324        return ex[indices[gn]].value
325
326class SetSig(GeneSetTrans):
327
328    def __init__(self, **kwargs):
329        self.no_unknowns = kwargs.pop("no_unknowns", False)
330        super(SetSig, self).__init__(**kwargs)
331
332    def build_feature(self, data, gs):
333
334        at = Orange.feature.Continuous(name=str(gs))
335        geneset = list(gs.genes)
336        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
337        indices = [ name_ind[gene] for gene in genes ]
338
339        def t(ex, w, gs=gs, data=data, indices=indices, takegenes=takegenes, geneset=geneset):
340            nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
341
342            domain = Orange.data.Domain([data.domain.attributes[i] for i in indices], data.domain.class_var)
343            datao = Orange.data.Table(domain, data)
344           
345            #convert the example to the same domain
346            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
347            example = Orange.data.Instance(domain, exvalues)
348
349            return setSig_example_geneset(example, datao, self.no_unknowns) #only this one is setsig specific
350     
351        at.get_value_from = t
352        return at
353
354class ParametrizedTransformation(GeneSetTrans):
355
356    def _get_par(self, datao):
357        """ Get parameters for a subset of data, that comprises only the gene set """
358        pass
359       
360    def _use_par(self, ex, constructt):
361        pass
362   
363    def build_feature(self, data, gs):
364
365        at = Orange.feature.Continuous(name=str(gs))
366
367        geneset = list(gs.genes)
368        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
369        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
370        datao = Orange.data.Table(domain, data)
371
372        constructt = self._get_par(datao)
373
374        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
375            nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
376         
377            #convert the example to the same domain
378            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
379            ex = Orange.data.Instance(domain, exvalues)
380
381            return self._use_par(ex, constructt)
382           
383        at.get_value_from = t
384        return at
385
386class PLS(ParametrizedTransformation):
387
388    def _get_par(self, datao):
389        return PLSCall(datao, nc=1, y=[datao.domain.class_var])
390       
391    def _use_par(self, ex, constructt):
392        ex = [ ex[i].value for i in range(len(ex.domain.attributes)) ]
393        xmean, W, P, _ = constructt
394        ex = ex - xmean # same input transformation
395
396        nc = W.shape[1]
397
398        TR = numpy.empty((1, nc))
399        XR = ex
400
401        dot = numpy.dot
402
403        for i in range(nc):
404           t = dot(XR, W[:,i].T)
405           XR = XR - t*numpy.array([P[:,i]])
406           TR[:,i] = t
407
408        return TR[0][0]
409       
410class PCA(ParametrizedTransformation):
411
412    def _get_par(self, datao):
413        return pca(datao)
414
415    def _use_par(self, arr, constructt):
416        arr = [ arr[i].value for i in range(len(arr.domain.attributes)) ]
417        evals, evect, xmean = constructt
418
419        arr = arr - xmean # same input transformation - a row in a matrix
420        ev0 = evect[0] #this is a row in a matrix - do a dot product
421        a = numpy.dot(arr, ev0)
422
423        return a
424
425class SimpleFun(GeneSetTrans):
426
427    def build_feature(self, data, gs):
428
429        at = Orange.feature.Continuous(name=str(gs))
430
431        def t(ex, w, gs=gs):
432            geneset = list(gs.genes)
433            nm2, name_ind2, genes2 = self._match_instance(ex, geneset)
434           
435            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
436            exvalues = filter(lambda x: x != "?", exvalues)
437
438            return self.fn(exvalues)
439     
440        at.get_value_from = t
441        return at
442
443class Mean(SimpleFun):
444
445    def __init__(self, **kwargs):
446       self.fn = numpy.mean
447       super(Mean, self).__init__(**kwargs)
448
449class Median(SimpleFun):
450
451    def __init__(self, **kwargs):
452       self.fn = numpy.median
453       super(Median, self).__init__(**kwargs)
454
455class GSA(GeneSetTrans):
456
457    def build_features(self, data, gene_sets):
458
459        attributes = []
460
461        def tscorec(data, at, cache=None):
462            ma = obiExpression.MA_t_test()(at,data)
463            return ma
464
465        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
466
467        def to_z_score(t):
468            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
469
470        zscores = map(to_z_score, tscores)
471
472        for gs in gene_sets:
473
474            at = Orange.feature.Continuous(name=str(gs))
475
476            geneset = list(gs.genes)
477            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
478            #take each gene only once
479            genes = set(genes)
480
481            D = numpy.mean([max(zscores[name_ind[g]],0) for g in genes]) \
482                + numpy.mean([min(zscores[name_ind[g]],0) for g in genes])
483
484            if D >= 0:
485                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] > 0.0 ]
486            else:
487                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] < 0.0 ]
488
489            def t(ex, w, consider_genes=consider_genes):
490                nm2, name_ind2, genes2 = self._match_instance(ex, consider_genes)
491             
492                #convert the example to the same domain
493                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
494                exvalues = filter(lambda x: x != "?", exvalues)
495             
496                return numpy.mean(exvalues)
497
498            at.get_value_from = t
499            attributes.append(at)
500
501        return attributes
502
503def tscorec(data, at, cache=None):
504    """ Cached attribute  tscore calculation """
505    if cache != None and at in cache: return cache[at]
506    ma = obiExpression.MA_t_test()(at,data)
507    if cache != None: cache[at] = ma
508    return ma
509
510def nth(l, n):
511    return [a[n] for a in l]
512
513def compute_corg(data, inds, tscorecache):
514    """
515    Compute CORG for this geneset specified with gene inds
516    in the example table. Output is the list of gene inds
517    in CORG.
518
519    """
520    #order member genes by their t-scores: decreasing, if av(t-score) >= 0,
521    #else increasing
522    tscores = [ tscorec(data, at, tscorecache) for at in inds ]
523    sortedinds = nth(sorted(zip(inds,tscores), key=lambda x: x[1], \
524        reverse=numpy.mean(tscores) >= 0), 0)
525
526    def S(corg):
527        """ Activity score separation - S(G) in
528        the article """
529        asv = Orange.feature.Continuous(name='AS')
530        asv.getValueFrom = lambda ex,rw: Orange.data.Value(asv, corgs_activity_score(ex, corg))
531        data2 = Orange.data.Table(Orange.data.Domain([asv], data.domain.classVar), data)
532        return abs(tscorec(data2, 0)) #FIXME absolute - nothing in the article abs()
533           
534    #greedily find CORGS procing the best separation
535    g = S(sortedinds[:1])
536    bg = 1
537    for a in range(2, len(sortedinds)+1):
538        tg = S(sortedinds[:a])
539        if tg > g:
540            g = tg
541            bg = a
542        else:
543            break
544       
545    return sortedinds[:a]
546
547class CORGs(ParametrizedTransformation):
548    """
549    WARNING: input has to be z_ij table! each gene needs to be normalized
550    (mean=0, stdev=1) for all samples.
551    """
552
553    def __call__(self, *args, **kwargs):
554        self.tscorecache = {} #reset a cache
555        return super(CORGs, self).__call__(*args, **kwargs)
556
557    def build_feature(self, data, gs):
558
559        at = Orange.feature.Continuous(name=str(gs))
560        geneset = list(gs.genes)
561
562        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
563        indices = compute_corg(data, [ name_ind[g] for g in genes ], self.tscorecache)
564
565        ind_names = dict( (a,b) for b,a in name_ind.items() )
566        selected_genes = sorted(set([to_geneset[ind_names[i]] for i in indices]))
567           
568        def t(ex, w, corg=selected_genes): #copy od the data
569            nm2, name_ind2, genes2 = self._match_instance(ex, corg, None)
570            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ]
571            return sum(v if v != '?' else 0.0 for v in exvalues)/len(corg)**0.5
572     
573        at.get_value_from = t
574        return at
575
576if __name__ == "__main__":
577
578    data = Orange.data.Table("iris")
579    gsets = obiGeneSets.collections({
580        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
581        "f3": ['sepal length', 'sepal width', 'petal length'],
582        "l3": ['sepal width', 'petal length', 'petal width'],
583        })
584    matcher = obiGene.matcher([])
585    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
586
587    """
588    data = Orange.data.Table("DLBCL_200a")
589    gsets = obiGeneSets.collections((("KEGG",),"9606"))
590    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
591    choosen_cv = None
592    """
593
594    def to_old_dic(d, data):
595        ar = defaultdict(list)
596        for ex1 in data:
597            ex = d(ex1)
598            for a,x in zip(d.attributes, ex):
599                ar[a.name].append(x.value)
600        return ar
601
602    def pp2(ar):
603        ol =  sorted(ar.items())
604        print '\n'.join([ a + ": " +str(b) for a,b in ol])
605
606    ass = Assess(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
607    ar = to_old_dic(ass.domain, data[:5])
608    pp2(ar)
Note: See TracBrowser for help on using the repository browser.