Mercurial > repos > bimib > marea

comparison Marea/ras_generator.py @ 47:3af9d394367c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bimib
date Wed, 19 Feb 2020 10:44:52 -0500
parents 5d5d01ef1d68
children 3b0e71e28c0b
comparison
equal deleted inserted replaced
46:5d5d01ef1d68 47:3af9d394367c
1 from __future__ import division 1 from __future__ import division
2 import sys 2 import sys
3 import pandas as pd 3 import pandas as pd
4 import itertools as it
5 import scipy.stats as st
6 import collections 4 import collections
7 import lxml.etree as ET
8 import pickle as pk 5 import pickle as pk
9 import math 6 import math
10 import os
11 import argparse 7 import argparse
12 from svglib.svglib import svg2rlg
13 from reportlab.graphics import renderPDF
14 8
15 ########################## argparse ########################################## 9 ########################## argparse ##########################################
16 10
17 def process_args(args): 11 def process_args(args):
18 parser = argparse.ArgumentParser(usage = '%(prog)s [options]', 12 parser = argparse.ArgumentParser(usage = '%(prog)s [options]',
24 choices = ['HMRcore', 'Recon', 'Custom'], 18 choices = ['HMRcore', 'Recon', 'Custom'],
25 help = 'chose which type of dataset you want use') 19 help = 'chose which type of dataset you want use')
26 parser.add_argument('-cr', '--custom', 20 parser.add_argument('-cr', '--custom',
27 type = str, 21 type = str,
28 help='your dataset if you want custom rules') 22 help='your dataset if you want custom rules')
29 parser.add_argument('-na', '--names',
30 type = str,
31 nargs = '+',
32 help = 'input names')
33 parser.add_argument('-n', '--none', 23 parser.add_argument('-n', '--none',
34 type = str, 24 type = str,
35 default = 'true', 25 default = 'true',
36 choices = ['true', 'false'], 26 choices = ['true', 'false'],
37 help = 'compute Nan values') 27 help = 'compute Nan values')
38 parser.add_argument('-pv' ,'--pValue',
39 type = float,
40 default = 0.05,
41 help = 'P-Value threshold (default: %(default)s)')
42 parser.add_argument('-fc', '--fChange',
43 type = float,
44 default = 1.5,
45 help = 'Fold-Change threshold (default: %(default)s)')
46 parser.add_argument('-td', '--tool_dir', 28 parser.add_argument('-td', '--tool_dir',
47 type = str, 29 type = str,
48 required = True, 30 required = True,
49 help = 'your tool directory') 31 help = 'your tool directory')
50 parser.add_argument('-op', '--option',
51 type = str,
52 choices = ['datasets', 'dataset_class', 'datasets_rasonly'],
53 help='dataset or dataset and class')
54 parser.add_argument('-ol', '--out_log', 32 parser.add_argument('-ol', '--out_log',
55 help = "Output log") 33 help = "Output log")
56 parser.add_argument('-ids', '--input_datas', 34 parser.add_argument('-id', '--input',
57 type = str,
58 nargs = '+',
59 help = 'input datasets')
60 parser.add_argument('-id', '--input_data',
61 type = str, 35 type = str,
62 help = 'input dataset') 36 help = 'input dataset')
63 parser.add_argument('-ic', '--input_class', 37 parser.add_argument('-ra', '--ras_output',
64 type = str,
65 help = 'sample group specification')
66 parser.add_argument('-cm', '--custom_map',
67 type = str,
68 help = 'custom map')
69 parser.add_argument('-yn', '--yes_no',
70 type = str, 38 type = str,
71 choices = ['yes', 'no'], 39 required = True,
72 help = 'if make or not custom map') 40 help = 'ras output')
73 parser.add_argument('-gs', '--generate_svg', 41
74 type = str,
75 default = 'true',
76 choices = ['true', 'false'],
77 help = 'generate svg map')
78 parser.add_argument('-gp', '--generate_pdf',
79 type = str,
80 default = 'true',
81 choices = ['true', 'false'],
82 help = 'generate pdf map')
83 parser.add_argument('-gr', '--generate_ras',
84 type = str,
85 default = 'true',
86 choices = ['true', 'false'],
87 help = 'generate reaction activity score')
88 parser.add_argument('-sr', '--single_ras_file',
89 type = str,
90 help = 'file that will contain ras')
91
92 args = parser.parse_args() 42 args = parser.parse_args()
93 return args 43 return args
94 44
95 ########################### warning ########################################### 45 ########################### warning ###########################################
96 46
294 return False 244 return False
295 l = l[3:] 245 l = l[3:]
296 else: 246 else:
297 return False 247 return False
298 return ris 248 return ris
299
300 ############################ map_methods ######################################
301
302 def fold_change(avg1, avg2):
303 if avg1 == 0 and avg2 == 0:
304 return 0
305 elif avg1 == 0:
306 return '-INF'
307 elif avg2 == 0:
308 return 'INF'
309 else:
310 return math.log(avg1 / avg2, 2)
311
312 def fix_style(l, col, width, dash):
313 tmp = l.split(';')
314 flag_col = False
315 flag_width = False
316 flag_dash = False
317 for i in range(len(tmp)):
318 if tmp[i].startswith('stroke:'):
319 tmp[i] = 'stroke:' + col
320 flag_col = True
321 if tmp[i].startswith('stroke-width:'):
322 tmp[i] = 'stroke-width:' + width
323 flag_width = True
324 if tmp[i].startswith('stroke-dasharray:'):
325 tmp[i] = 'stroke-dasharray:' + dash
326 flag_dash = True
327 if not flag_col:
328 tmp.append('stroke:' + col)
329 if not flag_width:
330 tmp.append('stroke-width:' + width)
331 if not flag_dash:
332 tmp.append('stroke-dasharray:' + dash)
333 return ';'.join(tmp)
334
335 def fix_map(d, core_map, threshold_P_V, threshold_F_C, max_F_C):
336 maxT = 12
337 minT = 2
338 grey = '#BEBEBE'
339 blue = '#0000FF'
340 red = '#E41A1C'
341 for el in core_map.iter():
342 el_id = str(el.get('id'))
343 if el_id.startswith('R_'):
344 tmp = d.get(el_id[2:])
345 if tmp != None:
346 p_val = tmp[0]
347 f_c = tmp[1]
348 if p_val < threshold_P_V:
349 if not isinstance(f_c, str):
350 if abs(f_c) < math.log(threshold_F_C, 2):
351 col = grey
352 width = str(minT)
353 else:
354 if f_c < 0:
355 col = blue
356 elif f_c > 0:
357 col = red
358 width = str(max((abs(f_c) * maxT) / max_F_C, minT))
359 else:
360 if f_c == '-INF':
361 col = blue
362 elif f_c == 'INF':
363 col = red
364 width = str(maxT)
365 dash = 'none'
366 else:
367 dash = '5,5'
368 col = grey
369 width = str(minT)
370 el.set('style', fix_style(el.get('style'), col, width, dash))
371 return core_map
372 249
373 ############################ make recon ####################################### 250 ############################ make recon #######################################
374 251
375 def check_and_doWord(l): 252 def check_and_doWord(l):
376 tmp = [] 253 tmp = []
613 490
614 ############################ resolve ########################################## 491 ############################ resolve ##########################################
615 492
616 def resolve(genes, rules, ids, resolve_none, name): 493 def resolve(genes, rules, ids, resolve_none, name):
617 resolve_rules = {} 494 resolve_rules = {}
618 names_array = []
619 not_found = [] 495 not_found = []
620 flag = False 496 flag = False
621 for key, value in genes.items(): 497 for key, value in genes.items():
622 tmp_resolve = [] 498 tmp_resolve = []
623 for i in range(len(rules)): 499 for i in range(len(rules)):
650 if not pd.isnull(classe): 526 if not pd.isnull(classe):
651 l = [] 527 l = []
652 for j in range(i, len(classes)): 528 for j in range(i, len(classes)):
653 if classes.iloc[j, 1] == classe: 529 if classes.iloc[j, 1] == classe:
654 pat_id = classes.iloc[j, 0] 530 pat_id = classes.iloc[j, 0]
655 tmp = resolve_rules.get(pat_id, None)
656 if tmp != None: 531 if tmp != None:
657 l.append(tmp) 532 l.append(tmp)
658 classes.iloc[j, 1] = None 533 classes.iloc[j, 1] = None
659 if l: 534 if l:
660 class_pat[classe] = list(map(list, zip(*l))) 535 class_pat[classe] = list(map(list, zip(*l)))
663 ', the class has been disregarded\n') 538 ', the class has been disregarded\n')
664 return class_pat 539 return class_pat
665 540
666 ############################ create_ras ####################################### 541 ############################ create_ras #######################################
667 542
668 def create_ras (resolve_rules, dataset_name, single_ras, rules, ids): 543 def create_ras (resolve_rules, dataset_name, rules, ids, file):
669 544
670 if resolve_rules == None: 545 if resolve_rules == None:
671 warning("Couldn't generate RAS for current dataset: " + dataset_name) 546 warning("Couldn't generate RAS for current dataset: " + dataset_name)
672 547
673 for geni in resolve_rules.values(): 548 for geni in resolve_rules.values():
678 output_ras = pd.DataFrame.from_dict(resolve_rules) 553 output_ras = pd.DataFrame.from_dict(resolve_rules)
679 554
680 output_ras.insert(0, 'Reactions', ids) 555 output_ras.insert(0, 'Reactions', ids)
681 output_to_csv = pd.DataFrame.to_csv(output_ras, sep = '\t', index = False) 556 output_to_csv = pd.DataFrame.to_csv(output_ras, sep = '\t', index = False)
682 557
683 if (single_ras): 558 text_file = open(file, "w")
684 args = process_args(sys.argv)
685 text_file = open(args.single_ras_file, "w")
686 else:
687 text_file = open("ras/Reaction_Activity_Score_Of_" + dataset_name + ".tsv", "w")
688 559
689 text_file.write(output_to_csv) 560 text_file.write(output_to_csv)
690 text_file.close() 561 text_file.close()
691 562
692 ############################ map ##############################################
693
694 def maps(core_map, class_pat, ids, threshold_P_V, threshold_F_C, create_svg, create_pdf):
695 args = process_args(sys.argv)
696 if (not class_pat) or (len(class_pat.keys()) < 2):
697 sys.exit('Execution aborted: classes provided for comparisons are ' +
698 'less than two\n')
699 for i, j in it.combinations(class_pat.keys(), 2):
700 tmp = {}
701 count = 0
702 max_F_C = 0
703 for l1, l2 in zip(class_pat.get(i), class_pat.get(j)):
704 try:
705 stat_D, p_value = st.ks_2samp(l1, l2)
706 avg = fold_change(sum(l1) / len(l1), sum(l2) / len(l2))
707 if not isinstance(avg, str):
708 if max_F_C < abs(avg):
709 max_F_C = abs(avg)
710 tmp[ids[count]] = [float(p_value), avg]
711 count += 1
712 except (TypeError, ZeroDivisionError):
713 count += 1
714 tab = 'result/' + i + '_vs_' + j + ' (Tabular Result).tsv'
715 tmp_csv = pd.DataFrame.from_dict(tmp, orient = "index")
716 tmp_csv = tmp_csv.reset_index()
717 header = ['ids', 'P_Value', 'Log2(fold change)']
718 tmp_csv.to_csv(tab, sep = '\t', index = False, header = header)
719
720 if create_svg or create_pdf:
721 if args.rules_selector == 'HMRcore' or (args.rules_selector == 'Custom'
722 and args.yes_no == 'yes'):
723 fix_map(tmp, core_map, threshold_P_V, threshold_F_C, max_F_C)
724 file_svg = 'result/' + i + '_vs_' + j + ' (SVG Map).svg'
725 with open(file_svg, 'wb') as new_map:
726 new_map.write(ET.tostring(core_map))
727
728
729 if create_pdf:
730 file_pdf = 'result/' + i + '_vs_' + j + ' (PDF Map).pdf'
731 renderPDF.drawToFile(svg2rlg(file_svg), file_pdf)
732
733 if not create_svg:
734 #Ho utilizzato il file svg per generare il pdf,
735 #ma l'utente non ne ha richiesto il ritorno, quindi
736 #lo elimino
737 os.remove('result/' + i + '_vs_' + j + ' (SVG Map).svg')
738
739 return None
740
741 ############################ MAIN ############################################# 563 ############################ MAIN #############################################
742 564
743 def main(): 565 def main():
744 args = process_args(sys.argv) 566 args = process_args(sys.argv)
745 567
746 create_svg = check_bool(args.generate_svg)
747 create_pdf = check_bool(args.generate_pdf)
748 generate_ras = check_bool(args.generate_ras)
749
750 os.makedirs('result')
751
752 if generate_ras:
753 os.makedirs('ras')
754
755 if args.rules_selector == 'HMRcore': 568 if args.rules_selector == 'HMRcore':
756 recon = pk.load(open(args.tool_dir + '/local/HMRcore_rules.p', 'rb')) 569 recon = pk.load(open(args.tool_dir + '/local/HMRcore_rules.p', 'rb'))
757 elif args.rules_selector == 'Recon': 570 elif args.rules_selector == 'Recon':
758 recon = pk.load(open(args.tool_dir + '/local/Recon_rules.p', 'rb')) 571 recon = pk.load(open(args.tool_dir + '/local/Recon_rules.p', 'rb'))
759 elif args.rules_selector == 'Custom': 572 elif args.rules_selector == 'Custom':
760 ids, rules, gene_in_rule = make_recon(args.custom) 573 ids, rules, gene_in_rule = make_recon(args.custom)
761 574
762 resolve_none = check_bool(args.none) 575 resolve_none = check_bool(args.none)
763 576
764 class_pat = {} 577
765 578 name = "RAS Dataset"
766 if args.option == 'datasets_rasonly': 579 dataset = read_dataset(args.input, "dataset")
767 name = "RAS Dataset" 580
768 dataset = read_dataset(args.input_datas[0],"dataset") 581 dataset.iloc[:, 0] = (dataset.iloc[:, 0]).astype(str)
769 582
770 dataset.iloc[:, 0] = (dataset.iloc[:, 0]).astype(str) 583 type_gene = gene_type(dataset.iloc[0, 0], name)
771
772 type_gene = gene_type(dataset.iloc[0, 0], name)
773
774 if args.rules_selector != 'Custom':
775 genes = data_gene(dataset, type_gene, name, None)
776 ids, rules = load_id_rules(recon.get(type_gene))
777 elif args.rules_selector == 'Custom':
778 genes = data_gene(dataset, type_gene, name, gene_in_rule)
779 584
780 resolve_rules, err = resolve(genes, rules, ids, resolve_none, name) 585 if args.rules_selector != 'Custom':
781 586 genes = data_gene(dataset, type_gene, name, None)
782 create_ras(resolve_rules, name, True, rules, ids) 587 ids, rules = load_id_rules(recon.get(type_gene))
783 588 elif args.rules_selector == 'Custom':
784 if err != None and err: 589 genes = data_gene(dataset, type_gene, name, gene_in_rule)
785 warning('Warning: gene\n' + str(err) + '\nnot found in class ' 590
786 + name + ', the expression level for this gene ' + 591 resolve_rules, err = resolve(genes, rules, ids, resolve_none, name)
787 'will be considered NaN\n') 592
788 593 create_ras(resolve_rules, name, rules, ids, args.ras_output)
789 print('execution succeded') 594
790 return None 595 if err != None and err:
791 596 warning('Warning: gene\n' + str(err) + '\nnot found in class '
792 597 + name + ', the expression level for this gene ' +
793 elif args.option == 'datasets': 598 'will be considered NaN\n')
794 num = 1 599
795 for i, j in zip(args.input_datas, args.names): 600
796
797 name = name_dataset(j, num)
798 dataset = read_dataset(i, name)
799
800 dataset.iloc[:, 0] = (dataset.iloc[:, 0]).astype(str)
801
802 type_gene = gene_type(dataset.iloc[0, 0], name)
803
804 if args.rules_selector != 'Custom':
805 genes = data_gene(dataset, type_gene, name, None)
806 ids, rules = load_id_rules(recon.get(type_gene))
807 elif args.rules_selector == 'Custom':
808 genes = data_gene(dataset, type_gene, name, gene_in_rule)
809
810
811 resolve_rules, err = resolve(genes, rules, ids, resolve_none, name)
812
813 if generate_ras:
814 create_ras(resolve_rules, name, False, rules, ids)
815
816 if err != None and err:
817 warning('Warning: gene\n' + str(err) + '\nnot found in class '
818 + name + ', the expression level for this gene ' +
819 'will be considered NaN\n')
820 if resolve_rules != None:
821 class_pat[name] = list(map(list, zip(*resolve_rules.values())))
822 num += 1
823 elif args.option == 'dataset_class':
824 name = 'RNAseq'
825 dataset = read_dataset(args.input_data, name)
826 dataset.iloc[:, 0] = (dataset.iloc[:, 0]).astype(str)
827 type_gene = gene_type(dataset.iloc[0, 0], name)
828 classes = read_dataset(args.input_class, 'class')
829 if not len(classes.columns) == 2:
830 warning('Warning: more than 2 columns in class file. Extra' +
831 'columns have been disregarded\n')
832 classes = classes.astype(str)
833 if args.rules_selector != 'Custom':
834 genes = data_gene(dataset, type_gene, name, None)
835 ids, rules = load_id_rules(recon.get(type_gene))
836 elif args.rules_selector == 'Custom':
837 genes = data_gene(dataset, type_gene, name, gene_in_rule)
838 resolve_rules, err = resolve(genes, rules, ids, resolve_none, name)
839 if err != None and err:
840 warning('Warning: gene\n'+str(err)+'\nnot found in class '
841 + name + ', the expression level for this gene ' +
842 'will be considered NaN\n')
843 if resolve_rules != None:
844 class_pat = split_class(classes, resolve_rules)
845 if generate_ras:
846 create_ras(resolve_rules, name, False, rules, ids)
847
848
849 if args.rules_selector == 'Custom':
850 if args.yes_no == 'yes':
851 try:
852 core_map = ET.parse(args.custom_map)
853 except (ET.XMLSyntaxError, ET.XMLSchemaParseError):
854 sys.exit('Execution aborted: custom map in wrong format')
855 elif args.yes_no == 'no':
856 core_map = ET.parse(args.tool_dir + '/local/HMRcoreMap.svg')
857 else:
858 core_map = ET.parse(args.tool_dir+'/local/HMRcoreMap.svg')
859
860 maps(core_map, class_pat, ids, args.pValue, args.fChange, create_svg, create_pdf)
861
862 print('Execution succeded') 601 print('Execution succeded')
863 602
864 return None 603 return None
865 604
866 ############################################################################### 605 ###############################################################################