source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1711:643d37885055

Revision 1711:643d37885055, 20.3 KB checked in by markotoplak, 20 months ago (diff)

obiGeneSetSig: speedup of ASSESS.

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        self.example_buffer = {}
220        self.attransv = 0
221        super(Assess, self).__init__(**kwargs)
222
223    def _ordered_and_lcor(self, ex, nm, name_ind, attrans, attransv):
224        """ Buffered! It should be computed only once per example. """ 
225        #name_ind and nm are always co-created, so I need to have only one as a key
226        key = (ex, nm, attransv)
227        if key not in self.example_buffer:
228            ex_atts = [ at.name for at in ex.domain.attributes ]
229            new_atts = [ name_ind[nm.umatch(an)] if nm.umatch(an) != None else None
230                for an in ex_atts ]
231
232            #new_atts: indices of genes in original data for that sample
233            #POSSIBLE REVERSE IMPLEMENTATION (slightly different
234            #for data from different chips):
235            #save pairs together and sort (or equiv. dictionary transformation)
236
237            indexes = filter(lambda x: x[0] != None, zip(new_atts, range(len(ex_atts))))
238
239            lcor = [ attrans[index_in_data](ex[index_in_ex].value) 
240                for index_in_data, index_in_ex in indexes if
241                ex[index_in_ex].value != '?' ]
242            #indexes in original lcor, sorted from higher to lower values
243            ordered = obiGsea.orderedPointersCorr(lcor)
244            rev2 = numpy.argsort(ordered)
245            self.example_buffer[key] = lcor,ordered,rev2
246        return self.example_buffer[key]
247
248    def build_features(self, data, gene_sets):
249
250        attributes = []
251
252        #attrans: { i_orig: ranking_function }
253        attrans = [ self.rankingf(iat, data) for iat, at in enumerate(data.domain.attributes) ]
254        attransv = self.attransv
255        self.attransv += 1
256
257        nm_all, _ =  self._mat_ni(data)
258
259        for gs in gene_sets:
260
261            at = Orange.feature.Continuous(name=str(gs))
262
263            geneset = list(gs.genes)
264            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
265            takegenes = [ geneset[i] for i in takegenes ]
266            genes = set(genes)
267
268            def t(ex, w, takegenes=takegenes, nm=nm, attrans=attrans, attransv=attransv):
269
270                nm2, name_ind2, genes2 = self._match_instance(ex, takegenes)
271                lcor, ordered, rev2 = self._ordered_and_lcor(ex, nm, name_ind, attrans, attransv)
272
273                #subset = list of indices, lcor = correlations, ordered = order
274                subset = [ name_ind2[g] for g in genes2 ]
275                return obiGsea.enrichmentScoreRanked(subset, lcor, ordered, rev2=rev2)[0] 
276
277            at.get_value_from = t
278            attributes.append(at)
279
280        return attributes
281   
282def setSig_example_geneset(ex, data, no_unknowns):
283    """ Gets learning data and example with the same domain, both
284    containing only genes from the gene set. """
285
286    distances = [ [], [] ]   
287
288    def pearson(ex1, ex2):
289        vals1 = ex1.native(0)[:-1]
290        vals2 = ex2.native(0)[:-1]
291
292        #leaves undefined elements out
293        if not no_unknowns:
294            common = [ True if v1 != "?" and v2 != "?" else False \
295                for v1,v2 in zip(vals1,vals2) ]
296            vals1 = [ v for v,c in zip(vals1, common) if c ]
297            vals2 = [ v for v,c in zip(vals2, common) if c ]
298
299        return numpy.corrcoef([vals1, vals2])[0,1]
300
301    def ttest(ex1, ex2):
302        try:
303            return stats.lttest_ind(ex1, ex2)[0]
304        except:
305            return 0.0
306   
307    #maps class value to its index
308    classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.class_var.values) ])
309 
310    #create distances to all learning data - save or other class
311    for c in data:
312        distances[classValueMap[c[-1].value]].append(pearson(c, ex))
313
314    return ttest(distances[0], distances[1])
315
316def mat_ni(data, matcher):
317    nm = matcher([at.name for at in data.domain.attributes])
318    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes))
319    return nm, name_ind
320
321def select_genesets(nm, gene_sets, min_size=3, max_size=1000, min_part=0.1):
322    """ Returns a list of gene sets that have sizes in limits """
323
324    def ok_sizes(gs):
325        """compares sizes of genesets to limitations"""
326        transl = filter(lambda x: x != None, [ nm.umatch(gene) for gene in gs.genes ])
327        if len(transl) >= min_size \
328            and len(transl) <= max_size \
329            and float(len(transl))/len(gs.genes) >= min_part:
330            return True
331        return False
332
333    return filter(ok_sizes, gene_sets) 
334
335def vou(ex, gn, indices):
336    """ returns the value or "?" for the given gene name gn"""
337    if gn not in indices:
338        return "?"
339    else:
340        return ex[indices[gn]].value
341
342class SetSig(GeneSetTrans):
343
344    def __init__(self, **kwargs):
345        self.no_unknowns = kwargs.pop("no_unknowns", False)
346        super(SetSig, self).__init__(**kwargs)
347
348    def build_feature(self, data, gs):
349
350        at = Orange.feature.Continuous(name=str(gs))
351        geneset = list(gs.genes)
352        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
353        indices = [ name_ind[gene] for gene in genes ]
354        takegenes = [ geneset[i] for i in takegenes ]
355
356        def t(ex, w, gs=gs, data=data, indices=indices, takegenes=takegenes):
357            nm2, name_ind2, genes2 = self._match_instance(ex, takegenes)
358
359            domain = Orange.data.Domain([data.domain.attributes[i] for i in indices], data.domain.class_var)
360            datao = Orange.data.Table(domain, data)
361           
362            #convert the example to the same domain
363            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
364            example = Orange.data.Instance(domain, exvalues)
365
366            return setSig_example_geneset(example, datao, self.no_unknowns) #only this one is setsig specific
367     
368        at.get_value_from = t
369        return at
370
371class ParametrizedTransformation(GeneSetTrans):
372
373    def _get_par(self, datao):
374        """ Get parameters for a subset of data, that comprises only the gene set """
375        pass
376       
377    def _use_par(self, ex, constructt):
378        pass
379   
380    def build_feature(self, data, gs):
381
382        at = Orange.feature.Continuous(name=str(gs))
383
384        geneset = list(gs.genes)
385        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
386        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
387        datao = Orange.data.Table(domain, data)
388        takegenes = [ geneset[i] for i in takegenes ]
389
390        constructt = self._get_par(datao)
391
392        def t(ex, w, constructt=constructt, takegenes=takegenes, domain=domain):
393            nm2, name_ind2, genes2 = self._match_instance(ex, takegenes)
394         
395            #convert the example to the same domain
396            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
397            ex = Orange.data.Instance(domain, exvalues)
398
399            return self._use_par(ex, constructt)
400           
401        at.get_value_from = t
402        return at
403
404class PLS(ParametrizedTransformation):
405
406    def _get_par(self, datao):
407        return PLSCall(datao, nc=1, y=[datao.domain.class_var])
408       
409    def _use_par(self, ex, constructt):
410        ex = [ ex[i].value for i in range(len(ex.domain.attributes)) ]
411        xmean, W, P, _ = constructt
412        ex = ex - xmean # same input transformation
413
414        nc = W.shape[1]
415
416        TR = numpy.empty((1, nc))
417        XR = ex
418
419        dot = numpy.dot
420
421        for i in range(nc):
422           t = dot(XR, W[:,i].T)
423           XR = XR - t*numpy.array([P[:,i]])
424           TR[:,i] = t
425
426        return TR[0][0]
427       
428class PCA(ParametrizedTransformation):
429
430    def _get_par(self, datao):
431        return pca(datao)
432
433    def _use_par(self, arr, constructt):
434        arr = [ arr[i].value for i in range(len(arr.domain.attributes)) ]
435        evals, evect, xmean = constructt
436
437        arr = arr - xmean # same input transformation - a row in a matrix
438        ev0 = evect[0] #this is a row in a matrix - do a dot product
439        a = numpy.dot(arr, ev0)
440
441        return a
442
443class SimpleFun(GeneSetTrans):
444
445    def build_feature(self, data, gs):
446
447        at = Orange.feature.Continuous(name=str(gs))
448
449        def t(ex, w, gs=gs):
450            geneset = list(gs.genes)
451            nm2, name_ind2, genes2 = self._match_instance(ex, geneset)
452           
453            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
454            exvalues = filter(lambda x: x != "?", exvalues)
455
456            return self.fn(exvalues)
457     
458        at.get_value_from = t
459        return at
460
461class Mean(SimpleFun):
462
463    def __init__(self, **kwargs):
464       self.fn = numpy.mean
465       super(Mean, self).__init__(**kwargs)
466
467class Median(SimpleFun):
468
469    def __init__(self, **kwargs):
470       self.fn = numpy.median
471       super(Median, self).__init__(**kwargs)
472
473class GSA(GeneSetTrans):
474
475    def build_features(self, data, gene_sets):
476
477        attributes = []
478
479        def tscorec(data, at, cache=None):
480            ma = obiExpression.MA_t_test()(at,data)
481            return ma
482
483        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
484
485        def to_z_score(t):
486            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
487
488        zscores = map(to_z_score, tscores)
489
490        for gs in gene_sets:
491
492            at = Orange.feature.Continuous(name=str(gs))
493
494            geneset = list(gs.genes)
495            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
496            #take each gene only once
497            genes = set(genes)
498
499            D = numpy.mean([max(zscores[name_ind[g]],0) for g in genes]) \
500                + numpy.mean([min(zscores[name_ind[g]],0) for g in genes])
501
502            if D >= 0:
503                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] > 0.0 ]
504            else:
505                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] < 0.0 ]
506
507            def t(ex, w, consider_genes=consider_genes):
508                nm2, name_ind2, genes2 = self._match_instance(ex, consider_genes)
509             
510                #convert the example to the same domain
511                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
512                exvalues = filter(lambda x: x != "?", exvalues)
513             
514                return numpy.mean(exvalues)
515
516            at.get_value_from = t
517            attributes.append(at)
518
519        return attributes
520
521def tscorec(data, at, cache=None):
522    """ Cached attribute  tscore calculation """
523    if cache != None and at in cache: return cache[at]
524    ma = obiExpression.MA_t_test()(at,data)
525    if cache != None: cache[at] = ma
526    return ma
527
528def nth(l, n):
529    return [a[n] for a in l]
530
531def compute_corg(data, inds, tscorecache):
532    """
533    Compute CORG for this geneset specified with gene inds
534    in the example table. Output is the list of gene inds
535    in CORG.
536
537    """
538    #order member genes by their t-scores: decreasing, if av(t-score) >= 0,
539    #else increasing
540    tscores = [ tscorec(data, at, tscorecache) for at in inds ]
541    sortedinds = nth(sorted(zip(inds,tscores), key=lambda x: x[1], \
542        reverse=numpy.mean(tscores) >= 0), 0)
543
544    def S(corg):
545        """ Activity score separation - S(G) in
546        the article """
547        asv = Orange.feature.Continuous(name='AS')
548        asv.getValueFrom = lambda ex,rw: Orange.data.Value(asv, corgs_activity_score(ex, corg))
549        data2 = Orange.data.Table(Orange.data.Domain([asv], data.domain.classVar), data)
550        return abs(tscorec(data2, 0)) #FIXME absolute - nothing in the article abs()
551           
552    #greedily find CORGS procing the best separation
553    g = S(sortedinds[:1])
554    bg = 1
555    for a in range(2, len(sortedinds)+1):
556        tg = S(sortedinds[:a])
557        if tg > g:
558            g = tg
559            bg = a
560        else:
561            break
562       
563    return sortedinds[:a]
564
565class CORGs(ParametrizedTransformation):
566    """
567    WARNING: input has to be z_ij table! each gene needs to be normalized
568    (mean=0, stdev=1) for all samples.
569    """
570
571    def __call__(self, *args, **kwargs):
572        self.tscorecache = {} #reset a cache
573        return super(CORGs, self).__call__(*args, **kwargs)
574
575    def build_feature(self, data, gs):
576
577        at = Orange.feature.Continuous(name=str(gs))
578        geneset = list(gs.genes)
579
580        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
581        indices = compute_corg(data, [ name_ind[g] for g in genes ], self.tscorecache)
582
583        ind_names = dict( (a,b) for b,a in name_ind.items() )
584        selected_genes = sorted(set([to_geneset[ind_names[i]] for i in indices]))
585           
586        def t(ex, w, corg=selected_genes): #copy od the data
587            nm2, name_ind2, genes2 = self._match_instance(ex, corg, None)
588            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ]
589            return sum(v if v != '?' else 0.0 for v in exvalues)/len(corg)**0.5
590     
591        at.get_value_from = t
592        return at
593
594if __name__ == "__main__":
595
596    data = Orange.data.Table("iris")
597    gsets = obiGeneSets.collections({
598        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
599        "f3": ['sepal length', 'sepal width', 'petal length'],
600        "l3": ['sepal width', 'petal length', 'petal width'],
601        })
602    matcher = obiGene.matcher([])
603    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
604
605    """
606    data = Orange.data.Table("DLBCL_200a")
607    gsets = obiGeneSets.collections((("KEGG",),"9606"))
608    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
609    choosen_cv = None
610    """
611
612    def to_old_dic(d, data):
613        ar = defaultdict(list)
614        for ex1 in data:
615            ex = d(ex1)
616            for a,x in zip(d.attributes, ex):
617                ar[a.name].append(x.value)
618        return ar
619
620    def pp2(ar):
621        ol =  sorted(ar.items())
622        print '\n'.join([ a + ": " +str(b) for a,b in ol])
623
624    ass = Assess(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
625    ar = to_old_dic(ass.domain, data[:5])
626    pp2(ar)
Note: See TracBrowser for help on using the repository browser.