source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1774:d0585dc25172

Revision 1774:d0585dc25172, 25.0 KB checked in by markotoplak, 12 months ago (diff)

obiGeneSetSig: added another method based on likelihood ratios by Su et al (2009).

RevLine 
[1632]1from __future__ import absolute_import
2
3import math
4from collections import defaultdict
5
6import scipy.stats
7
[1572]8import numpy
[1632]9
10import Orange, Orange.utils, statc
11
[1707]12if __name__ == "__main__":
13    __package__ = "Orange.bio"
14
[1632]15from .obiGsea import takeClasses
16from .obiAssess import pca, PLSCall, corgs_activity_score
17from . import obiExpression, obiGene, obiGeneSets, obiGsea, stats
[1597]18
[1607]19class GeneSetTrans(object):
[1572]20
[1607]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):
[1774]30        """
31        Return
32            - a gene matcher with the instance as a target
33            - { name: attribute indices } of an instance
34            - genes names on the data set that were matched by the gene set
35
36        If takegenes is a list of indices, use only genes from
37        the gene set with specified indices.
38        """
[1607]39        nm, name_ind = self._mat_ni(instance)
40        genes = [ nm.umatch(gene) for gene in geneset ]
41        if takegenes:
42            genes = [ genes[i] for i in takegenes ]
43        return nm, name_ind, genes
44
45    def _match_data(self, data, geneset, odic=False):
46        nm, name_ind = self._mat_ni(data)
47        genes = [ nm.umatch(gene) for gene in geneset ]
48        if odic:
49            to_geneset = dict(zip(genes, geneset))
50        takegenes = [ i for i,a in enumerate(genes) if a != None ]
51        genes = [ genes[i] for i in takegenes ]
52        if odic:
53            return nm, name_ind, genes, takegenes, to_geneset
54        else:
55            return nm, name_ind, genes, takegenes
56
[1691]57    def __init__(self, matcher=None, gene_sets=None, min_size=3, max_size=1000, min_part=0.1, class_values=None, cv=False):
[1607]58        self.matcher = matcher
59        self.gene_sets = gene_sets
60        self.min_size = min_size
61        self.max_size = max_size
62        self.min_part = min_part
63        self.class_values = class_values
64        self._cache = {}
[1691]65        self.cv = cv
[1607]66
67    def __call__(self, data, weight_id=None):
68
69        #selection of classes and gene sets
70        data = takeClasses(data, classValues=self.class_values)
71        nm,_ =  self._mat_ni(data)
72        gene_sets = select_genesets(nm, self.gene_sets, self.min_size, self.max_size, self.min_part)
73
74        #build a new domain
75        newfeatures = self.build_features(data, gene_sets)
76        newdomain = Orange.data.Domain(newfeatures, data.domain.class_var)
[1691]77
78        #build a data set with cross validation
79        if self.cv == False:
80            return Orange.data.Table(newdomain, data)
81        else:
82            # The domain has the transformer that is build on all samples,
83            # while the transformed data table uses cross-validation
84            # internally
[1772]85            if self.cv == True:
86                cvi = Orange.data.sample.SubsetIndicesCV(data, 5)
87            elif self.cv != False:
88                cvi = self.cv(data)
[1691]89            data_cv = [ [] for _ in range(len(data)) ]
[1772]90            for f in set(cvi):
[1691]91                learn = data.select(cvi, f, negate=True)
92                test = data.select(cvi, f)
93                lf = self.build_features(learn, gene_sets)
94                transd = Orange.data.Domain(lf, data.domain.class_var)
95                trans_test = Orange.data.Table(transd, test)
96                for ex, pos in \
97                    zip(trans_test, [ i for i,n in enumerate(cvi) if n == f ]):
98                    data_cv[pos] = ex.native(0)
99            return Orange.data.Table(newdomain, data_cv)
[1607]100
101    def build_features(self, data, gene_sets):
102        return [ self.build_feature(data, gs) for gs in gene_sets ]
103
104def normcdf(x, mi, st):
105    return 0.5*(2. - stats.erfcc((x - mi)/(st*math.sqrt(2))))
106
107class AT_edelmanParametric(object):
108
109    def __init__(self, **kwargs):
110        for a,b in kwargs.items():
111            setattr(self, a, b)
112
113    def __call__(self, nval):
114
115        if self.mi1 == None or self.mi2 == None or self.st1 == None or self.st2 == None:
116            return 0.0 
117
118        val = nval
119
120        try:
121            if val >= self.mi1:
122                p1 = 1 - normcdf(val, self.mi1, self.st1)
123            else:
124                p1 = normcdf(val, self.mi1, self.st1)
125
126            if val >= self.mi2:
127                p2 = 1 - normcdf(val, self.mi2, self.st2)
128            else:
129                p2 = normcdf(val, self.mi2, self.st2)
130
131            #print p1, p2
132            return math.log(p1/p2)
133        except:
134            #print p1, p2, "exception"
135            return 0
136
[1774]137def estimate_gaussian_per_class(data, i, a=None, b=None):
138    cv = data.domain.class_var
139
140    if a == None: a = cv.values[0]
141    if b == None: 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(a)
147    list2 = avWCVal(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 mi1, st1, mi2, st2
164
[1607]165class AT_edelmanParametricLearner(object):
166    """
167    Returns attribute transfromer for Edelman parametric measure for a
168    given attribute in the dataset.  Edelman et al, 06. Modified a bit.
169    """
170
171    def __init__(self, a=None, b=None):
172        """
173        a and b are choosen class values.
174        """
175        self.a = a
176        self.b = b
177
178    def __call__(self, i, data):
[1774]179        cv = data.domain.class_var
[1607]180
181        if self.a == None: self.a = cv.values[0]
182        if self.b == None: self.b = cv.values[1]
183
[1774]184        mi1, st1, mi2, st2 = estimate_gaussian_per_class(data, i, a=self.a, b=self.b)
[1607]185
186        return AT_edelmanParametric(mi1=mi1, mi2=mi2, st1=st1, st2=st2)
187
188class AT_loess(object):
189
190    def __init__(self, **kwargs):
191        for a,b in kwargs.items():
192            setattr(self, a, b)
193
194    def __call__(self, nval):
195
196        val = nval
197
198        def saveplog(a,b):
199            try:
200                return math.log(a/b)
201            except:
202                if a < b:
203                    return -10
204                else:
205                    return +10
206
207        try:
208            ocene = self.condprob(val)
209            if sum(ocene) < 0.01:
210                return 0.0
211            return saveplog(ocene[0], ocene[1])
212
213        except:
214            return 0.0
215
216class AT_loessLearner(object):
217
218    def __call__(self, i, data):
219        try:
220            ca = Orange.statistics.contingency.VarClass(data.domain.attributes[i], data)
221            a =  Orange.statistics.estimate.ConditionalLoess(ca, nPoints=5)
222            return AT_loess(condprob=a)
223        except:
224            return AT_loess(condprob=None)
225
226def nth(l, n):
227    return [a[n] for a in l]
228
229class Assess(GeneSetTrans):
230    """
231    Uses the underlying GSEA code to select genes.
232    Takes data and creates attribute transformations.
233    """
234
235    def __init__(self, rankingf=None, **kwargs):
236        self.rankingf = rankingf
237        if self.rankingf == None:
238            self.rankingf = AT_edelmanParametricLearner()
[1711]239        self.example_buffer = {}
240        self.attransv = 0
[1728]241        self.ignore_unmatchable_context = True
[1607]242        super(Assess, self).__init__(**kwargs)
243
[1711]244    def _ordered_and_lcor(self, ex, nm, name_ind, attrans, attransv):
245        """ Buffered! It should be computed only once per example. """ 
246        #name_ind and nm are always co-created, so I need to have only one as a key
247        key = (ex, nm, attransv)
248        if key not in self.example_buffer:
249            ex_atts = [ at.name for at in ex.domain.attributes ]
[1728]250            new_atts = [ name_ind[nm.umatch(an)] if nm.umatch(an) != None else (None if self.ignore_unmatchable_context else i)
251                for i,an in enumerate(ex_atts) ]
[1711]252
253            #new_atts: indices of genes in original data for that sample
254            #POSSIBLE REVERSE IMPLEMENTATION (slightly different
255            #for data from different chips):
256            #save pairs together and sort (or equiv. dictionary transformation)
257
258            indexes = filter(lambda x: x[0] != None, zip(new_atts, range(len(ex_atts))))
259
260            lcor = [ attrans[index_in_data](ex[index_in_ex].value) 
261                for index_in_data, index_in_ex in indexes if
262                ex[index_in_ex].value != '?' ]
[1725]263
264            indices_to_lcori = dict( (index_in_ex, i) for i,(_, index_in_ex) in enumerate(indexes) 
265                if ex[index_in_ex].value != '?')
266
[1711]267            #indexes in original lcor, sorted from higher to lower values
268            ordered = obiGsea.orderedPointersCorr(lcor)
269            rev2 = numpy.argsort(ordered)
[1725]270            self.example_buffer[key] = lcor, ordered, rev2, indices_to_lcori
[1711]271        return self.example_buffer[key]
272
[1607]273    def build_features(self, data, gene_sets):
274
275        attributes = []
276
277        #attrans: { i_orig: ranking_function }
278        attrans = [ self.rankingf(iat, data) for iat, at in enumerate(data.domain.attributes) ]
[1711]279        attransv = self.attransv
280        self.attransv += 1
[1607]281
282        nm_all, _ =  self._mat_ni(data)
283
284        for gs in gene_sets:
285
286            at = Orange.feature.Continuous(name=str(gs))
287
288            geneset = list(gs.genes)
289            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
[1711]290            takegenes = [ geneset[i] for i in takegenes ]
[1607]291            genes = set(genes)
292
[1711]293            def t(ex, w, takegenes=takegenes, nm=nm, attrans=attrans, attransv=attransv):
294                nm2, name_ind2, genes2 = self._match_instance(ex, takegenes)
[1725]295                lcor, ordered, rev2, indices_to_lcori = \
296                    self._ordered_and_lcor(ex, nm, name_ind, attrans, attransv)
[1728]297
298           
[1607]299                #subset = list of indices, lcor = correlations, ordered = order
[1725]300                #make it compatible with lcor, if some are missing in lcor
301                subset = filter(None,
302                    [ indices_to_lcori.get(name_ind2[g], None) for g in genes2 ] )
[1728]303
[1711]304                return obiGsea.enrichmentScoreRanked(subset, lcor, ordered, rev2=rev2)[0] 
[1607]305
306            at.get_value_from = t
307            attributes.append(at)
308
309        return attributes
310   
[1772]311def setSig_example_geneset(ex, data, no_unknowns, check_same=False):
[1583]312    """ Gets learning data and example with the same domain, both
313    containing only genes from the gene set. """
314
315    distances = [ [], [] ]   
316
[1707]317    def pearson(ex1, ex2):
318        vals1 = ex1.native(0)[:-1]
319        vals2 = ex2.native(0)[:-1]
[1583]320
[1772]321        if check_same and vals1 == vals2:
322            return 10 #they are the same
323
[1583]324        #leaves undefined elements out
[1707]325        if not no_unknowns:
326            common = [ True if v1 != "?" and v2 != "?" else False \
327                for v1,v2 in zip(vals1,vals2) ]
328            vals1 = [ v for v,c in zip(vals1, common) if c ]
329            vals2 = [ v for v,c in zip(vals2, common) if c ]
[1583]330
[1723]331        #statc correlation is from 5-10 times faster than numpy!
332        try:
333            return statc.pearsonr(vals1, vals2)[0]
334        except:
335            return numpy.corrcoef([vals1, vals2])[0,1] 
336       
[1583]337
338    def ttest(ex1, ex2):
339        try:
340            return stats.lttest_ind(ex1, ex2)[0]
341        except:
342            return 0.0
343   
344    #maps class value to its index
345    classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.class_var.values) ])
346 
347    #create distances to all learning data - save or other class
348    for c in data:
[1772]349        p = pearson(c, ex)
350        if p != 10:
351             distances[classValueMap[c[-1].value]].append(pearson(c, ex))
[1583]352
353    return ttest(distances[0], distances[1])
354
355def mat_ni(data, matcher):
[1774]356    """ Return (in a tuple):
357        - a gene matcher that matches to the attribute names of data
358        - a dictionary attribute names -> indices in the data set.
359    """
[1583]360    nm = matcher([at.name for at in data.domain.attributes])
361    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes))
362    return nm, name_ind
363
364def select_genesets(nm, gene_sets, min_size=3, max_size=1000, min_part=0.1):
365    """ Returns a list of gene sets that have sizes in limits """
366
367    def ok_sizes(gs):
368        """compares sizes of genesets to limitations"""
369        transl = filter(lambda x: x != None, [ nm.umatch(gene) for gene in gs.genes ])
370        if len(transl) >= min_size \
371            and len(transl) <= max_size \
372            and float(len(transl))/len(gs.genes) >= min_part:
373            return True
374        return False
375
376    return filter(ok_sizes, gene_sets) 
[1578]377
[1592]378def vou(ex, gn, indices):
379    """ returns the value or "?" for the given gene name gn"""
380    if gn not in indices:
381        return "?"
382    else:
383        return ex[indices[gn]].value
384
385class SetSig(GeneSetTrans):
386
[1707]387    def __init__(self, **kwargs):
388        self.no_unknowns = kwargs.pop("no_unknowns", False)
[1772]389        self.check_same = kwargs.pop("check_same", False)
[1707]390        super(SetSig, self).__init__(**kwargs)
391
[1594]392    def build_feature(self, data, gs):
[1592]393
[1594]394        at = Orange.feature.Continuous(name=str(gs))
[1707]395        geneset = list(gs.genes)
396        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
397        indices = [ name_ind[gene] for gene in genes ]
[1708]398        takegenes = [ geneset[i] for i in takegenes ]
[1572]399
[1708]400        def t(ex, w, gs=gs, data=data, indices=indices, takegenes=takegenes):
401            nm2, name_ind2, genes2 = self._match_instance(ex, takegenes)
[1593]402
[1707]403            domain = Orange.data.Domain([data.domain.attributes[i] for i in indices], data.domain.class_var)
[1594]404            datao = Orange.data.Table(domain, data)
405           
406            #convert the example to the same domain
407            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
408            example = Orange.data.Instance(domain, exvalues)
[1583]409
[1772]410            return setSig_example_geneset(example, datao, self.no_unknowns, check_same=self.check_same) #only this one is setsig specific
[1594]411     
412        at.get_value_from = t
413        return at
[1593]414
[1599]415class ParametrizedTransformation(GeneSetTrans):
[1595]416
[1599]417    def _get_par(self, datao):
[1690]418        """ Get parameters for a subset of data, that comprises only the gene set """
[1599]419        pass
420       
421    def _use_par(self, ex, constructt):
422        pass
423   
[1595]424    def build_feature(self, data, gs):
425
426        at = Orange.feature.Continuous(name=str(gs))
427
428        geneset = list(gs.genes)
[1598]429        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
[1595]430        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
431        datao = Orange.data.Table(domain, data)
[1709]432        takegenes = [ geneset[i] for i in takegenes ]
[1595]433
[1599]434        constructt = self._get_par(datao)
[1595]435
[1709]436        def t(ex, w, constructt=constructt, takegenes=takegenes, domain=domain):
437            nm2, name_ind2, genes2 = self._match_instance(ex, takegenes)
[1595]438         
439            #convert the example to the same domain
440            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
[1690]441            ex = Orange.data.Instance(domain, exvalues)
[1595]442
[1599]443            return self._use_par(ex, constructt)
[1772]444       
[1595]445        at.get_value_from = t
[1772]446        at.dbg = constructt #for debugging
447       
[1595]448        return at
449
[1599]450class PLS(ParametrizedTransformation):
[1593]451
[1599]452    def _get_par(self, datao):
453        return PLSCall(datao, nc=1, y=[datao.domain.class_var])
454       
455    def _use_par(self, ex, constructt):
[1690]456        ex = [ ex[i].value for i in range(len(ex.domain.attributes)) ]
[1599]457        xmean, W, P, _ = constructt
458        ex = ex - xmean # same input transformation
[1593]459
[1599]460        nc = W.shape[1]
[1593]461
[1599]462        TR = numpy.empty((1, nc))
463        XR = ex
[1593]464
[1599]465        dot = numpy.dot
[1593]466
[1599]467        for i in range(nc):
468           t = dot(XR, W[:,i].T)
469           XR = XR - t*numpy.array([P[:,i]])
470           TR[:,i] = t
[1593]471
[1599]472        return TR[0][0]
473       
474class PCA(ParametrizedTransformation):
[1593]475
[1599]476    def _get_par(self, datao):
477        return pca(datao)
478
479    def _use_par(self, arr, constructt):
[1690]480        arr = [ arr[i].value for i in range(len(arr.domain.attributes)) ]
[1599]481        evals, evect, xmean = constructt
482
483        arr = arr - xmean # same input transformation - a row in a matrix
484        ev0 = evect[0] #this is a row in a matrix - do a dot product
485        a = numpy.dot(arr, ev0)
486
487        return a
[1593]488
[1592]489class SimpleFun(GeneSetTrans):
490
[1594]491    def build_feature(self, data, gs):
[1592]492
[1594]493        at = Orange.feature.Continuous(name=str(gs))
[1592]494
[1594]495        def t(ex, w, gs=gs):
496            geneset = list(gs.genes)
[1598]497            nm2, name_ind2, genes2 = self._match_instance(ex, geneset)
[1594]498           
499            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
500            exvalues = filter(lambda x: x != "?", exvalues)
[1592]501
[1594]502            return self.fn(exvalues)
503     
504        at.get_value_from = t
505        return at
[1592]506
507class Mean(SimpleFun):
508
509    def __init__(self, **kwargs):
510       self.fn = numpy.mean
511       super(Mean, self).__init__(**kwargs)
512
513class Median(SimpleFun):
514
515    def __init__(self, **kwargs):
516       self.fn = numpy.median
517       super(Median, self).__init__(**kwargs)
[1572]518
[1597]519class GSA(GeneSetTrans):
520
521    def build_features(self, data, gene_sets):
522
523        attributes = []
524
525        def tscorec(data, at, cache=None):
526            ma = obiExpression.MA_t_test()(at,data)
527            return ma
528
529        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
530
531        def to_z_score(t):
532            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
533
534        zscores = map(to_z_score, tscores)
535
536        for gs in gene_sets:
537
538            at = Orange.feature.Continuous(name=str(gs))
539
540            geneset = list(gs.genes)
[1598]541            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
[1597]542            #take each gene only once
543            genes = set(genes)
544
545            D = numpy.mean([max(zscores[name_ind[g]],0) for g in genes]) \
546                + numpy.mean([min(zscores[name_ind[g]],0) for g in genes])
547
548            if D >= 0:
549                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] > 0.0 ]
550            else:
551                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] < 0.0 ]
552
553            def t(ex, w, consider_genes=consider_genes):
[1598]554                nm2, name_ind2, genes2 = self._match_instance(ex, consider_genes)
[1597]555             
556                #convert the example to the same domain
557                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
558                exvalues = filter(lambda x: x != "?", exvalues)
559             
560                return numpy.mean(exvalues)
561
562            at.get_value_from = t
563            attributes.append(at)
564
565        return attributes
566
[1600]567def tscorec(data, at, cache=None):
568    """ Cached attribute  tscore calculation """
569    if cache != None and at in cache: return cache[at]
570    ma = obiExpression.MA_t_test()(at,data)
571    if cache != None: cache[at] = ma
572    return ma
573
574def nth(l, n):
575    return [a[n] for a in l]
576
577def compute_corg(data, inds, tscorecache):
578    """
579    Compute CORG for this geneset specified with gene inds
580    in the example table. Output is the list of gene inds
581    in CORG.
582
583    """
584    #order member genes by their t-scores: decreasing, if av(t-score) >= 0,
585    #else increasing
586    tscores = [ tscorec(data, at, tscorecache) for at in inds ]
587    sortedinds = nth(sorted(zip(inds,tscores), key=lambda x: x[1], \
588        reverse=numpy.mean(tscores) >= 0), 0)
589
590    def S(corg):
591        """ Activity score separation - S(G) in
592        the article """
593        asv = Orange.feature.Continuous(name='AS')
594        asv.getValueFrom = lambda ex,rw: Orange.data.Value(asv, corgs_activity_score(ex, corg))
595        data2 = Orange.data.Table(Orange.data.Domain([asv], data.domain.classVar), data)
596        return abs(tscorec(data2, 0)) #FIXME absolute - nothing in the article abs()
597           
598    #greedily find CORGS procing the best separation
599    g = S(sortedinds[:1])
600    bg = 1
601    for a in range(2, len(sortedinds)+1):
602        tg = S(sortedinds[:a])
[1772]603        if tg > g: #improvement
[1600]604            g = tg
605            bg = a
606        else:
607            break
608       
[1772]609    return sortedinds[:bg] #FIXED - one too many was taken
[1600]610
611class CORGs(ParametrizedTransformation):
612    """
613    WARNING: input has to be z_ij table! each gene needs to be normalized
614    (mean=0, stdev=1) for all samples.
615    """
616
617    def __call__(self, *args, **kwargs):
618        self.tscorecache = {} #reset a cache
619        return super(CORGs, self).__call__(*args, **kwargs)
620
621    def build_feature(self, data, gs):
622
623        at = Orange.feature.Continuous(name=str(gs))
624        geneset = list(gs.genes)
625
626        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
627        indices = compute_corg(data, [ name_ind[g] for g in genes ], self.tscorecache)
628
629        ind_names = dict( (a,b) for b,a in name_ind.items() )
630        selected_genes = sorted(set([to_geneset[ind_names[i]] for i in indices]))
631           
632        def t(ex, w, corg=selected_genes): #copy od the data
633            nm2, name_ind2, genes2 = self._match_instance(ex, corg, None)
634            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ]
635            return sum(v if v != '?' else 0.0 for v in exvalues)/len(corg)**0.5
636     
637        at.get_value_from = t
638        return at
639
[1774]640def compute_llr(data, inds, cache):
641    """
642    Compute CORG for this geneset specified with gene inds
643    in the example table. Output is the list of gene inds
644    in CORG.
645    """
646    def gaussianc(data, at, cache=None):
647        """ Cached attribute  tscore calculation """
648        if cache != None and at in cache: return cache[at]
649        ma = estimate_gaussian_per_class(data, at)
650        if cache != None: cache[at] = ma
651        return ma
652
653    gf = [ gaussianc(data, at, cache) for at in inds ]
654    return gf
655
656""" To avoid scipy overhead """
657from math import pi
658_norm_pdf_C = math.sqrt(2*pi)
659_norm_pdf_logC = math.log(_norm_pdf_C)
660def _norm_logpdf(x, mi, std):
661    return -(x-mi)**2 / (2.0*std**2) - _norm_pdf_logC - math.log(std)
662
663def _llrlogratio(v, mi1, std1, mi2, std2):
664    if mi1 == None or std1 == None or mi2 == None or std2 == None:
665        return 0. #problem with estimation
666    #lpdf1 = scipy.stats.norm.logpdf(v, mi1, std1)
667    #lpdf2 = scipy.stats.norm.logpdf(v, mi2, std2)
668    lpdf1 = _norm_logpdf(v, mi1, std1) #avoids scipy's overhead
669    lpdf2 = _norm_logpdf(v, mi2, std2)
670    r = lpdf1 - lpdf2
671    return r
672
673class LLR(ParametrizedTransformation):
674    """
675    From Su et al: Accurate and Reliable Cancer Classification Base
676    on Probalistic Inference of Pathway Activity. Plos One, 2009.
677    """
678
679    def __init__(self, **kwargs):
680        self.normalize = kwargs.pop("normalize", False) #normalize final results
681        super(LLR, self).__init__(**kwargs)
682
683    def __call__(self, *args, **kwargs):
684        self._gauss_cache = {} #caching of gaussian estimates
685        self._normalizec = {}
686        return super(LLR, self).__call__(*args, **kwargs)
687
688    def build_feature(self, data, gs):
689
690        at = Orange.feature.Continuous(name=str(gs))
691        geneset = list(gs.genes)
692
693        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
694
695        gsi = [ name_ind[g] for g in genes ]
696        gausse = compute_llr(data, gsi, self._gauss_cache)
697        genes_gs = [ to_geneset[g] for g in genes ]
698
699        if self.normalize: # per (3) in the paper
700            #compute log ratios for all samples and genes from this gene set
701            for i, gene_gs, g in zip(gsi, genes_gs, gausse):
702                if gene_gs not in self._normalizec: #skip if computed already
703                    r = [ _llrlogratio(ex[i].value, *g) for ex in data ]
704                    self._normalizec[gene_gs] = (statc.mean(r), statc.std(r))
705
706        def t(ex, w, genes_gs=genes_gs, gausse=gausse, normalizec=self._normalizec):
707            nm2, name_ind2, genes2 = self._match_instance(ex, genes_gs, None)
708            gsvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ]
709
710            vals = [ _llrlogratio(v, *g) if v != "?" else 0.0 for v,g in zip(gsvalues, gausse) ]
711
712            if len(normalizec): #normalize according to (3)
713                vals2 = []
714                for v,g in zip(vals, genes_gs):
715                    m,s = self._normalizec[g]
716                    vals2.append((v-m)/s)
717                vals = vals2
718           
719            return sum(vals)
720
721     
722        at.get_value_from = t
723        return at
724
[1572]725if __name__ == "__main__":
726
727    data = Orange.data.Table("iris")
728    gsets = obiGeneSets.collections({
[1597]729        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
[1572]730        "f3": ['sepal length', 'sepal width', 'petal length'],
731        "l3": ['sepal width', 'petal length', 'petal width'],
732        })
[1587]733    matcher = obiGene.matcher([])
734    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
[1572]735
[1587]736    """
737    data = Orange.data.Table("DLBCL_200a")
738    gsets = obiGeneSets.collections((("KEGG",),"9606"))
739    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
740    choosen_cv = None
741    """
[1572]742
743    def to_old_dic(d, data):
744        ar = defaultdict(list)
745        for ex1 in data:
746            ex = d(ex1)
747            for a,x in zip(d.attributes, ex):
748                ar[a.name].append(x.value)
749        return ar
750
751    def pp2(ar):
752        ol =  sorted(ar.items())
753        print '\n'.join([ a + ": " +str(b) for a,b in ol])
754
[1607]755    ass = Assess(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
[1572]756    ar = to_old_dic(ass.domain, data[:5])
757    pp2(ar)
Note: See TracBrowser for help on using the repository browser.