source: orange-bioinformatics/obiGeneSetSig.py @ 1594:98ff35c02b40

Revision 1594:98ff35c02b40, 7.8 KB checked in by markotoplak, 2 years ago (diff)

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
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 PCA(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        evals, evect, xmean = pca(datao)
159        constructt = evals, evect, xmean
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            arr = numpy.array(exvalues[:-1])
171           
172            evals, evect, xmean = constructt
173
174            arr = arr - xmean # same input transformation - a row in a matrix
175            ev0 = evect[0] #this is a row in a matrix - do a dot product
176            a = numpy.dot(arr, ev0)
177
178            return a
179                 
180        at.get_value_from = t
181        return at
182
183class SimpleFun(GeneSetTrans):
184
185    def build_feature(self, data, gs):
186
187        at = Orange.feature.Continuous(name=str(gs))
188
189        def t(ex, w, gs=gs):
190            geneset = list(gs.genes)
191            nm2, name_ind2 = self._mat_ni(ex)
192            genes2 = [ nm2.umatch(gene) for gene in geneset ]
193           
194            exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
195            exvalues = filter(lambda x: x != "?", exvalues)
196
197            return self.fn(exvalues)
198     
199        at.get_value_from = t
200        return at
201
202class Mean(SimpleFun):
203
204    def __init__(self, **kwargs):
205       self.fn = numpy.mean
206       super(Mean, self).__init__(**kwargs)
207
208class Median(SimpleFun):
209
210    def __init__(self, **kwargs):
211       self.fn = numpy.median
212       super(Median, self).__init__(**kwargs)
213
214if __name__ == "__main__":
215
216    data = Orange.data.Table("iris")
217    gsets = obiGeneSets.collections({
218        "ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
219        "f3": ['sepal length', 'sepal width', 'petal length'],
220        "l3": ['sepal width', 'petal length', 'petal width'],
221        })
222    matcher = obiGene.matcher([])
223    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
224
225    """
226    data = Orange.data.Table("DLBCL_200a")
227    gsets = obiGeneSets.collections((("KEGG",),"9606"))
228    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
229    choosen_cv = None
230    """
231
232    def to_old_dic(d, data):
233        ar = defaultdict(list)
234        for ex1 in data:
235            ex = d(ex1)
236            for a,x in zip(d.attributes, ex):
237                ar[a.name].append(x.value)
238        return ar
239
240    def pp2(ar):
241        ol =  sorted(ar.items())
242        print '\n'.join([ a + ": " +str(b) for a,b in ol])
243
244    ass = Mean(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
245    ar = to_old_dic(ass.domain, data[:5])
246    pp2(ar)
Note: See TracBrowser for help on using the repository browser.