source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1691:d2aeb3889a90

Revision 1691:d2aeb3889a90, 19.4 KB checked in by markotoplak, 22 months ago (diff)

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