source: orange-bioinformatics/obiGeneSetSig.py @ 1583:2b90e7c95085

Revision 1583:2b90e7c95085, 5.2 KB checked in by markotoplak, 2 years ago (diff)

Fixed a bug in obiGene (umatch did not work in Match objects). Gene set signatures method now accept examples from different domains on the input.

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 pearsonr(vals1, vals2)
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 SetSig(object):
69
70    __new__ = Orange.misc._orange__new__(object)
71
72    def __init__(self, matcher, gene_sets, min_size=3, max_size=1000, min_part=0.1, class_values=None):
73        self.matcher = matcher
74        self.gene_sets = gene_sets
75        self.min_size = min_size
76        self.max_size = max_size
77        self.min_part = min_part
78        self.class_values = class_values
79
80    def __call__(self, data, weight_id=None):
81
82        data = obiGsea.takeClasses(data, classValues=self.class_values)
83        nm,_ =  mat_ni(data, matcher)
84        gene_sets = select_genesets(nm, self.gene_sets, self.min_size, self.max_size, self.min_part)
85
86        attributes = []
87
88        for gs in gene_sets:
89            at = Orange.feature.Continuous(name=gs.id)
90
91            def t(ex, w, gs=gs, data=data): #copy od the data
92                geneset = list(gs.genes)
93
94                nm, name_ind = mat_ni(data, matcher)
95                nm2, name_ind2 = mat_ni(ex, matcher)
96
97                genes = [ nm.umatch(gene) for gene in geneset ]
98                genes2 = [ nm2.umatch(gene) for gene in geneset ]
99
100                genes, genes2 = zip(*[ (g,g2) for g,g2 in zip(genes, genes2) if g != None])
101
102                domain = Orange.data.Domain([data.domain.attributes[name_ind[gene]] for gene in genes], data.domain.class_var)
103                datao = Orange.data.Table(domain, data)
104
105                def vou(ex, gn, indices):
106                    """ returns the value or unknown for the given gene name"""
107                    if gn not in indices:
108                        return "?"
109                    else:
110                        return ex[indices[gn]].value
111               
112                #convert the example to the same domain
113                exvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] + [ "?" ]
114                example = Orange.data.Instance(domain, exvalues)
115
116                return setSig_example_geneset(example, datao) #only this one is setsig specific
117         
118            at.get_value_from = t
119            attributes.append(at)
120       
121        newdomain = Orange.data.Domain(attributes, data.domain.class_var)
122        return Orange.data.Table(newdomain, data)
123
124if __name__ == "__main__":
125
126    data = Orange.data.Table("iris")
127    gsets = obiGeneSets.collections({
128        "ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'],
129        "f3": ['sepal length', 'sepal width', 'petal length'],
130        "l3": ['sepal width', 'petal length', 'petal width'],
131        })
132
133    fp = 120
134    ldata = Orange.data.Table(data.domain, data[:fp])
135    tdata = Orange.data.Table(data.domain, data[fp:])
136
137    matcher = obiGene.matcher([])
138
139    choosen_cv = ["Iris-setosa", "Iris-versicolor"]
140
141    def to_old_dic(d, data):
142        ar = defaultdict(list)
143        for ex1 in data:
144            ex = d(ex1)
145            for a,x in zip(d.attributes, ex):
146                ar[a.name].append(x.value)
147        return ar
148
149    def pp2(ar):
150        ol =  sorted(ar.items())
151        print '\n'.join([ a + ": " +str(b) for a,b in ol])
152
153    ass = SetSig(matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0)
154    ass = ass(data)
155    ar = to_old_dic(ass.domain, data[:5])
156    pp2(ar)
Note: See TracBrowser for help on using the repository browser.