source: orange-bioinformatics/obiGeneSetSig.py @ 1601:d3a65a71731e

Revision 1601:d3a65a71731e, 13.0 KB checked in by markotoplak, 2 years ago (diff)

Fixes supporting the move of (most of) Orange.misc to Orange.utils.

Line 
1import Orange
2import Orange.utils
3import obiGeneSets
4import obiGene
5import numpy
6from collections import defaultdict
7import stats
8from obiGsea import takeClasses
9from obiAssess import pca, PLSCall, corgs_activity_score
10import obiExpression
11import scipy.stats
12
13#STILL MISSING: Assess
14
15def setSig_example_geneset(ex, data):
16    """ Gets learning data and example with the same domain, both
17    containing only genes from the gene set. """
18
19    distances = [ [], [] ]   
20
21    def pearsonr(v1, v2):
22        return numpy.corrcoef([v1, v2])[0,1]
23
24    def pearson(ex1, ex2):
25        #leaves undefined elements out
26
27        attrs = range(len(ex1.domain.attributes))
28        vals1 = [ ex1[i].value for i in attrs ]
29        vals2 = [ ex2[i].value for i in attrs ]
30
31        common = [ True if v1 != "?" and v2 != "?" else False \
32            for v1,v2 in zip(vals1,vals2) ]
33        vals1 = [ v for v,c in zip(vals1, common) if c ]
34        vals2 = [ v for v,c in zip(vals2, common) if c ]
35
36        return numpy.corrcoef([vals1, vals2])[0,1]
37
38    def ttest(ex1, ex2):
39        try:
40            return stats.lttest_ind(ex1, ex2)[0]
41        except:
42            return 0.0
43   
44    #maps class value to its index
45    classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.class_var.values) ])
46 
47    #create distances to all learning data - save or other class
48    for c in data:
49        distances[classValueMap[c[-1].value]].append(pearson(c, ex))
50
51    return ttest(distances[0], distances[1])
52
53def mat_ni(data, matcher):
54    nm = matcher([at.name for at in data.domain.attributes])
55    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes))
56    return nm, name_ind
57
58def select_genesets(nm, gene_sets, min_size=3, max_size=1000, min_part=0.1):
59    """ Returns a list of gene sets that have sizes in limits """
60
61    def ok_sizes(gs):
62        """compares sizes of genesets to limitations"""
63        transl = filter(lambda x: x != None, [ nm.umatch(gene) for gene in gs.genes ])
64        if len(transl) >= min_size \
65            and len(transl) <= max_size \
66            and float(len(transl))/len(gs.genes) >= min_part:
67            return True
68        return False
69
70    return filter(ok_sizes, gene_sets) 
71
72class GeneSetTrans(object):
73
74    __new__ = Orange.utils._orange__new__(object)
75
76    def _mat_ni(self, data):
77        """ With cached gene matchers. """
78        if data.domain not in self._cache:
79            self._cache[data.domain] = mat_ni(data, self.matcher)
80        return self._cache[data.domain]
81
82    def _match_instance(self, instance, geneset, takegenes=None):
83        nm, name_ind = self._mat_ni(instance)
84        genes = [ nm.umatch(gene) for gene in geneset ]
85        if takegenes:
86            genes = [ genes[i] for i in takegenes ]
87        return nm, name_ind, genes
88
89    def _match_data(self, data, geneset, odic=False):
90        nm, name_ind = self._mat_ni(data)
91        genes = [ nm.umatch(gene) for gene in geneset ]
92        if odic:
93            to_geneset = dict(zip(genes, geneset))
94        takegenes = [ i for i,a in enumerate(genes) if a != None ]
95        genes = [ genes[i] for i in takegenes ]
96        if odic:
97            return nm, name_ind, genes, takegenes, to_geneset
98        else:
99            return nm, name_ind, genes, takegenes
100
101    def __init__(self, matcher=None, gene_sets=None, min_size=3, max_size=1000, min_part=0.1, class_values=None):
102        self.matcher = matcher
103        self.gene_sets = gene_sets
104        self.min_size = min_size
105        self.max_size = max_size
106        self.min_part = min_part
107        self.class_values = class_values
108        self._cache = {}
109
110    def __call__(self, data, weight_id=None):
111
112        #selection of classes and gene sets
113        data = takeClasses(data, classValues=self.class_values)
114        nm,_ =  self._mat_ni(data)
115        gene_sets = select_genesets(nm, self.gene_sets, self.min_size, self.max_size, self.min_part)
116
117        #build a new domain
118        newfeatures = self.build_features(data, gene_sets)
119        newdomain = Orange.data.Domain(newfeatures, data.domain.class_var)
120        return Orange.data.Table(newdomain, data)
121
122    def build_features(self, data, gene_sets):
123        return [ self.build_feature(data, gs) for gs in gene_sets ]
124
125def vou(ex, gn, indices):
126    """ returns the value or "?" for the given gene name gn"""
127    if gn not in indices:
128        return "?"
129    else:
130        return ex[indices[gn]].value
131
132class SetSig(GeneSetTrans):
133
134    def build_feature(self, data, gs):
135
136        at = Orange.feature.Continuous(name=str(gs))
137
138        def t(ex, w, gs=gs, data=data): #copy od the data
139            geneset = list(gs.genes)
140
141            nm, name_ind, genes, takegenes = self._match_data(data, geneset)
142            nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
143
144            domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
145            datao = Orange.data.Table(domain, data)
146           
147            #convert the example to the same domain
148            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
149            example = Orange.data.Instance(domain, exvalues)
150
151            return setSig_example_geneset(example, datao) #only this one is setsig specific
152     
153        at.get_value_from = t
154        return at
155
156class ParametrizedTransformation(GeneSetTrans):
157
158    def _get_par(self, datao):
159        pass
160       
161    def _use_par(self, ex, constructt):
162        pass
163   
164    def build_feature(self, data, gs):
165
166        at = Orange.feature.Continuous(name=str(gs))
167
168        geneset = list(gs.genes)
169        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
170        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
171        datao = Orange.data.Table(domain, data)
172
173        constructt = self._get_par(datao)
174
175        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
176            nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
177         
178            #convert the example to the same domain
179            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
180            ex = numpy.array(exvalues[:-1])
181
182            return self._use_par(ex, constructt)
183           
184        at.get_value_from = t
185        return at
186
187class PLS(ParametrizedTransformation):
188
189    def _get_par(self, datao):
190        return PLSCall(datao, nc=1, y=[datao.domain.class_var])
191       
192    def _use_par(self, ex, constructt):
193        xmean, W, P, _ = constructt
194        ex = ex - xmean # same input transformation
195
196        nc = W.shape[1]
197
198        TR = numpy.empty((1, nc))
199        XR = ex
200
201        dot = numpy.dot
202
203        for i in range(nc):
204           t = dot(XR, W[:,i].T)
205           XR = XR - t*numpy.array([P[:,i]])
206           TR[:,i] = t
207
208        return TR[0][0]
209       
210class PCA(ParametrizedTransformation):
211
212    def _get_par(self, datao):
213        return pca(datao)
214
215    def _use_par(self, arr, constructt):
216        evals, evect, xmean = constructt
217
218        arr = arr - xmean # same input transformation - a row in a matrix
219        ev0 = evect[0] #this is a row in a matrix - do a dot product
220        a = numpy.dot(arr, ev0)
221
222        return a
223
224class SimpleFun(GeneSetTrans):
225
226    def build_feature(self, data, gs):
227
228        at = Orange.feature.Continuous(name=str(gs))
229
230        def t(ex, w, gs=gs):
231            geneset = list(gs.genes)
232            nm2, name_ind2, genes2 = self._match_instance(ex, geneset)
233           
234            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
235            exvalues = filter(lambda x: x != "?", exvalues)
236
237            return self.fn(exvalues)
238     
239        at.get_value_from = t
240        return at
241
242class Mean(SimpleFun):
243
244    def __init__(self, **kwargs):
245       self.fn = numpy.mean
246       super(Mean, self).__init__(**kwargs)
247
248class Median(SimpleFun):
249
250    def __init__(self, **kwargs):
251       self.fn = numpy.median
252       super(Median, self).__init__(**kwargs)
253
254class GSA(GeneSetTrans):
255
256    def build_features(self, data, gene_sets):
257
258        attributes = []
259
260        def tscorec(data, at, cache=None):
261            ma = obiExpression.MA_t_test()(at,data)
262            return ma
263
264        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
265
266        def to_z_score(t):
267            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
268
269        zscores = map(to_z_score, tscores)
270
271        for gs in gene_sets:
272
273            at = Orange.feature.Continuous(name=str(gs))
274
275            geneset = list(gs.genes)
276            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
277            #take each gene only once
278            genes = set(genes)
279
280            D = numpy.mean([max(zscores[name_ind[g]],0) for g in genes]) \
281                + numpy.mean([min(zscores[name_ind[g]],0) for g in genes])
282
283            if D >= 0:
284                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] > 0.0 ]
285            else:
286                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] < 0.0 ]
287
288            def t(ex, w, consider_genes=consider_genes):
289                nm2, name_ind2, genes2 = self._match_instance(ex, consider_genes)
290             
291                #convert the example to the same domain
292                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
293                exvalues = filter(lambda x: x != "?", exvalues)
294             
295                return numpy.mean(exvalues)
296
297            at.get_value_from = t
298            attributes.append(at)
299
300        return attributes
301
302def tscorec(data, at, cache=None):
303    """ Cached attribute  tscore calculation """
304    if cache != None and at in cache: return cache[at]
305    ma = obiExpression.MA_t_test()(at,data)
306    if cache != None: cache[at] = ma
307    return ma
308
309def nth(l, n):
310    return [a[n] for a in l]
311
312def compute_corg(data, inds, tscorecache):
313    """
314    Compute CORG for this geneset specified with gene inds
315    in the example table. Output is the list of gene inds
316    in CORG.
317
318    """
319    #order member genes by their t-scores: decreasing, if av(t-score) >= 0,
320    #else increasing
321    tscores = [ tscorec(data, at, tscorecache) for at in inds ]
322    sortedinds = nth(sorted(zip(inds,tscores), key=lambda x: x[1], \
323        reverse=numpy.mean(tscores) >= 0), 0)
324
325    def S(corg):
326        """ Activity score separation - S(G) in
327        the article """
328        asv = Orange.feature.Continuous(name='AS')
329        asv.getValueFrom = lambda ex,rw: Orange.data.Value(asv, corgs_activity_score(ex, corg))
330        data2 = Orange.data.Table(Orange.data.Domain([asv], data.domain.classVar), data)
331        return abs(tscorec(data2, 0)) #FIXME absolute - nothing in the article abs()
332           
333    #greedily find CORGS procing the best separation
334    g = S(sortedinds[:1])
335    bg = 1
336    for a in range(2, len(sortedinds)+1):
337        tg = S(sortedinds[:a])
338        if tg > g:
339            g = tg
340            bg = a
341        else:
342            break
343       
344    return sortedinds[:a]
345
346class CORGs(ParametrizedTransformation):
347    """
348    WARNING: input has to be z_ij table! each gene needs to be normalized
349    (mean=0, stdev=1) for all samples.
350    """
351
352    def __call__(self, *args, **kwargs):
353        self.tscorecache = {} #reset a cache
354        return super(CORGs, self).__call__(*args, **kwargs)
355
356    def build_feature(self, data, gs):
357
358        at = Orange.feature.Continuous(name=str(gs))
359        geneset = list(gs.genes)
360
361        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
362        indices = compute_corg(data, [ name_ind[g] for g in genes ], self.tscorecache)
363
364        ind_names = dict( (a,b) for b,a in name_ind.items() )
365        selected_genes = sorted(set([to_geneset[ind_names[i]] for i in indices]))
366           
367        def t(ex, w, corg=selected_genes): #copy od the data
368            nm2, name_ind2, genes2 = self._match_instance(ex, corg, None)
369            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ]
370            return sum(v if v != '?' else 0.0 for v in exvalues)/len(corg)**0.5
371     
372        at.get_value_from = t
373        return at
374
375if __name__ == "__main__":
376
377    data = Orange.data.Table("iris")
378    gsets = obiGeneSets.collections({
379        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
380        "f3": ['sepal length', 'sepal width', 'petal length'],
381        "l3": ['sepal width', 'petal length', 'petal width'],
382        })
383    matcher = obiGene.matcher([])
384    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
385
386    """
387    data = Orange.data.Table("DLBCL_200a")
388    gsets = obiGeneSets.collections((("KEGG",),"9606"))
389    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
390    choosen_cv = None
391    """
392
393    def to_old_dic(d, data):
394        ar = defaultdict(list)
395        for ex1 in data:
396            ex = d(ex1)
397            for a,x in zip(d.attributes, ex):
398                ar[a.name].append(x.value)
399        return ar
400
401    def pp2(ar):
402        ol =  sorted(ar.items())
403        print '\n'.join([ a + ": " +str(b) for a,b in ol])
404
405    ass = CORGs(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
406    ar = to_old_dic(ass.domain, data[:5])
407    pp2(ar)
Note: See TracBrowser for help on using the repository browser.