source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1725:1c834422f29e

Revision 1725:1c834422f29e, 20.8 KB checked in by markotoplak, 17 months ago (diff)

obiGeneSetSig: assess now works even if gene matcher did not perfectly match the whole file.

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