Mercurial > repos > andrea.furlani > impc_tools

comparison impc_tool.py @ 0:4357848fb4e6 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 213f6eeb03f96bb13d0ace6e0c87e2562d37f728-dirty
author andrea.furlani
date Thu, 09 Jun 2022 10:51:23 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:4357848fb4e6
1 import sys
2 import requests
3 import pandas as pd
4 import urllib.request as url
5
6 impc_api_url = "https://www.gentar.org/impc-dev-api/"
7 impc_api_search_url = f"{impc_api_url}/genes"
8 impc_api_gene_bundle_url = f"{impc_api_url}/geneBundles"
9
10
11 def stop_err(msg):
12 sys.exit(msg)
13
14
15 def main():
16 inp = str(sys.argv[1])
17 query = str(sys.argv[3])
18
19 try:
20 if query == '7':
21 full_gene_table()
22 sys.exit(0)
23
24 if str(sys.argv[5])=="txt":
25 s = str(sys.argv[6])
26 if s == "t":
27 sep = "\t"
28 elif s == "s":
29 sep = " "
30 elif s in ",;.":
31 sep = s
32 else:
33 sys.exit("Separator not valid, please change it.")
34 inp = pd.read_csv(inp, header=None, delimiter=sep)
35 if len(inp.columns)==1:
36 inp = str(inp[0].values[0]).replace("'","")
37 else:
38 inp = inp.to_string(header=False, index=False).replace(" ",",")
39
40 if query == '8':
41 genes_in_pipeline(inp)
42 sys.exit(0)
43 elif query == '10': # it's here but not totally implemented
44 par_pip_ma(inp)
45 sys.exit(0)
46 elif query == '11': # it's here but not totally implemented
47 par_gen(inp)
48 sys.exit(0)
49 else:
50 tmp = inp.split(",")
51 final_list = []
52 sym_list = []
53 for i in tmp:
54 if 'MGI:' in i or 'MP:' in i:
55 final_list.append(i)
56 else:
57 sym_list.append(i)
58 del(i)
59 if len(sym_list) != 0:
60 sym_list = ",".join(sym_list)
61 biodbnet = f'https://biodbnet.abcc.ncifcrf.gov/webServices/rest.php/biodbnetRestApi.xml?method=db2db&format=row&input=genesymbol&inputValues={sym_list}&outputs=mgiid&taxonId=10090'
62 u = url.urlopen(biodbnet)
63 db = pd.read_xml(u, elems_only=True)
64 if len(db) == 0 and len(final_list) == 0:
65 stop_err("It was not possible to map the input.")
66 else:
67 for i in range(0,len(db)):
68 final_list.append(db['MGIID'][i][4:])
69
70 inp= ",".join(final_list)
71
72
73 if query == '1':
74 get_pheno(inp)
75 sys.exit(0)
76 elif query == '2':
77 get_genes(inp)
78 sys.exit(0)
79 elif query == '3':
80 gene_set(inp)
81 sys.exit(0)
82 elif query == '4':
83 extr_img(inp)
84 sys.exit(0)
85 elif query == '5':
86 parameters(inp)
87 sys.exit(0)
88 elif query == '6':
89 sign_par(inp)
90 sys.exit(0)
91 elif query == '9':
92 sign_mp(inp)
93 sys.exit(0)
94 else:
95 stop_err("Error, non-implemented query selected: " + query)
96 except Exception as ex:
97 stop_err('Error running get_pheno.py:\n' + str(ex))
98
99
100 # 1-Given a gene id, retrieve all the phenotypes related to it (id and name)
101 def get_pheno(inp):
102 head = sys.argv[4]
103 mgi_accession_id = inp
104
105 gene_url = f"{impc_api_search_url}/{mgi_accession_id}"
106 gene_data = requests.get(gene_url).json()
107
108 p_list = []
109 id_list = []
110
111 if gene_data['significantMpTerms'] == None:
112 stop_err("No significant MP terms found for this gene")
113 else:
114 for x in gene_data['significantMpTerms']:
115 p_list.append(x['mpTermId'])
116 id_list.append(x['mpTermName'])
117
118 df = pd.DataFrame()
119 df['MP term name'] = p_list
120 df['MP term id'] = id_list
121
122 if head == 'True':
123 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
124 else:
125 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
126
127
128 # 3-Extract all genes having a particular phenotype or a set of phenotypes (e.g. relevant to a disease)
129 def get_genes(inp):
130 head = sys.argv[4]
131 target_mp_terms = inp
132
133 ## All the data is paginated using the page and size parameters, by default the endpoint returns the first 20 hits
134 gene_by_phenotypes_query = f"{impc_api_search_url}/search/findAllBySignificantMpTermIdsContains?mpTermIds={target_mp_terms}&page=0&size=20"
135 genes_with_clinical_chemistry_phenotypes = requests.get(gene_by_phenotypes_query).json()
136 print(f"Genes with {target_mp_terms}: {genes_with_clinical_chemistry_phenotypes['page']['totalElements']}")
137 list_of_genes = pd.DataFrame(columns=['Gene accession id', 'Gene name', 'Gene bundle url'])
138 acc = []
139 name = []
140 url = []
141
142 for gene in genes_with_clinical_chemistry_phenotypes['_embedded']['genes']:
143 acc.append(gene['mgiAccessionId'])
144 name.append(gene['markerName'])
145 url.append(gene['_links']['geneBundle']['href'])
146
147 list_of_genes['Gene accession id'] = acc
148 list_of_genes['Gene name'] = name
149 list_of_genes['Gene bundle url'] = url
150
151 if head == 'True':
152 list_of_genes.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
153 else:
154 list_of_genes.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
155
156 # 4. Extract all phenotypes which are present in a particular gene set (e.g. genes together in a pathway)
157
158 def gene_set(inp):
159 head = sys.argv[4]
160 target_genes = inp
161
162 genes_in_gene_list_query = f"{impc_api_search_url}/search/findAllByMgiAccessionIdIn?mgiAccessionIds={target_genes}"
163
164 genes_in_gene_list = requests.get(genes_in_gene_list_query).json()
165 list_of_mp_terms_vs_gene_index = {}
166
167 for gene in genes_in_gene_list['_embedded']['genes']:
168 mp_terms = gene['significantMpTerms']
169 gene_acc_id = gene["mgiAccessionId"]
170 if mp_terms is None:
171 continue
172 for mp_term_name in mp_terms:
173 if mp_term_name['mpTermId'] not in list_of_mp_terms_vs_gene_index:
174 list_of_mp_terms_vs_gene_index[mp_term_name['mpTermId']] = {"mp_term": mp_term_name['mpTermId'], "genes": []}
175 list_of_mp_terms_vs_gene_index[mp_term_name['mpTermId']]["genes"].append(gene_acc_id)
176 genes_by_mp_term = list(list_of_mp_terms_vs_gene_index.values())
177
178 df = pd.DataFrame()
179 terms = []
180 genes = []
181 for i in genes_by_mp_term:
182 terms.append(i['mp_term'])
183 genes.append(i['genes'])
184
185 df['mp_term'] = terms
186 df['genes'] = genes
187
188 if head == 'True':
189 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
190 else:
191 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
192
193 # 7. Extract images with a particular phenotype or a set of phenotypes
194
195
196 def extr_img(inp):
197 head = sys.argv[4]
198 target_mp_terms = inp # ['MP:0002110', 'MP:0000559']
199
200 ## All the data is paginated using the page and size parameters, by default the endpoint returns the first 20 hits
201 gene_by_phenotypes_query = f"{impc_api_search_url}/search/findAllBySignificantMpTermIdsContains?mpTermIds={target_mp_terms}&page=0&size=20"
202 genes_with_morphology_mps = requests.get(gene_by_phenotypes_query).json()
203 list_of_gene_bundle_urls = [gene["_links"]["geneBundle"]['href'] for gene in
204 genes_with_morphology_mps['_embedded']['genes']]
205
206 gene_bundles = []
207 for gene_bundle_url in list_of_gene_bundle_urls:
208 gene_bundle = requests.get(gene_bundle_url).json()
209 gene_bundles.append(gene_bundle)
210
211 images_with_morphology_mps = []
212
213 ## Doing just the first 20 and filtering out fields on the images
214 display_fields = ['geneSymbol', 'parameterName', 'biologicalSampleGroup', 'colonyId', 'zygosity', 'sex',
215 'downloadUrl', 'externalSampleId', 'thumbnailUrl']
216
217
218 for gene_bundle in gene_bundles[:20]:
219 if len(gene_bundle) == 4:
220 continue
221 if gene_bundle["geneImages"] is not None:
222 images = gene_bundle["geneImages"]
223 for image in images:
224 display_image = {k: v for k, v in image.items() if k in display_fields}
225 images_with_morphology_mps.append(display_image)
226
227 images_table = []
228 print(f"Images related to phenotype {target_mp_terms}: {len(images_with_morphology_mps)}")
229 ## Displaying just the first 20 images
230 for i in images_with_morphology_mps[:20]:
231 row = [f"<img src='{i['thumbnailUrl']}' />"] + list(i.values())
232 images_table.append(row)
233
234 df = pd.DataFrame()
235 externalSampleId = []
236 geneSymbol = []
237 biologicalSampleGroup = []
238 sex = []
239 colonyId = []
240 zygosity = []
241 parameterName = []
242 downloadUrl = []
243 thumbnailUrl = []
244
245 for i in images_table:
246 externalSampleId.append(i[1])
247 geneSymbol.append(i[2])
248 biologicalSampleGroup.append(i[3])
249 sex.append(i[4])
250 colonyId.append(i[5])
251 zygosity.append(i[6])
252 parameterName.append(i[7])
253 downloadUrl.append(i[8])
254 thumbnailUrl.append(i[9])
255
256 df['externalSampleId'] = externalSampleId
257 df['geneSymbol'] = geneSymbol
258 df['biologicalSampleGroup'] = biologicalSampleGroup
259 df['sex'] = sex
260 df['colonyId'] = colonyId
261 df['zygosity'] = zygosity
262 df['parameterName'] = parameterName
263 df['downloadUrl'] = downloadUrl
264 df['thumbnailUrl'] = thumbnailUrl
265
266 if head == 'True':
267 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
268 else:
269 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
270
271 # 11- Which parameters have been measured for a particular knockout EASY
272
273
274 def parameters(inp):
275 head = sys.argv[4]
276 knockout = inp # "MGI:104636"
277 gene_info = requests.get(impc_api_search_url + "/" + knockout).json()
278
279 if gene_info['phenotypingDataAvailable']:
280 geneBundle = requests.get(gene_info['_links']['geneBundle']['href']).json()
281 gen_imgs = geneBundle['geneImages']
282 par_list = []
283 l = {}
284 for i in gen_imgs:
285 l = {"Parameter Name": i['parameterName']}
286 if l not in par_list:
287 par_list.append(l)
288 df = pd.DataFrame()
289 l = []
290
291 for i in par_list:
292 l.append(i['Parameter Name'])
293
294 df['Parameter'] = l
295 if head == 'True':
296 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
297 else:
298 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
299
300 else:
301 stop_err("No parameters available for this knockout gene")
302
303
304 # 12- Which parameters identified a significant finding for a particular knockout line (colony) EASY
305 def sign_par(inp):
306 head = sys.argv[4]
307 knockout = inp # "MGI:104636"
308
309 gene_info = requests.get(f"{impc_api_url}statisticalResults/search/findAllByMarkerAccessionIdIsAndSignificantTrue?mgiAccessionId=" + knockout).json()
310 gene_stats = gene_info['_embedded']['statisticalResults']
311
312 if len(gene_stats) == 0:
313 stop_err("No statistically relevant parameters found for this knockout gene")
314 else:
315 df = pd.DataFrame()
316 n = []
317 p = []
318 for g in gene_stats:
319 n.append(g['parameterName'])
320 p.append(g['pvalue'])
321
322 df['Parameter name'] = n
323 df['p-value'] = p
324 if head == 'True':
325 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
326 else:
327 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
328
329
330 # 13- List of genes names and ID measured in a pipeline
331 def genes_in_pipeline(inp):
332 head = sys.argv[4]
333 pip = inp
334
335 g_in_p_query = f"{impc_api_search_url}/search/findAllByTestedPipelineId?pipelineId={pip}&page=0&size=1000"
336 genes_in_pip = requests.get(g_in_p_query).json()
337 pages = genes_in_pip['page']['totalPages']
338 max_elem = genes_in_pip['page']['totalElements']
339
340 print(f"Genes with {pip}: {genes_in_pip['page']['totalElements']}")
341 d ={ }
342 list_d = []
343 list_of_genes = pd.DataFrame(columns=['Gene accession id', 'Gene name'])
344 acc = []
345 name = []
346
347 if max_elem > 1000:
348 g_in_p_query = genes_in_pip['_embedded']['genes']
349 for i in range(1,pages):
350 gl = requests.get(f'{impc_api_search_url}/search/findAllByTestedPipelineId?pipelineId={pip}&page={i}&size=1000').json()['_embedded']['genes']
351 g_in_p_query += gl
352 else:
353 g_in_p_query = genes_in_pip['_embedded']['genes']
354
355 for g in g_in_p_query:
356 d = {"Gene Accession ID": g['mgiAccessionId'], "Gene Name": g['markerName']}
357 list_d.append(d)
358
359 for i in list_d:
360 acc.append(i['Gene Accession ID'])
361 name.append(i['Gene Name'])
362
363 list_of_genes['Gene accession id'] = acc
364 list_of_genes['Gene name'] = name
365
366 if head == 'True':
367 list_of_genes.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
368 else:
369 list_of_genes.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
370
371
372 # 14- Extract all genes and corresponding phenotypes related to a particular organ system(eg: significatMPTerm)
373 def sign_mp(inp):
374 head = sys.argv[4]
375 mp_term = inp # ['MP:0005391']
376
377 gene_by_mpterm_query = f"{impc_api_search_url}/search/findAllBySignificantMpTermIdsContains?mpTermIds={mp_term}&size=1000"
378 genes_with_mpterm = requests.get(gene_by_mpterm_query).json()
379
380 pages = genes_with_mpterm['page']['totalPages']
381 genes_info = genes_with_mpterm['_embedded']['genes']
382
383 for pn in range(1,pages):
384 pq = f"{impc_api_search_url}/search/findAllBySignificantMpTermIdsContains?mpTermIds={mp_term}&page={pn}&size=1000"
385 g = requests.get(pq).json()['_embedded']['genes']
386 genes_info += g
387
388 list_d=[]
389 d={}
390 for g in genes_info:
391 names=[]
392 ids=[]
393 for s in g['significantMpTerms']:
394 names.append(s['mpTermName'])
395 ids.append(s['mpTermId'])
396 d={'Gene':g['mgiAccessionId'], 'mpTermId': ids, 'mpTermName':names}
397 list_d.append(d)
398
399
400 g = []
401 ids = []
402 names = []
403 for i in list_d:
404 g.append(i['Gene'])
405 ids.append(i['mpTermId'])
406 names.append(i['mpTermName'])
407
408 df = pd.DataFrame()
409 df['Gene Id']=g
410 df['Significant MP terms Ids']=ids
411 df['Significant MP terms Names']=names
412
413 if head == 'True':
414 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
415 else:
416 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
417
418
419 # 16- Full table of genes and all identified phenotypes
420
421 def full_gene_table():
422 head = sys.argv[4]
423 gene_list = requests.get(impc_api_search_url + '?page=0&size=1000').json()
424 pages = gene_list['page']['totalPages']
425 genes_info = gene_list['_embedded']['genes']
426
427 for pn in range(1,pages):
428 gp = requests.get(impc_api_search_url + f'?page={pn}&size=1000').json()['_embedded']['genes']
429 genes_info += gp
430
431 d = {}
432 list_d=[]
433
434 for i in genes_info:
435 l = []
436 if i['significantMpTerms'] is None:
437 d={"Gene": i['mgiAccessionId'], "Identified phenotypes": "None"}
438 else:
439 d = {"Gene": i['mgiAccessionId'], "Identified phenotypes": [sub['mpTermId'] for sub in i['significantMpTerms']]}
440 list_d.append(d)
441
442 df = pd.DataFrame()
443 g = []
444 p = []
445 for i in list_d:
446 g.append(i['Gene'])
447 p.append(i['Identified phenotypes'])
448
449 df['MGI id'] = g
450 df['MP term list'] = p
451
452 for i in range(0, len(df)):
453 if df['MP term list'][i] != "None":
454 df['MP term list'][i] = str(df['MP term list'][i])[1:-1].replace("'", "")
455
456 if str(sys.argv[1]) == 'True':
457 if head == 'True':
458 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
459 else:
460 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
461 else:
462 df = df[df['MP term list'] != "None"]
463 df.reset_index(drop=True, inplace=True)
464 if head == 'True':
465 df.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
466 else:
467 df.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
468
469 # Old method, chech which is faster
470 # max_elem = gene_list['page']['totalElements']
471 # d = {}
472 # list_d = []
473 # for i in range(0, pages):
474 # gl = requests.get(impc_api_search_url + '?page=' + str(i) + '&size=' + str(max_elem)).json()
475 # for g in gl['_embedded']['genes']:
476 # if g['significantMpTerms'] is None:
477 # d = {"Gene": g['mgiAccessionId'], "Identified phenotypes": "None"}
478 # else:
479 # d = {"Gene": g['mgiAccessionId'], "Identified phenotypes": [ sub['mpTermId'] for sub in g['significantMpTerms'] ]}
480 # list_d.append(d)
481
482
483
484
485 # 18- Extract measurements and analysis for a parameter or pipeline
486
487 def par_pip_ma(inp):
488 head = sys.argv[4]
489 id = inp
490
491 if id[0:4] == "IMPC":
492 par = True
493 ma_query = f"{impc_api_search_url}/search/findAllByTestedParameterId?parameterId={id}&page=0&size=1000"
494 else:
495 ma_query = f"{impc_api_search_url}/search/findAllByTestedPipelineId?pipelineId={id}&page=0&size=1000"
496 par = False
497
498 ma_in_pip = requests.get(ma_query).json()
499 pages = ma_in_pip['page']['totalPages']
500 max_elem = ma_in_pip['page']['totalElements']
501
502 print(f"Genes with {id}: {ma_in_pip['page']['totalElements']}")
503 d = {}
504 list_d = []
505 list_of_genes = pd.DataFrame(columns=['Measurements', 'Analysis'])
506 mes = []
507 an = []
508
509 if max_elem > 1000:
510
511 ma_in_pip = ma_in_pip['_embedded']['genes']
512 for pn in range(1, pages):
513 if par:
514 pip = requests.get(f"{impc_api_search_url}/search/findAllByTestedParameterId?parameterId={id}&page={pn}&size=1000").json()['_embedded']['genes']
515 else:
516 pip = requests.get(f"{impc_api_search_url}/search/findAllByTestedPipelineId?pipelineId={id}&page={pn}&size=1000").json()['_embedded']['genes']
517 ma_in_pip += pip
518
519 else:
520 ma_in_pip = ma_in_pip['_embedded']['genes']
521
522 for g in ma_in_pip:
523 d = {"Measurements": g[''], "Analysis": g['']}
524 list_d.append(d)
525
526 for i in list_d:
527 mes.append(i[''])
528 an.append(i[''])
529
530 list_of_genes['Analysis'] = an
531 list_of_genes['Measurements'] = mes
532
533 if head == 'True':
534 list_of_genes.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
535 else:
536 list_of_genes.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
537
538
539 # 19- Get all genes and measured values for a particular parameter
540 def par_gen(inp):
541 head = sys.argv[4]
542 id = inp
543
544 pa_query = f"{impc_api_search_url}/search/findAllByTestedParameterId?parameterId={id}&page=0&size=1000"
545
546 gm_par = requests.get(pa_query).json()
547 pages = gm_par['page']['totalPages']
548 max_elem = gm_par['page']['totalElements']
549
550 print(f"Genes with {id}: {gm_par['page']['totalElements']}")
551 d = {}
552 list_d = []
553 list_of_genes = pd.DataFrame(columns=['Genes', 'Measured Values'])
554 gen = []
555 mes = []
556
557 if max_elem > 1000:
558
559 gm_par = gm_par['_embedded']['genes']
560
561 for pn in range(1, pages):
562 pip = requests.get(f"{impc_api_search_url}/search/findAllByTestedParameterId?parameterId={id}&page={pn}&size=1000").json()['_embedded']['genes']
563 gm_par += pip
564
565 else:
566 gm_par = gm_par['_embedded']['genes']
567
568
569 for g in gm_par:
570 d = {"Genes": g['mgiAccessionId'], "Measured Values": g['']}
571 list_d.append(d)
572
573 for i in list_d:
574 gen.append(i['Genes'])
575 mes.append(i['Measured Values'])
576
577 list_of_genes['Genes'] = gen
578 list_of_genes['Measured Values'] = mes
579
580 if head == 'True':
581 list_of_genes.to_csv(sys.argv[2], header=True, index=False, sep="\t", index_label=False)
582 else:
583 list_of_genes.to_csv(sys.argv[2], header=False, index=False, sep="\t", index_label=False)
584
585
586 if __name__ == "__main__":
587 main()