source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1705:0b6781662f02

Revision 1705:0b6781662f02, 19.4 KB checked in by markotoplak, 21 months ago (diff)

obiGsea.takeClasses: now builds a class value with get_value_from.

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