source: orange-bioinformatics/orangecontrib/bio/obiGeneSetSig.py @ 1873:0810c5708cc5

Revision 1873:0810c5708cc5, 30.5 KB checked in by Ales Erjavec <ales.erjavec@…>, 6 months ago (diff)

Moved '_bioinformatics' into orangecontrib namespace.

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