source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1690:cba131a6a0c5

Revision 1690:cba131a6a0c5, 18.4 KB checked in by markotoplak, 22 months ago (diff)

obiGeneSetSig: small changes in the interface.

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