source: orange-bioinformatics/obiGeneSetSig.py @ 1597:574fe4f8d1d3

Revision 1597:574fe4f8d1d3, 11.3 KB checked in by markotoplak, 2 years ago (diff)

Added GSA to obiGeneSetSig.

Line 
1import Orange
2import Orange.misc
3import obiGeneSets
4import obiGene
5import numpy
6from collections import defaultdict
7import stats
8from obiGsea import takeClasses
9from obiAssess import pca, PLSCall
10import obiExpression
11import scipy.stats
12
13#STILL MISSING: Assess, CORGs
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.misc._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 __init__(self, matcher=None, gene_sets=None, min_size=3, max_size=1000, min_part=0.1, class_values=None):
83        self.matcher = matcher
84        self.gene_sets = gene_sets
85        self.min_size = min_size
86        self.max_size = max_size
87        self.min_part = min_part
88        self.class_values = class_values
89        self._cache = {}
90
91    def __call__(self, data, weight_id=None):
92
93        #selection of classes and gene sets
94        data = takeClasses(data, classValues=self.class_values)
95        nm,_ =  self._mat_ni(data)
96        gene_sets = select_genesets(nm, self.gene_sets, self.min_size, self.max_size, self.min_part)
97
98        #build a new domain
99        newfeatures = self.build_features(data, gene_sets)
100        newdomain = Orange.data.Domain(newfeatures, data.domain.class_var)
101        return Orange.data.Table(newdomain, data)
102
103    def build_features(self, data, gene_sets):
104        return [ self.build_feature(data, gs) for gs in gene_sets ]
105
106def vou(ex, gn, indices):
107    """ returns the value or "?" for the given gene name gn"""
108    if gn not in indices:
109        return "?"
110    else:
111        return ex[indices[gn]].value
112
113class SetSig(GeneSetTrans):
114
115    def build_feature(self, data, gs):
116
117        at = Orange.feature.Continuous(name=str(gs))
118
119        def t(ex, w, gs=gs, data=data): #copy od the data
120            geneset = list(gs.genes)
121
122            nm, name_ind = self._mat_ni(data)
123            nm2, name_ind2 = self._mat_ni(ex)
124
125            genes = [ nm.umatch(gene) for gene in geneset ]
126            genes2 = [ nm2.umatch(gene) for gene in geneset ]
127
128            takegenes = [ i for i,a in enumerate(genes) if a != None ]
129
130            genes = [ genes[i] for i in takegenes ]
131            genes2 = [ genes2[i] for i in takegenes ]
132
133            domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
134            datao = Orange.data.Table(domain, data)
135           
136            #convert the example to the same domain
137            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
138            example = Orange.data.Instance(domain, exvalues)
139
140            return setSig_example_geneset(example, datao) #only this one is setsig specific
141     
142        at.get_value_from = t
143        return at
144
145class PLS(GeneSetTrans):
146
147    def build_feature(self, data, gs):
148
149        at = Orange.feature.Continuous(name=str(gs))
150
151        geneset = list(gs.genes)
152
153        nm, name_ind = self._mat_ni(data)
154        genes = [ nm.umatch(gene) for gene in geneset ]
155        takegenes = [ i for i,a in enumerate(genes) if a != None ]
156        genes = [ genes[i] for i in takegenes ]
157
158        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
159
160        datao = Orange.data.Table(domain, data)
161
162        xmean, W, P, T = PLSCall(datao, nc=1, y=[datao.domain.class_var])
163        constructt = xmean, W, P
164
165        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
166
167            nm2, name_ind2 = self._mat_ni(ex)
168            genes2 = [ nm2.umatch(gene) for gene in geneset ]
169            genes2 = [ genes2[i] for i in takegenes ]
170         
171            #convert the example to the same domain
172            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
173           
174            ex = numpy.array(exvalues[:-1])
175
176            xmean, W, P = constructt
177            ex = ex - xmean # same input transformation
178
179            nc = W.shape[1]
180
181            TR = numpy.empty((1, nc))
182            XR = ex
183
184            dot = numpy.dot
185
186            for i in range(nc):
187               t = dot(XR, W[:,i].T)
188               XR = XR - t*numpy.array([P[:,i]])
189               TR[:,i] = t
190
191            return TR[0][0]
192
193        at.get_value_from = t
194        return at
195
196class PCA(GeneSetTrans):
197
198    def build_feature(self, data, gs):
199
200        at = Orange.feature.Continuous(name=str(gs))
201
202        geneset = list(gs.genes)
203
204        nm, name_ind = self._mat_ni(data)
205        genes = [ nm.umatch(gene) for gene in geneset ]
206        takegenes = [ i for i,a in enumerate(genes) if a != None ]
207        genes = [ genes[i] for i in takegenes ]
208
209        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
210
211        datao = Orange.data.Table(domain, data)
212
213        evals, evect, xmean = pca(datao)
214        constructt = evals, evect, xmean
215
216        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
217
218            nm2, name_ind2 = self._mat_ni(ex)
219            genes2 = [ nm2.umatch(gene) for gene in geneset ]
220            genes2 = [ genes2[i] for i in takegenes ]
221         
222            #convert the example to the same domain
223            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
224           
225            arr = numpy.array(exvalues[:-1])
226           
227            evals, evect, xmean = constructt
228
229            arr = arr - xmean # same input transformation - a row in a matrix
230            ev0 = evect[0] #this is a row in a matrix - do a dot product
231            a = numpy.dot(arr, ev0)
232
233            return a
234                 
235        at.get_value_from = t
236        return at
237
238class SimpleFun(GeneSetTrans):
239
240    def build_feature(self, data, gs):
241
242        at = Orange.feature.Continuous(name=str(gs))
243
244        def t(ex, w, gs=gs):
245            geneset = list(gs.genes)
246            nm2, name_ind2 = self._mat_ni(ex)
247            genes2 = [ nm2.umatch(gene) for gene in geneset ]
248           
249            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
250            exvalues = filter(lambda x: x != "?", exvalues)
251
252            return self.fn(exvalues)
253     
254        at.get_value_from = t
255        return at
256
257class Mean(SimpleFun):
258
259    def __init__(self, **kwargs):
260       self.fn = numpy.mean
261       super(Mean, self).__init__(**kwargs)
262
263class Median(SimpleFun):
264
265    def __init__(self, **kwargs):
266       self.fn = numpy.median
267       super(Median, self).__init__(**kwargs)
268
269class GSA(GeneSetTrans):
270
271    def build_features(self, data, gene_sets):
272
273        attributes = []
274
275        def tscorec(data, at, cache=None):
276            ma = obiExpression.MA_t_test()(at,data)
277            return ma
278
279        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
280
281        def to_z_score(t):
282            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
283
284        zscores = map(to_z_score, tscores)
285
286        nm, name_ind = self._mat_ni(data)
287
288        for gs in gene_sets:
289
290            at = Orange.feature.Continuous(name=str(gs))
291
292            geneset = list(gs.genes)
293
294            genes = [ nm.umatch(gene) for gene in geneset ]
295
296            to_geneset = dict(zip(genes, geneset))
297
298            takegenes = [ i for i,a in enumerate(genes) if a != None ]
299            genes = [ genes[i] for i in takegenes ]
300
301            #take each gene only once
302            genes = set(genes)
303
304            D = numpy.mean([max(zscores[name_ind[g]],0) for g in genes]) \
305                + numpy.mean([min(zscores[name_ind[g]],0) for g in genes])
306
307            if D >= 0:
308                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] > 0.0 ]
309            else:
310                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] < 0.0 ]
311
312            def t(ex, w, consider_genes=consider_genes):
313                #consider_genes included genes from the gene set that
314                #should be combined
315                nm2, name_ind2 = self._mat_ni(ex)
316                genes2 = [ nm2.umatch(gene) for gene in consider_genes ]
317             
318                #convert the example to the same domain
319                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
320                exvalues = filter(lambda x: x != "?", exvalues)
321             
322                return numpy.mean(exvalues)
323
324            at.get_value_from = t
325            attributes.append(at)
326
327        return attributes
328
329if __name__ == "__main__":
330
331    data = Orange.data.Table("iris")
332    gsets = obiGeneSets.collections({
333        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
334        "f3": ['sepal length', 'sepal width', 'petal length'],
335        "l3": ['sepal width', 'petal length', 'petal width'],
336        })
337    matcher = obiGene.matcher([])
338    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
339
340    """
341    data = Orange.data.Table("DLBCL_200a")
342    gsets = obiGeneSets.collections((("KEGG",),"9606"))
343    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
344    choosen_cv = None
345    """
346
347    def to_old_dic(d, data):
348        ar = defaultdict(list)
349        for ex1 in data:
350            ex = d(ex1)
351            for a,x in zip(d.attributes, ex):
352                ar[a.name].append(x.value)
353        return ar
354
355    def pp2(ar):
356        ol =  sorted(ar.items())
357        print '\n'.join([ a + ": " +str(b) for a,b in ol])
358
359    ass = GSA(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
360    ar = to_old_dic(ass.domain, data[:5])
361    pp2(ar)
Note: See TracBrowser for help on using the repository browser.