source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1775:023037d1e29a

Revision 1775:023037d1e29a, 25.3 KB checked in by markotoplak, 11 months ago (diff)

obiGeneSetSig: LLR uses combined standard deviation if a standard deviation for a class was zero (this can happen as an artifict of bootstrap sampling).

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