source: orange-bioinformatics/obiGeneSetSig.py @ 1599:d26971fdd077

Revision 1599:d26971fdd077, 10.5 KB checked in by markotoplak, 2 years ago (diff)

obiGeneSetSig: refactoring.

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 _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
302if __name__ == "__main__":
303
304    data = Orange.data.Table("iris")
305    gsets = obiGeneSets.collections({
306        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
307        "f3": ['sepal length', 'sepal width', 'petal length'],
308        "l3": ['sepal width', 'petal length', 'petal width'],
309        })
310    matcher = obiGene.matcher([])
311    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
312
313    """
314    data = Orange.data.Table("DLBCL_200a")
315    gsets = obiGeneSets.collections((("KEGG",),"9606"))
316    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
317    choosen_cv = None
318    """
319
320    def to_old_dic(d, data):
321        ar = defaultdict(list)
322        for ex1 in data:
323            ex = d(ex1)
324            for a,x in zip(d.attributes, ex):
325                ar[a.name].append(x.value)
326        return ar
327
328    def pp2(ar):
329        ol =  sorted(ar.items())
330        print '\n'.join([ a + ": " +str(b) for a,b in ol])
331
332    ass = PCA(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
333    ar = to_old_dic(ass.domain, data[:5])
334    pp2(ar)
Note: See TracBrowser for help on using the repository browser.