Mercurial > repos > bimib > marea

comparison marea-1.0.1/marea.py @ 15:d0e7f14b773f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Upload 1.0.1
author bimib
date Tue, 01 Oct 2019 06:03:12 -0400
parents
children
comparison
equal deleted inserted replaced
14:1a0c8c2780f2 15:d0e7f14b773f
1 from __future__ import division
2 import sys
3 import pandas as pd
4 import itertools as it
5 import scipy.stats as st
6 import collections
7 import lxml.etree as ET
8 import shutil
9 import pickle as pk
10 import math
11 import os
12 import argparse
13 from svglib.svglib import svg2rlg
14 from reportlab.graphics import renderPDF
15
16 ########################## argparse ##########################################
17
18 def process_args(args):
19 parser = argparse.ArgumentParser(usage = '%(prog)s [options]',
20 description = 'process some value\'s'+
21 ' genes to create a comparison\'s map.')
22 parser.add_argument('-rs', '--rules_selector',
23 type = str,
24 default = 'HMRcore',
25 choices = ['HMRcore', 'Recon', 'Custom'],
26 help = 'chose which type of dataset you want use')
27 parser.add_argument('-cr', '--custom',
28 type = str,
29 help='your dataset if you want custom rules')
30 parser.add_argument('-na', '--names',
31 type = str,
32 nargs = '+',
33 help = 'input names')
34 parser.add_argument('-n', '--none',
35 type = str,
36 default = 'true',
37 choices = ['true', 'false'],
38 help = 'compute Nan values')
39 parser.add_argument('-pv' ,'--pValue',
40 type = float,
41 default = 0.05,
42 help = 'P-Value threshold (default: %(default)s)')
43 parser.add_argument('-fc', '--fChange',
44 type = float,
45 default = 1.5,
46 help = 'Fold-Change threshold (default: %(default)s)')
47 parser.add_argument('-td', '--tool_dir',
48 type = str,
49 required = True,
50 help = 'your tool directory')
51 parser.add_argument('-op', '--option',
52 type = str,
53 choices = ['datasets', 'dataset_class'],
54 help='dataset or dataset and class')
55 parser.add_argument('-ol', '--out_log',
56 help = "Output log")
57 parser.add_argument('-ids', '--input_datas',
58 type = str,
59 nargs = '+',
60 help = 'input datasets')
61 parser.add_argument('-id', '--input_data',
62 type = str,
63 help = 'input dataset')
64 parser.add_argument('-ic', '--input_class',
65 type = str,
66 help = 'sample group specification')
67 parser.add_argument('-cm', '--custom_map',
68 type = str,
69 help = 'custom map')
70 parser.add_argument('-yn', '--yes_no',
71 type = str,
72 choices = ['yes', 'no'],
73 help = 'if make or not custom map')
74 parser.add_argument('-gs', '--generate_svg',
75 type = str,
76 default = 'true',
77 choices = ['true', 'false'],
78 help = 'generate svg map')
79 parser.add_argument('-gp', '--generate_pdf',
80 type = str,
81 default = 'true',
82 choices = ['true', 'false'],
83 help = 'generate pdf map')
84 parser.add_argument('-gr', '--generate_ras',
85 type = str,
86 default = 'true',
87 choices = ['true', 'false'],
88 help = 'generate reaction activity score')
89 args = parser.parse_args()
90 return args
91
92 ########################### warning ###########################################
93
94 def warning(s):
95 args = process_args(sys.argv)
96 with open(args.out_log, 'a') as log:
97 log.write(s)
98
99 ############################ dataset input ####################################
100
101 def read_dataset(data, name):
102 try:
103 dataset = pd.read_csv(data, sep = '\t', header = 0, engine='python')
104 except pd.errors.EmptyDataError:
105 sys.exit('Execution aborted: wrong format of ' + name + '\n')
106 if len(dataset.columns) < 2:
107 sys.exit('Execution aborted: wrong format of ' + name + '\n')
108 return dataset
109
110 ############################ dataset name #####################################
111
112 def name_dataset(name_data, count):
113 if str(name_data) == 'Dataset':
114 return str(name_data) + '_' + str(count)
115 else:
116 return str(name_data)
117
118 ############################ load id e rules ##################################
119
120 def load_id_rules(reactions):
121 ids, rules = [], []
122 for key, value in reactions.items():
123 ids.append(key)
124 rules.append(value)
125 return (ids, rules)
126
127 ############################ check_methods ####################################
128
129 def gene_type(l, name):
130 if check_hgnc(l):
131 return 'hugo_id'
132 elif check_ensembl(l):
133 return 'ensembl_gene_id'
134 elif check_symbol(l):
135 return 'symbol'
136 elif check_entrez(l):
137 return 'entrez_id'
138 else:
139 sys.exit('Execution aborted:\n' +
140 'gene ID type in ' + name + ' not supported. Supported ID'+
141 'types are: HUGO ID, Ensemble ID, HUGO symbol, Entrez ID\n')
142
143 def check_hgnc(l):
144 if len(l) > 5:
145 if (l.upper()).startswith('HGNC:'):
146 return l[5:].isdigit()
147 else:
148 return False
149 else:
150 return False
151
152 def check_ensembl(l):
153 if len(l) == 15:
154 if (l.upper()).startswith('ENS'):
155 return l[4:].isdigit()
156 else:
157 return False
158 else:
159 return False
160
161 def check_symbol(l):
162 if len(l) > 0:
163 if l[0].isalpha() and l[1:].isalnum():
164 return True
165 else:
166 return False
167 else:
168 return False
169
170 def check_entrez(l):
171 if len(l) > 0:
172 return l.isdigit()
173 else:
174 return False
175
176 def check_bool(b):
177 if b == 'true':
178 return True
179 elif b == 'false':
180 return False
181
182 ############################ resolve_methods ##################################
183
184 def replace_gene_value(l, d):
185 tmp = []
186 err = []
187 while l:
188 if isinstance(l[0], list):
189 tmp_rules, tmp_err = replace_gene_value(l[0], d)
190 tmp.append(tmp_rules)
191 err.extend(tmp_err)
192 else:
193 value = replace_gene(l[0], d)
194 tmp.append(value)
195 if value == None:
196 err.append(l[0])
197 l = l[1:]
198 return (tmp, err)
199
200 def replace_gene(l, d):
201 if l =='and' or l == 'or':
202 return l
203 else:
204 value = d.get(l, None)
205 if not(value == None or isinstance(value, (int, float))):
206 sys.exit('Execution aborted: ' + value + ' value not valid\n')
207 return value
208
209 def computes(val1, op, val2, cn):
210 if val1 != None and val2 != None:
211 if op == 'and':
212 return min(val1, val2)
213 else:
214 return val1 + val2
215 elif op == 'and':
216 if cn is True:
217 if val1 != None:
218 return val1
219 elif val2 != None:
220 return val2
221 else:
222 return None
223 else:
224 return None
225 else:
226 if val1 != None:
227 return val1
228 elif val2 != None:
229 return val2
230 else:
231 return None
232
233 def control(ris, l, cn):
234 if len(l) == 1:
235 if isinstance(l[0], (float, int)) or l[0] == None:
236 return l[0]
237 elif isinstance(l[0], list):
238 return control(None, l[0], cn)
239 else:
240 return False
241 elif len(l) > 2:
242 return control_list(ris, l, cn)
243 else:
244 return False
245
246 def control_list(ris, l, cn):
247 while l:
248 if len(l) == 1:
249 return False
250 elif (isinstance(l[0], (float, int)) or
251 l[0] == None) and l[1] in ['and', 'or']:
252 if isinstance(l[2], (float, int)) or l[2] == None:
253 ris = computes(l[0], l[1], l[2], cn)
254 elif isinstance(l[2], list):
255 tmp = control(None, l[2], cn)
256 if tmp is False:
257 return False
258 else:
259 ris = computes(l[0], l[1], tmp, cn)
260 else:
261 return False
262 l = l[3:]
263 elif l[0] in ['and', 'or']:
264 if isinstance(l[1], (float, int)) or l[1] == None:
265 ris = computes(ris, l[0], l[1], cn)
266 elif isinstance(l[1], list):
267 tmp = control(None,l[1], cn)
268 if tmp is False:
269 return False
270 else:
271 ris = computes(ris, l[0], tmp, cn)
272 else:
273 return False
274 l = l[2:]
275 elif isinstance(l[0], list) and l[1] in ['and', 'or']:
276 if isinstance(l[2], (float, int)) or l[2] == None:
277 tmp = control(None, l[0], cn)
278 if tmp is False:
279 return False
280 else:
281 ris = computes(tmp, l[1], l[2], cn)
282 elif isinstance(l[2], list):
283 tmp = control(None, l[0], cn)
284 tmp2 = control(None, l[2], cn)
285 if tmp is False or tmp2 is False:
286 return False
287 else:
288 ris = computes(tmp, l[1], tmp2, cn)
289 else:
290 return False
291 l = l[3:]
292 else:
293 return False
294 return ris
295
296 ############################ map_methods ######################################
297
298 def fold_change(avg1, avg2):
299 if avg1 == 0 and avg2 == 0:
300 return 0
301 elif avg1 == 0:
302 return '-INF'
303 elif avg2 == 0:
304 return 'INF'
305 else:
306 return math.log(avg1 / avg2, 2)
307
308 def fix_style(l, col, width, dash):
309 tmp = l.split(';')
310 flag_col = False
311 flag_width = False
312 flag_dash = False
313 for i in range(len(tmp)):
314 if tmp[i].startswith('stroke:'):
315 tmp[i] = 'stroke:' + col
316 flag_col = True
317 if tmp[i].startswith('stroke-width:'):
318 tmp[i] = 'stroke-width:' + width
319 flag_width = True
320 if tmp[i].startswith('stroke-dasharray:'):
321 tmp[i] = 'stroke-dasharray:' + dash
322 flag_dash = True
323 if not flag_col:
324 tmp.append('stroke:' + col)
325 if not flag_width:
326 tmp.append('stroke-width:' + width)
327 if not flag_dash:
328 tmp.append('stroke-dasharray:' + dash)
329 return ';'.join(tmp)
330
331 def fix_map(d, core_map, threshold_P_V, threshold_F_C, max_F_C):
332 maxT = 12
333 minT = 2
334 grey = '#BEBEBE'
335 blue = '#0000FF'
336 red = '#E41A1C'
337 for el in core_map.iter():
338 el_id = str(el.get('id'))
339 if el_id.startswith('R_'):
340 tmp = d.get(el_id[2:])
341 if tmp != None:
342 p_val = tmp[0]
343 f_c = tmp[1]
344 if p_val < threshold_P_V:
345 if not isinstance(f_c, str):
346 if abs(f_c) < math.log(threshold_F_C, 2):
347 col = grey
348 width = str(minT)
349 else:
350 if f_c < 0:
351 col = blue
352 elif f_c > 0:
353 col = red
354 width = str(max((abs(f_c) * maxT) / max_F_C, minT))
355 else:
356 if f_c == '-INF':
357 col = blue
358 elif f_c == 'INF':
359 col = red
360 width = str(maxT)
361 dash = 'none'
362 else:
363 dash = '5,5'
364 col = grey
365 width = str(minT)
366 el.set('style', fix_style(el.get('style'), col, width, dash))
367 return core_map
368
369 ############################ make recon #######################################
370
371 def check_and_doWord(l):
372 tmp = []
373 tmp_genes = []
374 count = 0
375 while l:
376 if count >= 0:
377 if l[0] == '(':
378 count += 1
379 tmp.append(l[0])
380 l.pop(0)
381 elif l[0] == ')':
382 count -= 1
383 tmp.append(l[0])
384 l.pop(0)
385 elif l[0] == ' ':
386 l.pop(0)
387 else:
388 word = []
389 while l:
390 if l[0] in [' ', '(', ')']:
391 break
392 else:
393 word.append(l[0])
394 l.pop(0)
395 word = ''.join(word)
396 tmp.append(word)
397 if not(word in ['or', 'and']):
398 tmp_genes.append(word)
399 else:
400 return False
401 if count == 0:
402 return (tmp, tmp_genes)
403 else:
404 return False
405
406 def brackets_to_list(l):
407 tmp = []
408 while l:
409 if l[0] == '(':
410 l.pop(0)
411 tmp.append(resolve_brackets(l))
412 else:
413 tmp.append(l[0])
414 l.pop(0)
415 return tmp
416
417 def resolve_brackets(l):
418 tmp = []
419 while l[0] != ')':
420 if l[0] == '(':
421 l.pop(0)
422 tmp.append(resolve_brackets(l))
423 else:
424 tmp.append(l[0])
425 l.pop(0)
426 l.pop(0)
427 return tmp
428
429 def priorityAND(l):
430 tmp = []
431 flag = True
432 while l:
433 if len(l) == 1:
434 if isinstance(l[0], list):
435 tmp.append(priorityAND(l[0]))
436 else:
437 tmp.append(l[0])
438 l = l[1:]
439 elif l[0] == 'or':
440 tmp.append(l[0])
441 flag = False
442 l = l[1:]
443 elif l[1] == 'or':
444 if isinstance(l[0], list):
445 tmp.append(priorityAND(l[0]))
446 else:
447 tmp.append(l[0])
448 tmp.append(l[1])
449 flag = False
450 l = l[2:]
451 elif l[1] == 'and':
452 tmpAnd = []
453 if isinstance(l[0], list):
454 tmpAnd.append(priorityAND(l[0]))
455 else:
456 tmpAnd.append(l[0])
457 tmpAnd.append(l[1])
458 if isinstance(l[2], list):
459 tmpAnd.append(priorityAND(l[2]))
460 else:
461 tmpAnd.append(l[2])
462 l = l[3:]
463 while l:
464 if l[0] == 'and':
465 tmpAnd.append(l[0])
466 if isinstance(l[1], list):
467 tmpAnd.append(priorityAND(l[1]))
468 else:
469 tmpAnd.append(l[1])
470 l = l[2:]
471 elif l[0] == 'or':
472 flag = False
473 break
474 if flag == True: #when there are only AND in list
475 tmp.extend(tmpAnd)
476 elif flag == False:
477 tmp.append(tmpAnd)
478 return tmp
479
480 def checkRule(l):
481 if len(l) == 1:
482 if isinstance(l[0], list):
483 if checkRule(l[0]) is False:
484 return False
485 elif len(l) > 2:
486 if checkRule2(l) is False:
487 return False
488 else:
489 return False
490 return True
491
492 def checkRule2(l):
493 while l:
494 if len(l) == 1:
495 return False
496 elif isinstance(l[0], list) and l[1] in ['and', 'or']:
497 if checkRule(l[0]) is False:
498 return False
499 if isinstance(l[2], list):
500 if checkRule(l[2]) is False:
501 return False
502 l = l[3:]
503 elif l[1] in ['and', 'or']:
504 if isinstance(l[2], list):
505 if checkRule(l[2]) is False:
506 return False
507 l = l[3:]
508 elif l[0] in ['and', 'or']:
509 if isinstance(l[1], list):
510 if checkRule(l[1]) is False:
511 return False
512 l = l[2:]
513 else:
514 return False
515 return True
516
517 def do_rules(rules):
518 split_rules = []
519 err_rules = []
520 tmp_gene_in_rule = []
521 for i in range(len(rules)):
522 tmp = list(rules[i])
523 if tmp:
524 tmp, tmp_genes = check_and_doWord(tmp)
525 tmp_gene_in_rule.extend(tmp_genes)
526 if tmp is False:
527 split_rules.append([])
528 err_rules.append(rules[i])
529 else:
530 tmp = brackets_to_list(tmp)
531 if checkRule(tmp):
532 split_rules.append(priorityAND(tmp))
533 else:
534 split_rules.append([])
535 err_rules.append(rules[i])
536 else:
537 split_rules.append([])
538 if err_rules:
539 warning('Warning: wrong format rule in ' + str(err_rules) + '\n')
540 return (split_rules, list(set(tmp_gene_in_rule)))
541
542 def make_recon(data):
543 try:
544 import cobra as cb
545 import warnings
546 with warnings.catch_warnings():
547 warnings.simplefilter('ignore')
548 recon = cb.io.read_sbml_model(data)
549 react = recon.reactions
550 rules = [react[i].gene_reaction_rule for i in range(len(react))]
551 ids = [react[i].id for i in range(len(react))]
552 except cb.io.sbml3.CobraSBMLError:
553 try:
554 data = (pd.read_csv(data, sep = '\t', dtype = str, engine='python')).fillna('')
555 if len(data.columns) < 2:
556 sys.exit('Execution aborted: wrong format of '+
557 'custom datarules\n')
558 if not len(data.columns) == 2:
559 warning('Warning: more than 2 columns in custom datarules.\n' +
560 'Extra columns have been disregarded\n')
561 ids = list(data.iloc[:, 0])
562 rules = list(data.iloc[:, 1])
563 except pd.errors.EmptyDataError:
564 sys.exit('Execution aborted: wrong format of custom datarules\n')
565 except pd.errors.ParserError:
566 sys.exit('Execution aborted: wrong format of custom datarules\n')
567 split_rules, tmp_genes = do_rules(rules)
568 gene_in_rule = {}
569 for i in tmp_genes:
570 gene_in_rule[i] = 'ok'
571 return (ids, split_rules, gene_in_rule)
572
573 ############################ gene #############################################
574
575 def data_gene(gene, type_gene, name, gene_custom):
576 args = process_args(sys.argv)
577 for i in range(len(gene)):
578 tmp = gene.iloc[i, 0]
579 if tmp.startswith(' ') or tmp.endswith(' '):
580 gene.iloc[i, 0] = (tmp.lstrip()).rstrip()
581 gene_dup = [item for item, count in
582 collections.Counter(gene[gene.columns[0]]).items() if count > 1]
583 pat_dup = [item for item, count in
584 collections.Counter(list(gene.columns)).items() if count > 1]
585 if gene_dup:
586 if gene_custom == None:
587 if args.rules_selector == 'HMRcore':
588 gene_in_rule = pk.load(open(args.tool_dir +
589 '/local/HMRcore_genes.p', 'rb'))
590 elif args.rules_selector == 'Recon':
591 gene_in_rule = pk.load(open(args.tool_dir +
592 '/local/Recon_genes.p', 'rb'))
593 gene_in_rule = gene_in_rule.get(type_gene)
594 else:
595 gene_in_rule = gene_custom
596 tmp = []
597 for i in gene_dup:
598 if gene_in_rule.get(i) == 'ok':
599 tmp.append(i)
600 if tmp:
601 sys.exit('Execution aborted because gene ID '
602 +str(tmp)+' in '+name+' is duplicated\n')
603 if pat_dup:
604 warning('Warning: duplicated label\n' + str(pat_dup) + 'in ' + name +
605 '\n')
606 return (gene.set_index(gene.columns[0])).to_dict()
607
608 ############################ resolve ##########################################
609
610 def resolve(genes, rules, ids, resolve_none, name):
611 resolve_rules = {}
612 not_found = []
613 flag = False
614 for key, value in genes.items():
615 tmp_resolve = []
616 for i in range(len(rules)):
617 tmp = rules[i]
618 if tmp:
619 tmp, err = replace_gene_value(tmp, value)
620 if err:
621 not_found.extend(err)
622 ris = control(None, tmp, resolve_none)
623 if ris is False or ris == None:
624 tmp_resolve.append(None)
625 else:
626 tmp_resolve.append(ris)
627 flag = True
628 else:
629 tmp_resolve.append(None)
630 resolve_rules[key] = tmp_resolve
631 if flag is False:
632 warning('Warning: no computable score (due to missing gene values)' +
633 'for class ' + name + ', the class has been disregarded\n')
634 return (None, None)
635 return (resolve_rules, list(set(not_found)))
636
637 ############################ split class ######################################
638
639 def split_class(classes, resolve_rules):
640 class_pat = {}
641 for i in range(len(classes)):
642 classe = classes.iloc[i, 1]
643 if not pd.isnull(classe):
644 l = []
645 for j in range(i, len(classes)):
646 if classes.iloc[j, 1] == classe:
647 pat_id = classes.iloc[j, 0]
648 tmp = resolve_rules.get(pat_id, None)
649 if tmp != None:
650 l.append(tmp)
651 classes.iloc[j, 1] = None
652 if l:
653 class_pat[classe] = list(map(list, zip(*l)))
654 else:
655 warning('Warning: no sample found in class ' + classe +
656 ', the class has been disregarded\n')
657 return class_pat
658
659 ############################ create_ras #######################################
660
661 def create_ras (resolve_rules, dataset_name):
662
663 if resolve_rules == None:
664 warning("Couldn't generate RAS for current dataset: " + dataset_name)
665
666 for geni in resolve_rules.values():
667 for i, valori in enumerate(geni):
668 if valori == None:
669 geni[i] = 'None'
670
671 output_ras = pd.DataFrame.from_dict(resolve_rules)
672 output_to_csv = pd.DataFrame.to_csv(output_ras, sep = '\t', index = False)
673
674 text_file = open("ras/Reaction_Activity_Score_Of_" + dataset_name + ".tsv", "w")
675 text_file.write(output_to_csv)
676 text_file.close()
677
678 ############################ map ##############################################
679
680 def maps(core_map, class_pat, ids, threshold_P_V, threshold_F_C, create_svg, create_pdf):
681 args = process_args(sys.argv)
682 if (not class_pat) or (len(class_pat.keys()) < 2):
683 sys.exit('Execution aborted: classes provided for comparisons are ' +
684 'less than two\n')
685 for i, j in it.combinations(class_pat.keys(), 2):
686 tmp = {}
687 count = 0
688 max_F_C = 0
689 for l1, l2 in zip(class_pat.get(i), class_pat.get(j)):
690 try:
691 stat_D, p_value = st.ks_2samp(l1, l2)
692 avg = fold_change(sum(l1) / len(l1), sum(l2) / len(l2))
693 if not isinstance(avg, str):
694 if max_F_C < abs(avg):
695 max_F_C = abs(avg)
696 tmp[ids[count]] = [float(p_value), avg]
697 count += 1
698 except (TypeError, ZeroDivisionError):
699 count += 1
700 tab = 'result/' + i + '_vs_' + j + ' (Tabular Result).tsv'
701 tmp_csv = pd.DataFrame.from_dict(tmp, orient = "index")
702 tmp_csv = tmp_csv.reset_index()
703 header = ['ids', 'P_Value', 'Average']
704 tmp_csv.to_csv(tab, sep = '\t', index = False, header = header)
705
706 if create_svg or create_pdf:
707 if args.rules_selector == 'HMRcore' or (args.rules_selector == 'Custom'
708 and args.yes_no == 'yes'):
709 fix_map(tmp, core_map, threshold_P_V, threshold_F_C, max_F_C)
710 file_svg = 'result/' + i + '_vs_' + j + ' (SVG Map).svg'
711 with open(file_svg, 'wb') as new_map:
712 new_map.write(ET.tostring(core_map))
713
714
715 if create_pdf:
716 file_pdf = 'result/' + i + '_vs_' + j + ' (PDF Map).pdf'
717 renderPDF.drawToFile(svg2rlg(file_svg), file_pdf)
718
719 if not create_svg:
720 #Ho utilizzato il file svg per generare il pdf,
721 #ma l'utente non ne ha richiesto il ritorno, quindi
722 #lo elimino
723 os.remove('result/' + i + '_vs_' + j + ' (SVG Map).svg')
724
725 return None
726
727 ############################ MAIN #############################################
728
729 def main():
730 args = process_args(sys.argv)
731
732 create_svg = check_bool(args.generate_svg)
733 create_pdf = check_bool(args.generate_pdf)
734 generate_ras = check_bool(args.generate_ras)
735
736 os.makedirs('result')
737
738 if generate_ras:
739 os.makedirs('ras')
740
741 if args.rules_selector == 'HMRcore':
742 recon = pk.load(open(args.tool_dir + '/local/HMRcore_rules.p', 'rb'))
743 elif args.rules_selector == 'Recon':
744 recon = pk.load(open(args.tool_dir + '/local/Recon_rules.p', 'rb'))
745 elif args.rules_selector == 'Custom':
746 ids, rules, gene_in_rule = make_recon(args.custom)
747
748 resolve_none = check_bool(args.none)
749
750 class_pat = {}
751
752 if args.option == 'datasets':
753 num = 1
754 for i, j in zip(args.input_datas, args.names):
755
756 name = name_dataset(j, num)
757 dataset = read_dataset(i, name)
758
759 dataset.iloc[:, 0] = (dataset.iloc[:, 0]).astype(str)
760
761 type_gene = gene_type(dataset.iloc[0, 0], name)
762
763 if args.rules_selector != 'Custom':
764 genes = data_gene(dataset, type_gene, name, None)
765 ids, rules = load_id_rules(recon.get(type_gene))
766 elif args.rules_selector == 'Custom':
767 genes = data_gene(dataset, type_gene, name, gene_in_rule)
768
769 resolve_rules, err = resolve(genes, rules, ids, resolve_none, name)
770
771 if generate_ras:
772 create_ras(resolve_rules, name)
773
774
775 if err != None and err:
776 warning('Warning: gene\n' + str(err) + '\nnot found in class '
777 + name + ', the expression level for this gene ' +
778 'will be considered NaN\n')
779 if resolve_rules != None:
780 class_pat[name] = list(map(list, zip(*resolve_rules.values())))
781 num += 1
782 elif args.option == 'dataset_class':
783 name = 'RNAseq'
784 dataset = read_dataset(args.input_data, name)
785 dataset.iloc[:, 0] = (dataset.iloc[:, 0]).astype(str)
786 type_gene = gene_type(dataset.iloc[0, 0], name)
787 classes = read_dataset(args.input_class, 'class')
788 if not len(classes.columns) == 2:
789 warning('Warning: more than 2 columns in class file. Extra' +
790 'columns have been disregarded\n')
791 classes = classes.astype(str)
792 if args.rules_selector != 'Custom':
793 genes = data_gene(dataset, type_gene, name, None)
794 ids, rules = load_id_rules(recon.get(type_gene))
795 elif args.rules_selector == 'Custom':
796 genes = data_gene(dataset, type_gene, name, gene_in_rule)
797 resolve_rules, err = resolve(genes, rules, ids, resolve_none, name)
798 if err != None and err:
799 warning('Warning: gene\n'+str(err)+'\nnot found in class '
800 + name + ', the expression level for this gene ' +
801 'will be considered NaN\n')
802 if resolve_rules != None:
803 class_pat = split_class(classes, resolve_rules)
804
805 if args.rules_selector == 'Custom':
806 if args.yes_no == 'yes':
807 try:
808 core_map = ET.parse(args.custom_map)
809 except (ET.XMLSyntaxError, ET.XMLSchemaParseError):
810 sys.exit('Execution aborted: custom map in wrong format')
811 elif args.yes_no == 'no':
812 core_map = ET.parse(args.tool_dir + '/local/HMRcoreMap.svg')
813 else:
814 core_map = ET.parse(args.tool_dir+'/local/HMRcoreMap.svg')
815
816 maps(core_map, class_pat, ids, args.pValue, args.fChange, create_svg, create_pdf)
817
818 print('Execution succeded')
819
820 return None
821
822 ###############################################################################
823
824 if __name__ == "__main__":
825 main()