source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1708:c442e936f1ef

Revision 1708:c442e936f1ef, 19.6 KB checked in by markotoplak, 20 months ago (diff)

Another speedup of SetSig. Now it is 30% faster on DLBCL.tab with no_unknowns=True.

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        takegenes = [ geneset[i] for i in takegenes ]
339
340        def t(ex, w, gs=gs, data=data, indices=indices, takegenes=takegenes):
341            nm2, name_ind2, genes2 = self._match_instance(ex, takegenes)
342
343            domain = Orange.data.Domain([data.domain.attributes[i] for i in indices], data.domain.class_var)
344            datao = Orange.data.Table(domain, data)
345           
346            #convert the example to the same domain
347            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
348            example = Orange.data.Instance(domain, exvalues)
349
350            return setSig_example_geneset(example, datao, self.no_unknowns) #only this one is setsig specific
351     
352        at.get_value_from = t
353        return at
354
355class ParametrizedTransformation(GeneSetTrans):
356
357    def _get_par(self, datao):
358        """ Get parameters for a subset of data, that comprises only the gene set """
359        pass
360       
361    def _use_par(self, ex, constructt):
362        pass
363   
364    def build_feature(self, data, gs):
365
366        at = Orange.feature.Continuous(name=str(gs))
367
368        geneset = list(gs.genes)
369        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
370        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
371        datao = Orange.data.Table(domain, data)
372
373        constructt = self._get_par(datao)
374
375        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
376            nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
377         
378            #convert the example to the same domain
379            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
380            ex = Orange.data.Instance(domain, exvalues)
381
382            return self._use_par(ex, constructt)
383           
384        at.get_value_from = t
385        return at
386
387class PLS(ParametrizedTransformation):
388
389    def _get_par(self, datao):
390        return PLSCall(datao, nc=1, y=[datao.domain.class_var])
391       
392    def _use_par(self, ex, constructt):
393        ex = [ ex[i].value for i in range(len(ex.domain.attributes)) ]
394        xmean, W, P, _ = constructt
395        ex = ex - xmean # same input transformation
396
397        nc = W.shape[1]
398
399        TR = numpy.empty((1, nc))
400        XR = ex
401
402        dot = numpy.dot
403
404        for i in range(nc):
405           t = dot(XR, W[:,i].T)
406           XR = XR - t*numpy.array([P[:,i]])
407           TR[:,i] = t
408
409        return TR[0][0]
410       
411class PCA(ParametrizedTransformation):
412
413    def _get_par(self, datao):
414        return pca(datao)
415
416    def _use_par(self, arr, constructt):
417        arr = [ arr[i].value for i in range(len(arr.domain.attributes)) ]
418        evals, evect, xmean = constructt
419
420        arr = arr - xmean # same input transformation - a row in a matrix
421        ev0 = evect[0] #this is a row in a matrix - do a dot product
422        a = numpy.dot(arr, ev0)
423
424        return a
425
426class SimpleFun(GeneSetTrans):
427
428    def build_feature(self, data, gs):
429
430        at = Orange.feature.Continuous(name=str(gs))
431
432        def t(ex, w, gs=gs):
433            geneset = list(gs.genes)
434            nm2, name_ind2, genes2 = self._match_instance(ex, geneset)
435           
436            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
437            exvalues = filter(lambda x: x != "?", exvalues)
438
439            return self.fn(exvalues)
440     
441        at.get_value_from = t
442        return at
443
444class Mean(SimpleFun):
445
446    def __init__(self, **kwargs):
447       self.fn = numpy.mean
448       super(Mean, self).__init__(**kwargs)
449
450class Median(SimpleFun):
451
452    def __init__(self, **kwargs):
453       self.fn = numpy.median
454       super(Median, self).__init__(**kwargs)
455
456class GSA(GeneSetTrans):
457
458    def build_features(self, data, gene_sets):
459
460        attributes = []
461
462        def tscorec(data, at, cache=None):
463            ma = obiExpression.MA_t_test()(at,data)
464            return ma
465
466        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
467
468        def to_z_score(t):
469            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
470
471        zscores = map(to_z_score, tscores)
472
473        for gs in gene_sets:
474
475            at = Orange.feature.Continuous(name=str(gs))
476
477            geneset = list(gs.genes)
478            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
479            #take each gene only once
480            genes = set(genes)
481
482            D = numpy.mean([max(zscores[name_ind[g]],0) for g in genes]) \
483                + numpy.mean([min(zscores[name_ind[g]],0) for g in genes])
484
485            if D >= 0:
486                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] > 0.0 ]
487            else:
488                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] < 0.0 ]
489
490            def t(ex, w, consider_genes=consider_genes):
491                nm2, name_ind2, genes2 = self._match_instance(ex, consider_genes)
492             
493                #convert the example to the same domain
494                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
495                exvalues = filter(lambda x: x != "?", exvalues)
496             
497                return numpy.mean(exvalues)
498
499            at.get_value_from = t
500            attributes.append(at)
501
502        return attributes
503
504def tscorec(data, at, cache=None):
505    """ Cached attribute  tscore calculation """
506    if cache != None and at in cache: return cache[at]
507    ma = obiExpression.MA_t_test()(at,data)
508    if cache != None: cache[at] = ma
509    return ma
510
511def nth(l, n):
512    return [a[n] for a in l]
513
514def compute_corg(data, inds, tscorecache):
515    """
516    Compute CORG for this geneset specified with gene inds
517    in the example table. Output is the list of gene inds
518    in CORG.
519
520    """
521    #order member genes by their t-scores: decreasing, if av(t-score) >= 0,
522    #else increasing
523    tscores = [ tscorec(data, at, tscorecache) for at in inds ]
524    sortedinds = nth(sorted(zip(inds,tscores), key=lambda x: x[1], \
525        reverse=numpy.mean(tscores) >= 0), 0)
526
527    def S(corg):
528        """ Activity score separation - S(G) in
529        the article """
530        asv = Orange.feature.Continuous(name='AS')
531        asv.getValueFrom = lambda ex,rw: Orange.data.Value(asv, corgs_activity_score(ex, corg))
532        data2 = Orange.data.Table(Orange.data.Domain([asv], data.domain.classVar), data)
533        return abs(tscorec(data2, 0)) #FIXME absolute - nothing in the article abs()
534           
535    #greedily find CORGS procing the best separation
536    g = S(sortedinds[:1])
537    bg = 1
538    for a in range(2, len(sortedinds)+1):
539        tg = S(sortedinds[:a])
540        if tg > g:
541            g = tg
542            bg = a
543        else:
544            break
545       
546    return sortedinds[:a]
547
548class CORGs(ParametrizedTransformation):
549    """
550    WARNING: input has to be z_ij table! each gene needs to be normalized
551    (mean=0, stdev=1) for all samples.
552    """
553
554    def __call__(self, *args, **kwargs):
555        self.tscorecache = {} #reset a cache
556        return super(CORGs, self).__call__(*args, **kwargs)
557
558    def build_feature(self, data, gs):
559
560        at = Orange.feature.Continuous(name=str(gs))
561        geneset = list(gs.genes)
562
563        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
564        indices = compute_corg(data, [ name_ind[g] for g in genes ], self.tscorecache)
565
566        ind_names = dict( (a,b) for b,a in name_ind.items() )
567        selected_genes = sorted(set([to_geneset[ind_names[i]] for i in indices]))
568           
569        def t(ex, w, corg=selected_genes): #copy od the data
570            nm2, name_ind2, genes2 = self._match_instance(ex, corg, None)
571            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ]
572            return sum(v if v != '?' else 0.0 for v in exvalues)/len(corg)**0.5
573     
574        at.get_value_from = t
575        return at
576
577if __name__ == "__main__":
578
579    data = Orange.data.Table("iris")
580    gsets = obiGeneSets.collections({
581        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
582        "f3": ['sepal length', 'sepal width', 'petal length'],
583        "l3": ['sepal width', 'petal length', 'petal width'],
584        })
585    matcher = obiGene.matcher([])
586    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
587
588    """
589    data = Orange.data.Table("DLBCL_200a")
590    gsets = obiGeneSets.collections((("KEGG",),"9606"))
591    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
592    choosen_cv = None
593    """
594
595    def to_old_dic(d, data):
596        ar = defaultdict(list)
597        for ex1 in data:
598            ex = d(ex1)
599            for a,x in zip(d.attributes, ex):
600                ar[a.name].append(x.value)
601        return ar
602
603    def pp2(ar):
604        ol =  sorted(ar.items())
605        print '\n'.join([ a + ": " +str(b) for a,b in ol])
606
607    ass = Assess(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
608    ar = to_old_dic(ass.domain, data[:5])
609    pp2(ar)
Note: See TracBrowser for help on using the repository browser.