source: orange-bioinformatics/Orange/bioinformatics/obiDifscale.py @ 1625:cefeb35cbfc9

Revision 1625:cefeb35cbfc9, 12.6 KB checked in by mitar, 2 years ago (diff)

Moving files around.

Line 
1import random
2from math import log
3from operator import itemgetter
4
5import numpy
6
7import Orange
8import obiGEO
9from obiExpression import ExpressionSignificance_Test
10
11
12# Utility functiions
13
14log2 = lambda x: log(x, 2.)
15
16def my_ratio(x, y):
17    """ compute the log-ratio """
18    return log2(x/y)
19
20def sign(x):
21    return cmp(x, 0)
22
23def common_domain(data1, data2):
24    """Use only attributes that are in both data sets"""
25    atts = sorted(set([a.name for a in data1.domain.attributes]).intersection(
26               [a.name for a in data2.domain.attributes]))
27    new1 = Orange.data.Table(Orange.data.Domain(atts + [data1.domain.classVar], data1.domain), data1)
28    new2 = Orange.data.Table(Orange.data.Domain(atts + [data2.domain.classVar], data2.domain), data2)
29    return new1, new2
30
31def common_genes(data1, data2, att='gene'):
32    common_set = list(set(ex[att].value for ex in data1).intersection(ex[att].value for ex in data2))
33    print len(set(ex[att].value for ex in data1))
34    print len(common_set)
35    return data1.filter(**{att: common_set}), data2.filter(**{att: common_set})
36
37# Normalization
38
39def quantilenorm(data):
40    """normalization of microarray data to obtain the same distribution in all chips"""
41    chips = data.domain.attributes
42    genes = [str(d["gene"]) for d in data]
43    d_sort = {}
44    for c in chips:
45        dc = [(float(d[c]), str(d["gene"])) for d in data]
46        dc.sort()
47        d_sort[c.name] = dc
48    d_sort2 = dict([(str(d["gene"]),{}) for d in data])
49    for i in range(len(data)):
50        genes_c = [(d_sort[c.name][i][1], c.name) for c in chips]
51        mean_row = numpy.mean([d_sort[c.name][i][0] for c in chips])
52        for g, d in genes_c:
53            d_sort2.get(g, {}).update(dict([(d, mean_row)]))
54    data_norm = Orange.data.Table(data.domain)
55    for i, d in enumerate(data):
56        g = str(d["gene"])
57        ex = [d_sort2[g][c.name] for c in chips]
58        data_norm.append(ex)
59        data_norm[i]["gene"] = g
60    return data_norm
61
62def get_mediannorm(modeldata):
63    """returns the function for median scale normalization"""
64    globalmedian = numpy.median([float(d[c]) for d in modeldata for c in modeldata.domain.attributes])
65    def normalize(data):
66        chips = data.domain.attributes
67        medians = [numpy.median([float(d[c]) for d in data]) for c in chips]
68        data_norm = Orange.data.Table(data.domain)
69        data_norm.domain.add_metas(data.domain.get_metas())
70        for i, d in enumerate(data):
71            ex = [(d[c.name].value*globalmedian)/medians[k] for k,c in enumerate(chips)]
72            data_norm.append(ex)
73            for m in data.domain.get_metas():
74                data_norm[i][m] = data[i][m]
75        return data_norm
76    return normalize
77
78def medianscalenorm(data, newsamples=None, intersection=True):
79    """normalization of microarray data to center the distributions in all chips on the same global median"""
80    if not newsamples:
81        n = get_mediannorm(data)
82        return n(data), None
83    else:
84        if intersection:
85            idata, inewsamples = common_genes(data, newsamples)
86        n = get_mediannorm(idata)
87        return n(idata), n(inewsamples)
88
89def normalize(data1, data2=None, type='median'):
90    if type == 'quantile':
91        ndata1 = quantilenorm(data1)
92        return ndata1, medianscalenorm(ndata1, data2)[1] if data2 else None
93    elif type == 'median':
94        return medianscalenorm(data1, data2)
95    else:
96        return Error
97
98
99# Gene filtering
100
101def compute_diff_pairs(at_a, d):
102    """ computes the pairwise differences between the replicates """
103    differ = []
104    for i in range(len (at_a)-1):
105        for j in range(i+1, len(at_a)):
106            differ.append(log2(d[at_a[i]])-log2(d[at_a[j]]))
107    return differ
108
109def costruct_series(data, attr_set, differences=False):
110    """ Constructs the time series by averaging the replicates of the same time point """
111    serie = dict([(str(d["gene"]), [numpy.mean([d[at] for at in at_samples])
112                    for at_name, at_samples in attr_set]) for d in data])
113    if differences:
114        """store the differences between replicates while creating the time series"""
115        differences = []
116        r = [[differences.extend(compute_diff_pairs(at_samples, d))
117              for at_name, at_samples in attr_set] for d in data]
118        return serie, differences
119    else:
120        return serie
121
122def costruct_control_series(serie, num_points, control):
123    """ Creation of the control series (for each gene) as a constant profile equal to expression at time 0 """
124    ##print "control series..."
125    serie_0 = {}
126    if control == "avg":
127        serie_0.update(dict([(g,[numpy.mean(serie[g]) for i in range(num_points)]) for g in serie]))
128    elif control == "t0":
129        serie_0.update(dict([(g,[serie[g][0] for i in range(num_points)]) for g in serie]))
130    return serie_0
131
132def compute_area(vals, t, baseline=0.):
133    return sum(a_pair((t[i], vals[i]), (t[i+1], vals[i+1]), baseline) \
134        for i in xrange(len(vals)-1))
135
136def a_pair(p1, p2, baseline):
137    """Area under the line bounded by a pair of two points (x,y) with
138    respect to baseline. Both parts around the diagonal are
139    positive."""
140    x1,y1 = p1
141    x2,y2 = p2
142    x2 = x2-x1 #same start
143    x1 = 0.0
144    a = y1-baseline
145    b = y2-baseline
146    if a*b >= 0: #both on one side
147        return abs(x2*(b+a)/2.0)
148    else:
149        xp = -a * x2 / float(y2-y1)
150        return (abs(xp * a) + abs((x2-xp)*b)) / 2.0
151
152def uniform_time_scale(attr_set):
153    """ Obtains time points with a unique measure (the lowest among [min,h,d]) present in data"""
154    ord_measures = ["min","h","d"]
155    converter = {"min":{"h":60, "d":24*60}, "h": {"d":24}}
156    measures = list(set([t.split(" ")[1] for t,s in attr_set]))
157    if len(measures) == 1:
158        time_points = [float(t.split(" ")[0]) for t,s in attr_set]
159    else:
160        first_measure = min([(ord_measures.index(m),m) for m in measures])[1]
161        time_points = [float(t.split(" ")[0]) for t, s in attr_set
162                       if t.split(" ")[1] == first_measure]
163        time_points.extend([float(t.split(" ")[0]) * converter[first_measure][t.split(" ")[1]]
164                            for t, s in attr_set if t.split(" ")[1] != first_measure])
165    return time_points
166
167def AREA(data, attr_set, control='t0', weighted=False, auto=False, perc=99):
168    """ AREA (Area Under the Curve) filtering method """
169    from matplotlib.mlab import prctile
170    if weighted:
171        time_points = uniform_time_scale(attr_set)
172    else:
173        time_points = range(len(attr_set))
174
175    # Monte Carlo approach to create the null distribution of the areas
176    if auto: # null distribution
177        serie, differences = costruct_series(data, attr_set, auto)
178        serie_campionata = {}
179        area = []
180        for i in range(20000):
181            serie_campionata = random.sample(differences, len(time_points))
182            area.append(compute_area(serie_campionata, time_points))
183        area_threshold = prctile(area, perc)
184    else:
185        serie = costruct_series(data, attr_set, auto)
186
187    serie_0 = costruct_control_series(serie, len(time_points), control)
188
189    # Gene filtering
190    areas = []
191    for g in serie:
192        diff_s = [log2(serie[g][i]) - log2(serie_0[g][i]) for i in range(len(time_points))]
193        area_diff = compute_area(diff_s, time_points);
194        if not auto or area_diff > area_threshold:
195            areas.append((g, area_diff))
196    return areas
197
198class ExpressionSignificance_AREA(ExpressionSignificance_Test):
199    def __call__(self, target=None):
200        attr_set = {}
201        for a in self.data.domain.attributes:
202            attr_set[a.attributes['time']] = attr_set.get(a.attributes['time'], []) + [a.name]
203        scores = AREA(self.data, sorted(attr_set.items()))
204        gene2ind = dict((g, i) for i,g in enumerate(ex['gene'].value for ex in self.data))
205        return [(gene2ind[g], s) for g, s in scores]
206
207def FC(data, attr_set, control='t0', thr=2, auto=False, p_thr=0.2):
208    """ Gene filtering based on the number of FC of all time points with the control series > thr """
209    serie = costruct_series(data, attr_set, False)
210    num_points = len(serie.values()[0])
211    serie_0 = costruct_control_series(serie, num_points, control)
212    fc = [(g, len([0 for i in range(num_points) if abs(my_ratio(s[i], serie_0[g][i])) >= thr]))
213          for g, s in serie.items()]
214    if auto:
215        thr_points = round(p_thr*num_points)
216        fc = [(g, v) for g, v in fc if v >= thr_points]
217    return fc
218
219class ExpressionSignificance_FCts(ExpressionSignificance_Test):
220    def __call__(self, target=None):
221        attr_set = {}
222        for a in self.data.domain.attributes:
223            attr_set[a.attributes['time']] = attr_set.get(a.attributes['time'], []) + [a.name]
224        scores = FC(self.data, sorted(attr_set.items()))
225        return zip(self.keys, map(itemgetter(1), scores))
226
227def spearmanr_filter(data, limit=1000):
228    """ Spearman ranks gene filtering """
229    from scipy.stats import spearmanr
230    time = [a.attributes['time'] for a in data.domain.attributes]
231    exs = sorted(data, reverse=True, key=lambda ex: abs(spearmanr([a.value for a in ex], time)[0]))
232    return [str(ex['gene']) for ex in exs[:limit]]
233
234
235def signed_PCA(data):
236    pca = Orange.projection.pca.Pca(data, standardize=False, max_components=1)
237    classifier = lambda X: [x[0].value for x in pca(X)]
238    predictions = classifier(data)
239    classes = [ex.getclass().value for ex in data]
240    n = 0
241    for i1,c1 in enumerate(classes):
242        for i2,c2 in enumerate(classes[:i1]):
243            n += cmp(c1,c2) * cmp(predictions[i1], predictions[i2])
244    if n < 0:
245        def invert(X):
246            y = classifier(X)
247            return -y if type(y) == float else [-x for x in y]
248        return invert
249    else:
250        return classifier
251
252signed_PCA.name = 'PCA'
253
254
255def conttime(data, d):
256    for a in data.domain.attributes:
257        a.attributes['time'] = d[a.attributes['time']]
258
259def conv(attr_set, ticks=True):
260    """Obtain time points with a unique measure (the lowest among [min,h,d]) present in data"""
261    ord_measures = ["min","h","d"]
262    converter = {"min":{"h":60, "d":24*60}, "h": {"d":24}}
263    measures = list(set([t.split(" ")[1] for t in attr_set]))
264    if len(measures) == 1:
265        time_points = [(t, float(t.split(" ")[0])) for t in attr_set]
266    else:
267        first_measure = min([(ord_measures.index(m),m) for m in measures])[1]
268        time_points = [(t, float(t.split(" ")[0])) for t in attr_set if t.split(" ")[1] == first_measure]
269        time_points.extend([(t, float(t.split(" ")[0]) * converter[first_measure][t.split(" ")[1]])
270                            for t in attr_set if t.split(" ")[1] != first_measure])
271    time_points.sort(key=itemgetter(1))
272    if ticks:
273        time_points = [(t[0],float(i)) for i,t in enumerate(time_points)]
274    return dict(time_points)
275
276
277def get_projections(data1, data2=None):
278    labels1 = list(a.attributes['time'] for a in data1.domain.attributes)
279    tdata1 = obiGEO.transpose(data1)
280    if data2:
281        labels2 = list('[%s]' % a.attributes['time'] for a in data2.domain.attributes)
282        tdata2 = obiGEO.transpose(data2)
283        tdata1, tdata2 = common_domain(tdata1, tdata2)
284        classifier = signed_PCA(tdata1)
285        proj1 = classifier(tdata1)
286        proj2 = classifier(tdata2)
287    else:
288        classifier = signed_PCA(tdata1)
289        proj1 = classifier(tdata1)
290        proj2, labels2 = [], []
291    return proj1, labels1, proj2, labels2
292
293
294############
295
296if False and __name__ == '__main__':
297    import obiGEO
298    # Data set 1
299    data1 = obiGEO.GDS('GDS2666').getdata(report_genes=True, transpose=False)
300    labels1 = list(a.attributes['time'] for a in data1.domain.attributes)
301    attr_set = list(set(a.attributes['time'] for a in data1.domain.attributes))
302    convd = conv(attr_set)
303    conttime(data1, convd)
304    # Data set 2
305    data2 = obiGEO.GDS('GDS2667').getdata(report_genes=True, transpose=False)
306    labels2 = list(a.attributes['time'] for a in data2.domain.attributes)
307    attr_set = list(set(a.attributes['time'] for a in data2.domain.attributes))
308    convd = conv(attr_set)
309    conttime(data2, convd)
310    # Normalize data set 1
311    ndata1, _ = normalize(data1, type='quantile')
312    # Filtering
313    attr_set = {}
314    for a in ndata1.domain.attributes:
315        attr_set[a.attributes['time']] = attr_set.get(a.attributes['time'], []) + [a.name]
316    scores = AREA(ndata1, sorted(attr_set.items()))
317    genes = map(itemgetter(0), sorted(scores, key=itemgetter(1))[:1000])
318    fdata1 = ndata1.filter(gene=genes)
319    # Rescale data set 2
320    ttrain, ttest = normalize(fdata1, data2)
321    # Model construction and prediction
322    # train = obiGEO.transpose(ttrain)
323    # test = obiGEO.transpose(ttest)
324    # cdtrain, cdtest = common_domain(train, test)
325    # classifier = signed_PCA(cdtrain)
326    # proj1 = classifier(cdtrain)
327    # proj2 = classifier(cdtest)
328
329
330
Note: See TracBrowser for help on using the repository browser.