source: orange-bioinformatics/obiGeneSetSig.py @ 1587:c29c74349ecc

Revision 1587:c29c74349ecc, 5.5 KB checked in by markotoplak, 2 years ago (diff)

Added caching of gene matchers.

Line 
1import Orange
2import obiAssess
3import Orange.misc
4import obiGeneSets
5import obiGene
6import numpy
7from collections import defaultdict
8import stats
9import obiGsea
10import scipy.stats
11
12
13
14def setSig_example_geneset(ex, data):
15    """ Gets learning data and example with the same domain, both
16    containing only genes from the gene set. """
17
18    distances = [ [], [] ]   
19
20    def pearsonr(v1, v2):
21        return numpy.corrcoef([v1, v2])[0,1]
22
23    def pearson(ex1, ex2):
24        #leaves undefined elements out
25
26        attrs = range(len(ex1.domain.attributes))
27        vals1 = [ ex1[i].value for i in attrs ]
28        vals2 = [ ex2[i].value for i in attrs ]
29
30        common = [ True if v1 != "?" and v2 != "?" else False \
31            for v1,v2 in zip(vals1,vals2) ]
32        vals1 = [ v for v,c in zip(vals1, common) if c ]
33        vals2 = [ v for v,c in zip(vals2, common) if c ]
34
35        return numpy.corrcoef([vals1, vals2])[0,1]
36
37    def ttest(ex1, ex2):
38        try:
39            return stats.lttest_ind(ex1, ex2)[0]
40        except:
41            return 0.0
42   
43    #maps class value to its index
44    classValueMap = dict( [ (val,i) for i,val in enumerate(data.domain.class_var.values) ])
45 
46    #create distances to all learning data - save or other class
47    for c in data:
48        distances[classValueMap[c[-1].value]].append(pearson(c, ex))
49
50    return ttest(distances[0], distances[1])
51
52def mat_ni(data, matcher):
53    nm = matcher([at.name for at in data.domain.attributes])
54    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes))
55    return nm, name_ind
56
57def select_genesets(nm, gene_sets, min_size=3, max_size=1000, min_part=0.1):
58    """ Returns a list of gene sets that have sizes in limits """
59
60    def ok_sizes(gs):
61        """compares sizes of genesets to limitations"""
62        transl = filter(lambda x: x != None, [ nm.umatch(gene) for gene in gs.genes ])
63        if len(transl) >= min_size \
64            and len(transl) <= max_size \
65            and float(len(transl))/len(gs.genes) >= min_part:
66            return True
67        return False
68
69    return filter(ok_sizes, gene_sets) 
70
71class SetSig(object):
72
73    __new__ = Orange.misc._orange__new__(object)
74
75    def __init__(self, matcher, gene_sets, min_size=3, max_size=1000, min_part=0.1, class_values=None):
76        self.matcher = matcher
77        self.gene_sets = gene_sets
78        self.min_size = min_size
79        self.max_size = max_size
80        self.min_part = min_part
81        self.class_values = class_values
82        self._cache = {}
83
84    def _mat_ni(self, data):
85        """ With cached gene matchers. """
86        if data.domain not in self._cache:
87            self._cache[data.domain] = mat_ni(data, self.matcher)
88        return self._cache[data.domain]
89
90    def __call__(self, data, weight_id=None):
91
92        data = obiGsea.takeClasses(data, classValues=self.class_values)
93        nm,_ =  self._mat_ni(data)
94        gene_sets = select_genesets(nm, self.gene_sets, self.min_size, self.max_size, self.min_part)
95
96        attributes = []
97
98        for gs in gene_sets:
99            at = Orange.feature.Continuous(name=str(gs))
100
101            def t(ex, w, gs=gs, data=data): #copy od the data
102                geneset = list(gs.genes)
103
104                nm, name_ind = self._mat_ni(data)
105                nm2, name_ind2 = self._mat_ni(ex)
106
107                genes = [ nm.umatch(gene) for gene in geneset ]
108                genes2 = [ nm2.umatch(gene) for gene in geneset ]
109
110                genes, genes2 = zip(*[ (g,g2) for g,g2 in zip(genes, genes2) if g != None])
111
112                domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
113                datao = Orange.data.Table(domain, data)
114
115                def vou(ex, gn, indices):
116                    """ returns the value or unknown for the given gene name"""
117                    if gn not in indices:
118                        return "?"
119                    else:
120                        return ex[indices[gn]].value
121               
122                #convert the example to the same domain
123                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
124                example = Orange.data.Instance(domain, exvalues)
125
126                return setSig_example_geneset(example, datao) #only this one is setsig specific
127         
128            at.get_value_from = t
129            attributes.append(at)
130       
131        newdomain = Orange.data.Domain(attributes, data.domain.class_var)
132        return Orange.data.Table(newdomain, data)
133
134if __name__ == "__main__":
135
136
137    data = Orange.data.Table("iris")
138    gsets = obiGeneSets.collections({
139        "ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
140        "f3": ['sepal length', 'sepal width', 'petal length'],
141        "l3": ['sepal width', 'petal length', 'petal width'],
142        })
143    matcher = obiGene.matcher([])
144    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
145
146    """
147    data = Orange.data.Table("DLBCL_200a")
148    gsets = obiGeneSets.collections((("KEGG",),"9606"))
149    matcher = obiGene.matcher([obiGene.GMKEGG("hsa")])
150    choosen_cv = None
151    """
152
153
154    def to_old_dic(d, data):
155        ar = defaultdict(list)
156        for ex1 in data:
157            ex = d(ex1)
158            for a,x in zip(d.attributes, ex):
159                ar[a.name].append(x.value)
160        return ar
161
162    def pp2(ar):
163        ol =  sorted(ar.items())
164        print '\n'.join([ a + ": " +str(b) for a,b in ol])
165
166    ass = SetSig(matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
167    ass = ass(data)
168    ar = to_old_dic(ass.domain, data[:5])
169    pp2(ar)
Note: See TracBrowser for help on using the repository browser.