source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1728:a7d42a8d1ed2

Revision 1728:a7d42a8d1ed2, 20.9 KB checked in by markotoplak, 14 months ago (diff)

Added an option in Assess to not ignore umatchable context.

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