source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1786:df95834ddac7

Revision 1786:df95834ddac7, 29.4 KB checked in by markotoplak, 11 months ago (diff)

obiGeneSetSig: SPCA with threshold selection by permutation.

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