source: orange-bioinformatics/obiGeneSetSig.py @ 1593:f888bd8a5ae2

Revision 1593:f888bd8a5ae2, 8.3 KB checked in by markotoplak, 2 years ago (diff)

Added PCA transformation to obiGeneSetSig.

Line 
1import Orange
2import obiAssess
3import Orange.misc
4import obiGeneSets
5import obiGene
6import numpy
7from collections import defaultdict
8import stats
9import obiGsea
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 = obiGsea.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
99def vou(ex, gn, indices):
100    """ returns the value or "?" for the given gene name gn"""
101    if gn not in indices:
102        return "?"
103    else:
104        return ex[indices[gn]].value
105
106class SetSig(GeneSetTrans):
107
108    def build_features(self, data, gene_sets):
109
110        attributes = []
111
112        for gs in gene_sets:
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            attributes.append(at)
140
141        return attributes
142
143from obiAssess import pca
144
145class PCA(GeneSetTrans):
146
147    def build_features(self, data, gene_sets):
148
149        attributes = []
150
151        for gs in gene_sets:
152            at = Orange.feature.Continuous(name=str(gs))
153
154            geneset = list(gs.genes)
155
156            nm, name_ind = self._mat_ni(data)
157            genes = [ nm.umatch(gene) for gene in geneset ]
158            takegenes = [ i for i,a in enumerate(genes) if a != None ]
159            genes = [ genes[i] for i in takegenes ]
160   
161            domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
162
163            datao = Orange.data.Table(domain, data)
164
165            evals, evect, xmean = pca(datao)
166            constructt = evals, evect, xmean
167
168            def t(ex, w, geneset=geneset, constructt=constructt, takegenes=takegenes, domain=domain):
169
170                nm2, name_ind2 = self._mat_ni(ex)
171                genes2 = [ nm2.umatch(gene) for gene in geneset ]
172                genes2 = [ genes2[i] for i in takegenes ]
173             
174                #convert the example to the same domain
175                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
176               
177                arr = numpy.array(exvalues[:-1])
178               
179                evals, evect, xmean = constructt
180
181                arr = arr - xmean # same input transformation - a row in a matrix
182                ev0 = evect[0] #this is a row in a matrix - do a dot product
183                a = numpy.dot(arr, ev0)
184
185                return a
186                     
187            at.get_value_from = t
188            attributes.append(at)
189
190        return attributes
191
192class SimpleFun(GeneSetTrans):
193
194    def build_features(self, data, gene_sets):
195
196        attributes = []
197
198        for gs in gene_sets:
199            at = Orange.feature.Continuous(name=str(gs))
200
201            def t(ex, w, gs=gs):
202                geneset = list(gs.genes)
203                nm2, name_ind2 = self._mat_ni(ex)
204                genes2 = [ nm2.umatch(gene) for gene in geneset ]
205               
206                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
207                exvalues = filter(lambda x: x != "?", exvalues)
208
209                return self.fn(exvalues)
210         
211            at.get_value_from = t
212            attributes.append(at)
213
214        return attributes
215
216class Mean(SimpleFun):
217
218    def __init__(self, **kwargs):
219       self.fn = numpy.mean
220       super(Mean, self).__init__(**kwargs)
221
222class Median(SimpleFun):
223
224    def __init__(self, **kwargs):
225       self.fn = numpy.median
226       super(Median, self).__init__(**kwargs)
227
228
229
230if __name__ == "__main__":
231
232    data = Orange.data.Table("iris")
233    gsets = obiGeneSets.collections({
234        "ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
235        "f3": ['sepal length', 'sepal width', 'petal length'],
236        "l3": ['sepal width', 'petal length', 'petal width'],
237        })
238    matcher = obiGene.matcher([])
239    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
240
241    """
242    data = Orange.data.Table("DLBCL_200a")
243    gsets = obiGeneSets.collections((("KEGG",),"9606"))
244    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
245    choosen_cv = None
246    """
247
248    def to_old_dic(d, data):
249        ar = defaultdict(list)
250        for ex1 in data:
251            ex = d(ex1)
252            for a,x in zip(d.attributes, ex):
253                ar[a.name].append(x.value)
254        return ar
255
256    def pp2(ar):
257        ol =  sorted(ar.items())
258        print '\n'.join([ a + ": " +str(b) for a,b in ol])
259
260    ass = PCA(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
261    ar = to_old_dic(ass.domain, data[:5])
262    pp2(ar)
Note: See TracBrowser for help on using the repository browser.