source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1636:10d234fdadb9

Revision 1636:10d234fdadb9, 18.2 KB checked in by mitar, 2 years ago (diff)

Restructuring because we will not be using namespaces.

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