source: orange-bioinformatics/obiGeneSetSig.py @ 1595:8ef1b76b2e6d

Revision 1595:8ef1b76b2e6d, 9.3 KB checked in by markotoplak, 2 years ago (diff)

Added PLS 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
10
11def setSig_example_geneset(ex, data):
12    """ Gets learning data and example with the same domain, both
13    containing only genes from the gene set. """
14
15    distances = [ [], [] ]   
16
17    def pearsonr(v1, v2):
18        return numpy.corrcoef([v1, v2])[0,1]
19
20    def pearson(ex1, ex2):
21        #leaves undefined elements out
22
23        attrs = range(len(ex1.domain.attributes))
24        vals1 = [ ex1[i].value for i in attrs ]
25        vals2 = [ ex2[i].value for i in attrs ]
26
27        common = [ True if v1 != "?" and v2 != "?" else False \
28            for v1,v2 in zip(vals1,vals2) ]
29        vals1 = [ v for v,c in zip(vals1, common) if c ]
30        vals2 = [ v for v,c in zip(vals2, common) if c ]
31
32        return numpy.corrcoef([vals1, vals2])[0,1]
33
34    def ttest(ex1, ex2):
35        try:
36            return stats.lttest_ind(ex1, ex2)[0]
37        except:
38            return 0.0
39   
40    #maps class value to its index
41    classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.class_var.values) ])
42 
43    #create distances to all learning data - save or other class
44    for c in data:
45        distances[classValueMap[c[-1].value]].append(pearson(c, ex))
46
47    return ttest(distances[0], distances[1])
48
49def mat_ni(data, matcher):
50    nm = matcher([at.name for at in data.domain.attributes])
51    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes))
52    return nm, name_ind
53
54def select_genesets(nm, gene_sets, min_size=3, max_size=1000, min_part=0.1):
55    """ Returns a list of gene sets that have sizes in limits """
56
57    def ok_sizes(gs):
58        """compares sizes of genesets to limitations"""
59        transl = filter(lambda x: x != None, [ nm.umatch(gene) for gene in gs.genes ])
60        if len(transl) >= min_size \
61            and len(transl) <= max_size \
62            and float(len(transl))/len(gs.genes) >= min_part:
63            return True
64        return False
65
66    return filter(ok_sizes, gene_sets) 
67
68class GeneSetTrans(object):
69
70    __new__ = Orange.misc._orange__new__(object)
71
72    def _mat_ni(self, data):
73        """ With cached gene matchers. """
74        if data.domain not in self._cache:
75            self._cache[data.domain] = mat_ni(data, self.matcher)
76        return self._cache[data.domain]
77
78    def __init__(self, matcher=None, gene_sets=None, min_size=3, max_size=1000, min_part=0.1, class_values=None):
79        self.matcher = matcher
80        self.gene_sets = gene_sets
81        self.min_size = min_size
82        self.max_size = max_size
83        self.min_part = min_part
84        self.class_values = class_values
85        self._cache = {}
86
87    def __call__(self, data, weight_id=None):
88
89        #selection of classes and gene sets
90        data = takeClasses(data, classValues=self.class_values)
91        nm,_ =  self._mat_ni(data)
92        gene_sets = select_genesets(nm, self.gene_sets, self.min_size, self.max_size, self.min_part)
93
94        #build a new domain
95        newfeatures = self.build_features(data, gene_sets)
96        newdomain = Orange.data.Domain(newfeatures, data.domain.class_var)
97        return Orange.data.Table(newdomain, data)
98
99    def build_features(self, data, gene_sets):
100        return [ self.build_feature(data, gs) for gs in gene_sets ]
101
102def vou(ex, gn, indices):
103    """ returns the value or "?" for the given gene name gn"""
104    if gn not in indices:
105        return "?"
106    else:
107        return ex[indices[gn]].value
108
109class SetSig(GeneSetTrans):
110
111    def build_feature(self, data, gs):
112
113        at = Orange.feature.Continuous(name=str(gs))
114
115        def t(ex, w, gs=gs, data=data): #copy od the data
116            geneset = list(gs.genes)
117
118            nm, name_ind = self._mat_ni(data)
119            nm2, name_ind2 = self._mat_ni(ex)
120
121            genes = [ nm.umatch(gene) for gene in geneset ]
122            genes2 = [ nm2.umatch(gene) for gene in geneset ]
123
124            takegenes = [ i for i,a in enumerate(genes) if a != None ]
125
126            genes = [ genes[i] for i in takegenes ]
127            genes2 = [ genes2[i] for i in takegenes ]
128
129            domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
130            datao = Orange.data.Table(domain, data)
131           
132            #convert the example to the same domain
133            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
134            example = Orange.data.Instance(domain, exvalues)
135
136            return setSig_example_geneset(example, datao) #only this one is setsig specific
137     
138        at.get_value_from = t
139        return at
140
141class PLS(GeneSetTrans):
142
143    def build_feature(self, data, gs):
144
145        at = Orange.feature.Continuous(name=str(gs))
146
147        geneset = list(gs.genes)
148
149        nm, name_ind = self._mat_ni(data)
150        genes = [ nm.umatch(gene) for gene in geneset ]
151        takegenes = [ i for i,a in enumerate(genes) if a != None ]
152        genes = [ genes[i] for i in takegenes ]
153
154        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
155
156        datao = Orange.data.Table(domain, data)
157
158        xmean, W, P, T = PLSCall(datao, nc=1, y=[datao.domain.class_var])
159        constructt = xmean, W, P
160
161        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
162
163            nm2, name_ind2 = self._mat_ni(ex)
164            genes2 = [ nm2.umatch(gene) for gene in geneset ]
165            genes2 = [ genes2[i] for i in takegenes ]
166         
167            #convert the example to the same domain
168            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
169           
170            ex = numpy.array(exvalues[:-1])
171
172            xmean, W, P = constructt
173            ex = ex - xmean # same input transformation
174
175            nc = W.shape[1]
176
177            TR = numpy.empty((1, nc))
178            XR = ex
179
180            dot = numpy.dot
181
182            for i in range(nc):
183               t = dot(XR, W[:,i].T)
184               XR = XR - t*numpy.array([P[:,i]])
185               TR[:,i] = t
186
187            return TR[0][0]
188
189        at.get_value_from = t
190        return at
191
192class PCA(GeneSetTrans):
193
194    def build_feature(self, data, gs):
195
196        at = Orange.feature.Continuous(name=str(gs))
197
198        geneset = list(gs.genes)
199
200        nm, name_ind = self._mat_ni(data)
201        genes = [ nm.umatch(gene) for gene in geneset ]
202        takegenes = [ i for i,a in enumerate(genes) if a != None ]
203        genes = [ genes[i] for i in takegenes ]
204
205        domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
206
207        datao = Orange.data.Table(domain, data)
208
209        evals, evect, xmean = pca(datao)
210        constructt = evals, evect, xmean
211
212        def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
213
214            nm2, name_ind2 = self._mat_ni(ex)
215            genes2 = [ nm2.umatch(gene) for gene in geneset ]
216            genes2 = [ genes2[i] for i in takegenes ]
217         
218            #convert the example to the same domain
219            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
220           
221            arr = numpy.array(exvalues[:-1])
222           
223            evals, evect, xmean = constructt
224
225            arr = arr - xmean # same input transformation - a row in a matrix
226            ev0 = evect[0] #this is a row in a matrix - do a dot product
227            a = numpy.dot(arr, ev0)
228
229            return a
230                 
231        at.get_value_from = t
232        return at
233
234class SimpleFun(GeneSetTrans):
235
236    def build_feature(self, data, gs):
237
238        at = Orange.feature.Continuous(name=str(gs))
239
240        def t(ex, w, gs=gs):
241            geneset = list(gs.genes)
242            nm2, name_ind2 = self._mat_ni(ex)
243            genes2 = [ nm2.umatch(gene) for gene in geneset ]
244           
245            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
246            exvalues = filter(lambda x: x != "?", exvalues)
247
248            return self.fn(exvalues)
249     
250        at.get_value_from = t
251        return at
252
253class Mean(SimpleFun):
254
255    def __init__(self, **kwargs):
256       self.fn = numpy.mean
257       super(Mean, self).__init__(**kwargs)
258
259class Median(SimpleFun):
260
261    def __init__(self, **kwargs):
262       self.fn = numpy.median
263       super(Median, self).__init__(**kwargs)
264
265if __name__ == "__main__":
266
267    data = Orange.data.Table("iris")
268    gsets = obiGeneSets.collections({
269        "ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
270        "f3": ['sepal length', 'sepal width', 'petal length'],
271        "l3": ['sepal width', 'petal length', 'petal width'],
272        })
273    matcher = obiGene.matcher([])
274    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
275
276    """
277    data = Orange.data.Table("DLBCL_200a")
278    gsets = obiGeneSets.collections((("KEGG",),"9606"))
279    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
280    choosen_cv = None
281    """
282
283    def to_old_dic(d, data):
284        ar = defaultdict(list)
285        for ex1 in data:
286            ex = d(ex1)
287            for a,x in zip(d.attributes, ex):
288                ar[a.name].append(x.value)
289        return ar
290
291    def pp2(ar):
292        ol =  sorted(ar.items())
293        print '\n'.join([ a + ": " +str(b) for a,b in ol])
294
295    ass = PLS(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
296    ar = to_old_dic(ass.domain, data[:5])
297    pp2(ar)
Note: See TracBrowser for help on using the repository browser.