source: orange-bioinformatics/obiGeneSetSig.py @ 1598:fda2dc9e724e

Revision 1598:fda2dc9e724e, 10.9 KB checked in by markotoplak, 2 years ago (diff)

Some refactoring of 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 _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 PLS(GeneSetTrans):
157
158    def build_feature(self, data, gs):
159
160        at = Orange.feature.Continuous(name=str(gs))
161
162        geneset = list(gs.genes)
163        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
164        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
165        datao = Orange.data.Table(domain, data)
166
167        constructt = PLSCall(datao, nc=1, y=[datao.domain.class_var])
168
169        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
170            nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
171         
172            #convert the example to the same domain
173            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
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        nm, name_ind, genes, takegenes = self._match_data(data, geneset)
204        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
205        datao = Orange.data.Table(domain, data)
206
207        constructt = pca(datao)
208
209        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
210            nm2, name_ind2, genes2 = self._match_instance(ex, geneset, takegenes)
211         
212            #convert the example to the same domain
213            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
214           
215            arr = numpy.array(exvalues[:-1])
216           
217            evals, evect, xmean = constructt
218
219            arr = arr - xmean # same input transformation - a row in a matrix
220            ev0 = evect[0] #this is a row in a matrix - do a dot product
221            a = numpy.dot(arr, ev0)
222
223            return a
224                 
225        at.get_value_from = t
226        return at
227
228class SimpleFun(GeneSetTrans):
229
230    def build_feature(self, data, gs):
231
232        at = Orange.feature.Continuous(name=str(gs))
233
234        def t(ex, w, gs=gs):
235            geneset = list(gs.genes)
236            nm2, name_ind2, genes2 = self._match_instance(ex, geneset)
237           
238            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
239            exvalues = filter(lambda x: x != "?", exvalues)
240
241            return self.fn(exvalues)
242     
243        at.get_value_from = t
244        return at
245
246class Mean(SimpleFun):
247
248    def __init__(self, **kwargs):
249       self.fn = numpy.mean
250       super(Mean, self).__init__(**kwargs)
251
252class Median(SimpleFun):
253
254    def __init__(self, **kwargs):
255       self.fn = numpy.median
256       super(Median, self).__init__(**kwargs)
257
258class GSA(GeneSetTrans):
259
260    def build_features(self, data, gene_sets):
261
262        attributes = []
263
264        def tscorec(data, at, cache=None):
265            ma = obiExpression.MA_t_test()(at,data)
266            return ma
267
268        tscores = [ tscorec(data, at) for at in data.domain.attributes ]
269
270        def to_z_score(t):
271            return float(scipy.stats.norm.ppf(scipy.stats.t.cdf(t, len(data)-2)))
272
273        zscores = map(to_z_score, tscores)
274
275        for gs in gene_sets:
276
277            at = Orange.feature.Continuous(name=str(gs))
278
279            geneset = list(gs.genes)
280            nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True)
281            #take each gene only once
282            genes = set(genes)
283
284            D = numpy.mean([max(zscores[name_ind[g]],0) for g in genes]) \
285                + numpy.mean([min(zscores[name_ind[g]],0) for g in genes])
286
287            if D >= 0:
288                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] > 0.0 ]
289            else:
290                consider_genes = [ to_geneset[g] for g in genes if zscores[name_ind[g]] < 0.0 ]
291
292            def t(ex, w, consider_genes=consider_genes):
293                nm2, name_ind2, genes2 = self._match_instance(ex, consider_genes)
294             
295                #convert the example to the same domain
296                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
297                exvalues = filter(lambda x: x != "?", exvalues)
298             
299                return numpy.mean(exvalues)
300
301            at.get_value_from = t
302            attributes.append(at)
303
304        return attributes
305
306if __name__ == "__main__":
307
308    data = Orange.data.Table("iris")
309    gsets = obiGeneSets.collections({
310        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
311        "f3": ['sepal length', 'sepal width', 'petal length'],
312        "l3": ['sepal width', 'petal length', 'petal width'],
313        })
314    matcher = obiGene.matcher([])
315    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
316
317    """
318    data = Orange.data.Table("DLBCL_200a")
319    gsets = obiGeneSets.collections((("KEGG",),"9606"))
320    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
321    choosen_cv = None
322    """
323
324    def to_old_dic(d, data):
325        ar = defaultdict(list)
326        for ex1 in data:
327            ex = d(ex1)
328            for a,x in zip(d.attributes, ex):
329                ar[a.name].append(x.value)
330        return ar
331
332    def pp2(ar):
333        ol =  sorted(ar.items())
334        print '\n'.join([ a + ": " +str(b) for a,b in ol])
335
336    ass = GSA(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
337    ar = to_old_dic(ass.domain, data[:5])
338    pp2(ar)
Note: See TracBrowser for help on using the repository browser.