source: orange-bioinformatics/_bioinformatics/obiGeneSetSig.py @ 1802:959e64f489d2

Revision 1802:959e64f489d2, 29.5 KB checked in by markotoplak, 11 months ago (diff)

obiGeneSetSig: Prevent crashing of LLR if the standard deviation is zero.

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