source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1723:753e1d554488

Revision 1723:753e1d554488, 20.4 KB checked in by markotoplak, 18 months ago (diff)

SetSig uses statc for pearsonr. Around 10x speedup compared to numpy.

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