source: cobios/datasets.py @ 321

Last change on this file since 321 was 4, checked in by jjbot, 11 years ago

Initial import of new code base.

File size: 13.6 KB
Line 
1large_scale = ['tandem affinity purification','coimmunoprecipitation',\
2   'affinity chromatography technology',\
3   'ubiquitin reconstruction','two hybrid','mass spectrometry studies of complexes']
4
5small_scale_label = ['protease accessibility laddering',
6'experimental interaction detection','electron tomography','copurification',\
7'anti tag coimmunoprecipitation','cosedimentation in solution','cross-linking study',\
8'electron microscopy','colocalization by immunostaining','far western blotting',\
9'light scattering','molecular sieving','other biochemical technologies',\
10'pull down','protein kinase assay','comigration in sds page','cosedimentation',\
11'confirmational text mining','filter binding','biophysical','cosedimentation through density gradient',\
12'protein cross-linking with a bifunctional reagent','fluorescence technology','nuclear magnetic resonance',\
13'phage display','surface plasmon resonance','x-ray crystallography','biochemical',\
14'electrophoretic mobility shift assay']
15
16
17
18large_scale_complex = ['affinity chromatography technology','tandem affinity purification','coimmunoprecipitation',\
19'mass spectrometry studies of complexes'] 
20
21small_scale_complex_label = ['confirmational text mining','pull down','anti tag coimmunoprecipitation',\
22'molecular sieving','copurification','x-ray crystallography','cosedimentation through density gradient']
23
24import numpy
25import random
26from scipy.stats import pearsonr
27
28def load_pp(d,neg):
29    ssid = d.term[_.name.within(*small_scale_label)].term_id()
30    lsid = d.term[_.name.within(*large_scale)].term_id()
31    cssid = d.term[_.name.within(*small_scale_complex_label)].term_id()
32    clsid = d.term[_.name.within(*large_scale_complex)].term_id()
33    proteins = d.taxon_name[_.name == "yeast"].bioentrys.require_type("gene").children.require_type("CDS")
34    complex = d.taxon_name[_.name == "yeast"].bioentrys.require_type("complex")
35    protdata = proteins.bioentry_id()
36   
37    protrel = d.biorelationship[_.source_bioentry_id.within(*protdata) & _.target_bioentry_id.within(*protdata) & _.type_id.within(*(numpy.hstack((ssid,lsid))))].realize()
38    complexrel = d.biorelationship[_.source_bioentry_id.within(*complex.bioentry_id()) & _.target_bioentry_id.within(*protdata) & _.type_id.within(*(numpy.hstack((cssid,clsid))))].realize()
39
40    sp = protrel[_.type_id.within(*ssid)]()
41    ls = protrel[_.type_id.within(*lsid)]()
42    csp = protrel[_.type_id.within(*cssid)]()
43    cls = protrel[_.type_id.within(*clsid)]()
44
45    print "LS",len(ls)
46    print "SP",len(sp)
47    print "CSP",len(csp)
48    print "CLS",len(cls)
49
50
51    sps = []
52    for s in sp:
53        key = (s.source_bioentry_id,s.target_bioentry_id)
54        if(key[0] > key[1]):
55            key = (key[1],key[0])
56        sps.append(key)
57    sps = set(sps)
58
59
60    cspmap = dict()
61    for x in csp:
62        if(x.source_bioentry_id in cspmap):
63            cspmap[x.source_bioentry_id].append(x.target_bioentry_id)
64        else:
65            cspmap[x.source_bioentry_id] = [x.target_bioentry_id]
66       
67
68    for k,v in cspmap.iteritems():
69        for i in v:
70            for j in v:
71                if(i < j):
72                    sps.add((i,j))
73
74    r,nu = sample_edges_from_full_graph(protdata,neg,sps)
75
76    desc = list(sps) + list(r)
77    label = numpy.array([1] * len(sps) + [0] * len(r),dtype=bool)
78    feat = numpy.zeros((len(desc),len(lsid)),dtype="uint32")
79   
80    for i,d in enumerate(desc):
81       if(d[0] > d[1]):
82          d = (d[1],d[0])
83          desc[i] = d
84   
85    cslmap = dict()
86    csltype = dict()
87    for x in cls:
88        if(x.source_bioentry_id in cslmap):
89            cslmap[x.source_bioentry_id].append(x.target_bioentry_id)
90        else:
91            cslmap[x.source_bioentry_id] = [x.target_bioentry_id]
92       
93        if(x.type_id in csltype):
94            csltype[x.type_id].add(x.source_bioentry_id)
95        else:
96            csltype[x.type_id] = set([x.source_bioentry_id])
97
98
99    print "C1",len(csltype)
100    print "C2",len(cslmap)
101    desc_index = dict([(x,i) for i,x in enumerate(desc)])
102   
103    for featid,l in enumerate(lsid):
104        print l
105        tp = ls[ls.type_id == l]
106        for x in tp:
107            key = (x.source_bioentry_id,x.target_bioentry_id)
108            if(key[0] > key[1]):
109                key = (key[1],key[0])
110
111            if(key in desc_index):
112                idx = desc_index[key]
113                feat[idx,featid] += 1
114
115        if(l in csltype):
116            for c in csltype[l]:
117                v = cslmap[c]
118                for i in v:
119                    for j in v:
120                        if(i < j):
121                            key = (i,j)
122                            if(key in desc_index):
123                                idx =  desc_index[key]
124                                feat[idx,featid] += 1
125   
126    desc = numpy.array(desc)
127   
128   
129    desc_index = dict()
130    for t in ls:
131        key = (t.source_bioentry_id, t.target_bioentry_id)
132        if(key[0] > key[1]):
133            key = (key[1],key[0])
134        if(key in desc_index):
135            desc_index[key].append(t.type_id)
136        else:
137            desc_index[key] = [t.type_id]
138    print "F",len(desc_index)
139   
140    for k,v in csltype.iteritems():
141        for c in v:
142            bi = cslmap[c]
143            for i in bi:
144                for j in bi:
145                    if(i < j):
146                        key = (i,j)
147                        if(key in desc_index):
148                            desc_index[key].append(k)
149                        else:
150                            desc_index[key] = [k]
151    print "L",len(desc_index)
152
153    feattest_index = dict([(l,i) for i,l in enumerate(lsid)])
154    feattest = numpy.zeros((len(desc_index),len(lsid)))
155    testdesc = numpy.zeros((len(desc_index),2))
156
157    for i,(k,v) in enumerate(desc_index.iteritems()):
158        testdesc[i,:] = k
159        for t in v:
160            feattest[i,feattest_index[t]] += 1
161   
162    #to genes
163    g_to_p = proteins.parents_rel.require_type("gene_part")()
164
165    g_to_p_map = dict([(x.target_bioentry_id,x.source_bioentry_id) for x in g_to_p])
166
167    tmp = desc.ravel()
168    for i,j in enumerate(tmp):
169        tmp[i] = g_to_p_map[j]
170    tmp = testdesc.ravel()
171    for i,j in enumerate(tmp):
172        tmp[i] = g_to_p_map[j]
173
174
175    return ((label,desc,feat),(testdesc,feattest))
176
177
178def load_da(d):
179   genes = d.taxon_name[_.name == "yeast"].bioentrys.require_type("gene",'unknown')
180   
181   g_to_o = genes.children_rel.require_type("gene_part")()
182   g_to_o_map = dict([(x.target_bioentry_id,x.source_bioentry_id) for x in g_to_o])
183
184   tfrel = d.bioentry.require_type("CDS").children_rel.require_type("transfac regulation")()
185
186   
187   desc = numpy.zeros((len(tfrel),2))
188   for i,tf in enumerate(tfrel):
189       desc[i,0] = g_to_o_map[tf.source_bioentry_id]
190       desc[i,1] = g_to_o_map[tf.target_bioentry_id]
191
192
193   descset = set([(j[0],j[1]) for j in desc])
194
195   
196
197   tfrel2 = d.bioentry.require_type("CDS").children_rel.require_type("transcription factor regulation")
198   val1 = tfrel2.value.require_type("binding-site p-value")()
199   val2 = tfrel2.value.require_type("conservation cutoff")()
200   val1m = dict([(x.biorelationship_id,x.value) for x in val1])
201   val2m = dict([(x.biorelationship_id,x.value) for x in val2])
202
203   
204   desctest = numpy.zeros((len(tfrel2),2))
205   feattest = numpy.zeros((len(tfrel2),2))
206   tfrel2 = tfrel2()
207
208   for i,tf2 in enumerate(tfrel2):
209       desctest[i,:] = (g_to_o_map[tf2.source_bioentry_id],g_to_o_map[tf2.target_bioentry_id])
210       feattest[i,:] = (val1m[tf2.biorelationship_id],val2m[tf2.biorelationship_id])
211
212   dm = dict()
213   for i,j in enumerate(desctest):
214       key = (j[0],j[1])
215       if(key in descset):
216          dm[key] = feattest[i,:]
217
218   feattrain = numpy.zeros((len(descset),2))
219   desctrain = numpy.zeros((len(descset),2))
220   for i,k in enumerate(descset):
221       desctrain[i,:] = k
222       if(k in dm):
223           feattrain[i,:] = dm[k]
224
225   
226   r,nu = sample_edges_from_full_graph(genes.bioentry_id(),100000,set(map(tuple,desctrain.tolist())))
227   desctestmap = dict([((x[0],x[1]),i) for i,x in enumerate(desctest)])
228   feattrainneg = numpy.zeros((len(r),2))
229   for i,x in enumerate(r):
230       key = (x[0],x[1])
231       if(key in desctestmap):
232            feattrainneg[i,:] = feattest[desctestmap[key],:]
233 
234
235   r = numpy.array(list(r))
236   label = numpy.array([1] * len(desctrain) + [0] * len(r))
237   desctrain = numpy.vstack((desctrain,numpy.array(r)))
238   feattrain = numpy.vstack((feattrain,feattrainneg))
239
240   return ((label,desctrain,feattrain),(desctest,feattest))
241
242
243def load_cor(d,dx):
244   genes = d.taxon_name[_.name == "yeast"].bioentrys.require_type("gene",'unknown')
245   datasets = (+genes.data.datasets[_.description == 'Microarray dataset'].dataset_id)()
246
247   gval = genes()
248
249   gmap = dict()
250   for i,g in enumerate(gval):
251      gmap[g.bioentry_id] = i
252
253   val = numpy.zeros((len(gval),len(gval)),dtype="float32")
254   for ds in datasets[dx:(dx+1)]:
255        data = genes.data[_.dataset_id == ds]()
256
257        maxlen = max(map(len,data[:,-1]))
258        print "ML",maxlen
259        for d in data:
260            if(len(d.data) < maxlen):
261                t = numpy.zeros((maxlen,),dtype="float32")
262                t[:len(d.data)] = d.data
263                t[len(d.data):] = numpy.NaN
264            else:
265                t = numpy.array(d.data,dtype="float32")
266            t[t<-1000] = numpy.NaN
267            d[0] = gmap[d.bioentry_id]
268            d[-1] = t
269            d[-2] = ~numpy.isnan(t)
270            d[-3] = d[-2].all()
271            v = d[-1][d[-2]]
272            v = v - (v.sum() / len(v))
273            v2 = (v * v).sum()
274            d[-4] = (v,v2)
275
276
277
278
279        for i,left in enumerate(data):
280            print i
281            leftidx = left[0]
282            print leftidx
283            leftval = left[-1]
284            r = left[-2]
285            val[leftidx,leftidx] = 1
286            prelr = leftval[r]
287            pn = len(prelr)
288            prelr = prelr - (prelr.sum() / pn)
289            prelr2 = (prelr * prelr).sum()
290            pc = left[-3]
291
292            for right in data[(i+1):]:
293                rightidx = right[0]
294
295               
296                if(right[-3] == False):
297                    r2 = (r & right[-2])
298                    lr = leftval[r2]
299                    n = len(lr)
300                    lr = lr - (lr.sum()/n)
301                    lr2 = (lr * lr).sum()
302                else:
303                    r2 = r
304                    lr = prelr
305                    lr2 = prelr2
306                    n = pn
307
308                if(n < 10):
309                    cor = 0
310                else:
311                    if(pc == False):
312                        rr = right[-1][r2]
313                        rr = rr - (rr.sum()/n)
314                        rr2 = (rr * rr).sum()
315                    else:
316                        rr,rr2 = right[-4]
317
318               
319                    r_num = (lr * rr).sum()
320                    r_den = lr2 * rr2
321                    cor = ((r_num * r_num)/ r_den)
322                val[leftidx,rightidx] = cor
323                val[rightidx,leftidx] = cor
324       
325        val = numpy.sqrt(val)
326   return (gmap,val)
327
328
329def load_ko(d):
330    q = d.dataset[_.name == "p-value"].children_rel()
331    datarankmap = dict([(i.target_dataset_id,i.rank) for i in q])
332    bc = d.bioentry_qualifier.require_type("knockout")()
333    databiomap = dict([(i.dataset_id,i.bioentry_id) for i in bc])
334
335
336   
337    rankbio = numpy.zeros((max(datarankmap.values()),)) - 1;
338    for k,v in databiomap.iteritems():
339        rankbio[datarankmap[k]] = v
340   
341    l = d.dataset[_.name == "p-value"].dataset_id()
342    v = d.bioentry_qualifier_data[_.dataset_id == l]()
343
344    return (v.bioentry_id,rankbio,v[:,-1])
345   
346
347
348
349#FIXME warning: only use for sparse samples. Large samples (e.g. 80%+) on large dataset could take a long time!
350def sample_edges_from_full_graph(nodes,number,exclude=None):
351   source = list(nodes)
352   res = set()
353   nodes_used = set()
354
355   source_max_id = len(source) - 1
356   while(len(res) < number):
357      s = random.randint(0,source_max_id)
358      t = random.randint(0,source_max_id)
359      if(nodes[s] > nodes[t]):
360        s,t = t,s
361      elif(s==t):
362         continue
363      pair = (nodes[s],nodes[t])
364      if((not exclude is None) and pair in exclude):
365         continue
366      else:
367         nodes_used.add(nodes[s])
368         nodes_used.add(nodes[t])
369         res.add(pair)
370     
371   return (res,list(nodes_used))
372
373def all_edges_from_full_graph(nodes,exclude=None):
374   source = list(nodes)
375   res = []
376   for i in range(len(source)):
377      for j in range(i+1,len(source)):
378         s,t = source[i],source[j]
379         if(getattr(s,s._pid) > getattr(t,t._pid)):
380            s,t = t,s
381         pair = (s,t)
382         if(not exclude is None and pair in exclude):
383            continue
384         else:
385            res.append(pair)
386   return set(res)
387
388
389
390   
391#    SELECT bioentry12S.bioentry_id
392#    FROM
393#        ((( taxon_name AS taxon_name10S CROSS JOIN
394#            taxon_bioentry AS taxon_bioentry11S) CROSS JOIN
395#            bioentry AS bioentry12S) CROSS JOIN
396#            term AS term13S)
397#    WHERE (taxon_name10S.name)=('yeast') AND (taxon_name10S.taxon_id)=(taxon_bioentry11S.taxon_id) AND
398#          (taxon_bioentry11S.bioentry_id)=(bioentry12S.bioentry_id) AND (bioentry12S.type_id)=(term13S.term_id) AND
399#          (term13S.name) IN ('CDS'))))
400#    AND ((biorelationship9S.target_bioentry_id) IN ((
401#   
402#    SELECT bioentry12S.bioentry_id
403#    FROM
404#        (((taxon_name AS taxon_name10S CROSS JOIN
405#           taxon_bioentry AS taxon_bioentry11S) CROSS JOIN
406#           bioentry AS bioentry12S) CROSS JOIN
407#           term AS term13S)
408#    WHERE (taxon_name10S.name)=('yeast') AND (taxon_name10S.taxon_id)=(taxon_bioentry11S.taxon_id) AND
409#          (taxon_bioentry11S.bioentry_id)=(bioentry12S.bioentry_id) AND (bioentry12S.type_id)=(term13S.term_id) AND
410#          (term13S.name) IN ('CDS'))))) AND (biorelationship9S.type_id))!=(30613)
411#    return protrel
412
413
Note: See TracBrowser for help on using the repository browser.