|
0
|
1 from __future__ import division
|
|
|
2 import os
|
|
|
3 import sys
|
|
|
4 import pandas as pd
|
|
|
5 import collections
|
|
|
6 import pickle as pk
|
|
|
7 import argparse
|
|
|
8 from sklearn.cluster import KMeans
|
|
13
|
9 import matplotlib
|
|
14
|
10 # Force matplotlib to not use any Xwindows backend.
|
|
|
11 matplotlib.use('Agg')
|
|
0
|
12 import matplotlib.pyplot as plt
|
|
|
13
|
|
|
14 ########################## argparse ###########################################
|
|
|
15
|
|
|
16 def process_args(args):
|
|
|
17 parser = argparse.ArgumentParser(usage = '%(prog)s [options]',
|
|
|
18 description = 'process some value\'s' +
|
|
|
19 ' genes to create class.')
|
|
|
20 parser.add_argument('-rs', '--rules_selector',
|
|
|
21 type = str,
|
|
|
22 default = 'HMRcore',
|
|
|
23 choices = ['HMRcore', 'Recon', 'Custom'],
|
|
|
24 help = 'chose which type of dataset you want use')
|
|
|
25 parser.add_argument('-cr', '--custom',
|
|
|
26 type = str,
|
|
|
27 help='your dataset if you want custom rules')
|
|
|
28 parser.add_argument('-ch', '--cond_hier',
|
|
|
29 type = str,
|
|
|
30 default = 'no',
|
|
|
31 choices = ['no', 'yes'],
|
|
|
32 help = 'chose if you wanna hierical dendrogram')
|
|
|
33 parser.add_argument('-lk', '--k_min',
|
|
|
34 type = int,
|
|
|
35 help = 'min number of cluster')
|
|
|
36 parser.add_argument('-uk', '--k_max',
|
|
|
37 type = int,
|
|
|
38 help = 'max number of cluster')
|
|
|
39 parser.add_argument('-li', '--linkage',
|
|
|
40 type = str,
|
|
|
41 choices = ['single', 'complete', 'average'],
|
|
|
42 help='linkage hierarchical cluster')
|
|
|
43 parser.add_argument('-d', '--data',
|
|
|
44 type = str,
|
|
|
45 required = True,
|
|
|
46 help = 'input dataset')
|
|
|
47 parser.add_argument('-n', '--none',
|
|
|
48 type = str,
|
|
|
49 default = 'true',
|
|
|
50 choices = ['true', 'false'],
|
|
|
51 help = 'compute Nan values')
|
|
|
52 parser.add_argument('-td', '--tool_dir',
|
|
|
53 type = str,
|
|
|
54 required = True,
|
|
|
55 help = 'your tool directory')
|
|
|
56 parser.add_argument('-na', '--name',
|
|
|
57 type = str,
|
|
|
58 help = 'name of dataset')
|
|
|
59 parser.add_argument('-de', '--dendro',
|
|
|
60 help = "Dendrogram out")
|
|
|
61 parser.add_argument('-ol', '--out_log',
|
|
|
62 help = "Output log")
|
|
|
63 parser.add_argument('-el', '--elbow',
|
|
|
64 help = "Out elbow")
|
|
|
65 args = parser.parse_args()
|
|
|
66 return args
|
|
|
67
|
|
|
68 ########################### warning ###########################################
|
|
|
69
|
|
|
70 def warning(s):
|
|
|
71 args = process_args(sys.argv)
|
|
|
72 with open(args.out_log, 'a') as log:
|
|
|
73 log.write(s)
|
|
|
74
|
|
|
75 ############################ dataset input ####################################
|
|
|
76
|
|
|
77 def read_dataset(data, name):
|
|
|
78 try:
|
|
|
79 dataset = pd.read_csv(data, sep = '\t', header = 0)
|
|
|
80 except pd.errors.EmptyDataError:
|
|
|
81 sys.exit('Execution aborted: wrong format of '+name+'\n')
|
|
|
82 if len(dataset.columns) < 2:
|
|
|
83 sys.exit('Execution aborted: wrong format of '+name+'\n')
|
|
|
84 return dataset
|
|
|
85
|
|
|
86 ############################ dataset name #####################################
|
|
|
87
|
|
|
88 def name_dataset(name_data, count):
|
|
|
89 if str(name_data) == 'Dataset':
|
|
|
90 return str(name_data) + '_' + str(count)
|
|
|
91 else:
|
|
|
92 return str(name_data)
|
|
|
93
|
|
|
94 ############################ load id e rules ##################################
|
|
|
95
|
|
|
96 def load_id_rules(reactions):
|
|
|
97 ids, rules = [], []
|
|
|
98 for key, value in reactions.items():
|
|
|
99 ids.append(key)
|
|
|
100 rules.append(value)
|
|
|
101 return (ids, rules)
|
|
|
102
|
|
|
103 ############################ check_methods ####################################
|
|
|
104
|
|
|
105 def gene_type(l, name):
|
|
|
106 if check_hgnc(l):
|
|
|
107 return 'hugo_id'
|
|
|
108 elif check_ensembl(l):
|
|
|
109 return 'ensembl_gene_id'
|
|
|
110 elif check_symbol(l):
|
|
|
111 return 'symbol'
|
|
|
112 elif check_entrez(l):
|
|
|
113 return 'entrez_id'
|
|
|
114 else:
|
|
|
115 sys.exit('Execution aborted:\n' +
|
|
|
116 'gene ID type in ' + name + ' not supported. Supported ID' +
|
|
|
117 'types are: HUGO ID, Ensemble ID, HUGO symbol, Entrez ID\n')
|
|
|
118
|
|
|
119 def check_hgnc(l):
|
|
|
120 if len(l) > 5:
|
|
|
121 if (l.upper()).startswith('HGNC:'):
|
|
|
122 return l[5:].isdigit()
|
|
|
123 else:
|
|
|
124 return False
|
|
|
125 else:
|
|
|
126 return False
|
|
|
127
|
|
|
128 def check_ensembl(l):
|
|
|
129 if len(l) == 15:
|
|
|
130 if (l.upper()).startswith('ENS'):
|
|
|
131 return l[4:].isdigit()
|
|
|
132 else:
|
|
|
133 return False
|
|
|
134 else:
|
|
|
135 return False
|
|
|
136
|
|
|
137 def check_symbol(l):
|
|
|
138 if len(l) > 0:
|
|
|
139 if l[0].isalpha() and l[1:].isalnum():
|
|
|
140 return True
|
|
|
141 else:
|
|
|
142 return False
|
|
|
143 else:
|
|
|
144 return False
|
|
|
145
|
|
|
146 def check_entrez(l):
|
|
|
147 if len(l) > 0:
|
|
|
148 return l.isdigit()
|
|
|
149 else:
|
|
|
150 return False
|
|
|
151
|
|
|
152 def check_bool(b):
|
|
|
153 if b == 'true':
|
|
|
154 return True
|
|
|
155 elif b == 'false':
|
|
|
156 return False
|
|
|
157
|
|
|
158 ############################ make recon #######################################
|
|
|
159
|
|
|
160 def check_and_doWord(l):
|
|
|
161 tmp = []
|
|
|
162 tmp_genes = []
|
|
|
163 count = 0
|
|
|
164 while l:
|
|
|
165 if count >= 0:
|
|
|
166 if l[0] == '(':
|
|
|
167 count += 1
|
|
|
168 tmp.append(l[0])
|
|
|
169 l.pop(0)
|
|
|
170 elif l[0] == ')':
|
|
|
171 count -= 1
|
|
|
172 tmp.append(l[0])
|
|
|
173 l.pop(0)
|
|
|
174 elif l[0] == ' ':
|
|
|
175 l.pop(0)
|
|
|
176 else:
|
|
|
177 word = []
|
|
|
178 while l:
|
|
|
179 if l[0] in [' ', '(', ')']:
|
|
|
180 break
|
|
|
181 else:
|
|
|
182 word.append(l[0])
|
|
|
183 l.pop(0)
|
|
|
184 word = ''.join(word)
|
|
|
185 tmp.append(word)
|
|
|
186 if not(word in ['or', 'and']):
|
|
|
187 tmp_genes.append(word)
|
|
|
188 else:
|
|
|
189 return False
|
|
|
190 if count == 0:
|
|
|
191 return (tmp, tmp_genes)
|
|
|
192 else:
|
|
|
193 return False
|
|
|
194
|
|
|
195 def brackets_to_list(l):
|
|
|
196 tmp = []
|
|
|
197 while l:
|
|
|
198 if l[0] == '(':
|
|
|
199 l.pop(0)
|
|
|
200 tmp.append(resolve_brackets(l))
|
|
|
201 else:
|
|
|
202 tmp.append(l[0])
|
|
|
203 l.pop(0)
|
|
|
204 return tmp
|
|
|
205
|
|
|
206 def resolve_brackets(l):
|
|
|
207 tmp = []
|
|
|
208 while l[0] != ')':
|
|
|
209 if l[0] == '(':
|
|
|
210 l.pop(0)
|
|
|
211 tmp.append(resolve_brackets(l))
|
|
|
212 else:
|
|
|
213 tmp.append(l[0])
|
|
|
214 l.pop(0)
|
|
|
215 l.pop(0)
|
|
|
216 return tmp
|
|
|
217
|
|
|
218 def priorityAND(l):
|
|
|
219 tmp = []
|
|
|
220 flag = True
|
|
|
221 while l:
|
|
|
222 if len(l) == 1:
|
|
|
223 if isinstance(l[0], list):
|
|
|
224 tmp.append(priorityAND(l[0]))
|
|
|
225 else:
|
|
|
226 tmp.append(l[0])
|
|
|
227 l = l[1:]
|
|
|
228 elif l[0] == 'or':
|
|
|
229 tmp.append(l[0])
|
|
|
230 flag = False
|
|
|
231 l = l[1:]
|
|
|
232 elif l[1] == 'or':
|
|
|
233 if isinstance(l[0], list):
|
|
|
234 tmp.append(priorityAND(l[0]))
|
|
|
235 else:
|
|
|
236 tmp.append(l[0])
|
|
|
237 tmp.append(l[1])
|
|
|
238 flag = False
|
|
|
239 l = l[2:]
|
|
|
240 elif l[1] == 'and':
|
|
|
241 tmpAnd = []
|
|
|
242 if isinstance(l[0], list):
|
|
|
243 tmpAnd.append(priorityAND(l[0]))
|
|
|
244 else:
|
|
|
245 tmpAnd.append(l[0])
|
|
|
246 tmpAnd.append(l[1])
|
|
|
247 if isinstance(l[2], list):
|
|
|
248 tmpAnd.append(priorityAND(l[2]))
|
|
|
249 else:
|
|
|
250 tmpAnd.append(l[2])
|
|
|
251 l = l[3:]
|
|
|
252 while l:
|
|
|
253 if l[0] == 'and':
|
|
|
254 tmpAnd.append(l[0])
|
|
|
255 if isinstance(l[1], list):
|
|
|
256 tmpAnd.append(priorityAND(l[1]))
|
|
|
257 else:
|
|
|
258 tmpAnd.append(l[1])
|
|
|
259 l = l[2:]
|
|
|
260 elif l[0] == 'or':
|
|
|
261 flag = False
|
|
|
262 break
|
|
|
263 if flag == True: #se ci sono solo AND nella lista
|
|
|
264 tmp.extend(tmpAnd)
|
|
|
265 elif flag == False:
|
|
|
266 tmp.append(tmpAnd)
|
|
|
267 return tmp
|
|
|
268
|
|
|
269 def checkRule(l):
|
|
|
270 if len(l) == 1:
|
|
|
271 if isinstance(l[0], list):
|
|
|
272 if checkRule(l[0]) is False:
|
|
|
273 return False
|
|
|
274 elif len(l) > 2:
|
|
|
275 if checkRule2(l) is False:
|
|
|
276 return False
|
|
|
277 else:
|
|
|
278 return False
|
|
|
279 return True
|
|
|
280
|
|
|
281 def checkRule2(l):
|
|
|
282 while l:
|
|
|
283 if len(l) == 1:
|
|
|
284 return False
|
|
|
285 elif isinstance(l[0], list) and l[1] in ['and', 'or']:
|
|
|
286 if checkRule(l[0]) is False:
|
|
|
287 return False
|
|
|
288 if isinstance(l[2], list):
|
|
|
289 if checkRule(l[2]) is False:
|
|
|
290 return False
|
|
|
291 l = l[3:]
|
|
|
292 elif l[1] in ['and', 'or']:
|
|
|
293 if isinstance(l[2], list):
|
|
|
294 if checkRule(l[2]) is False:
|
|
|
295 return False
|
|
|
296 l = l[3:]
|
|
|
297 elif l[0] in ['and', 'or']:
|
|
|
298 if isinstance(l[1], list):
|
|
|
299 if checkRule(l[1]) is False:
|
|
|
300 return False
|
|
|
301 l = l[2:]
|
|
|
302 else:
|
|
|
303 return False
|
|
|
304 return True
|
|
|
305
|
|
|
306 def do_rules(rules):
|
|
|
307 split_rules = []
|
|
|
308 err_rules = []
|
|
|
309 tmp_gene_in_rule = []
|
|
|
310 for i in range(len(rules)):
|
|
|
311 tmp = list(rules[i])
|
|
|
312 if tmp:
|
|
|
313 tmp, tmp_genes = check_and_doWord(tmp)
|
|
|
314 tmp_gene_in_rule.extend(tmp_genes)
|
|
|
315 if tmp is False:
|
|
|
316 split_rules.append([])
|
|
|
317 err_rules.append(rules[i])
|
|
|
318 else:
|
|
|
319 tmp = brackets_to_list(tmp)
|
|
|
320 if checkRule(tmp):
|
|
|
321 split_rules.append(priorityAND(tmp))
|
|
|
322 else:
|
|
|
323 split_rules.append([])
|
|
|
324 err_rules.append(rules[i])
|
|
|
325 else:
|
|
|
326 split_rules.append([])
|
|
|
327 if err_rules:
|
|
|
328 warning('Warning: wrong format rule in ' + str(err_rules) + '\n')
|
|
|
329 return (split_rules, list(set(tmp_gene_in_rule)))
|
|
|
330
|
|
|
331 def make_recon(data):
|
|
|
332 try:
|
|
|
333 import cobra as cb
|
|
|
334 import warnings
|
|
|
335 with warnings.catch_warnings():
|
|
|
336 warnings.simplefilter('ignore')
|
|
|
337 recon = cb.io.read_sbml_model(data)
|
|
|
338 react = recon.reactions
|
|
|
339 rules = [react[i].gene_reaction_rule for i in range(len(react))]
|
|
|
340 ids = [react[i].id for i in range(len(react))]
|
|
|
341 except cb.io.sbml3.CobraSBMLError:
|
|
|
342 try:
|
|
|
343 data = (pd.read_csv(data, sep = '\t', dtype = str)).fillna('')
|
|
|
344 if len(data.columns) < 2:
|
|
|
345 sys.exit('Execution aborted: wrong format of ' +
|
|
|
346 'custom GPR rules\n')
|
|
|
347 if not len(data.columns) == 2:
|
|
|
348 warning('WARNING: more than 2 columns in custom GPR rules.\n' +
|
|
|
349 'Extra columns have been disregarded\n')
|
|
|
350 ids = list(data.iloc[:, 0])
|
|
|
351 rules = list(data.iloc[:, 1])
|
|
|
352 except pd.errors.EmptyDataError:
|
|
|
353 sys.exit('Execution aborted: wrong format of custom GPR rules\n')
|
|
|
354 except pd.errors.ParserError:
|
|
|
355 sys.exit('Execution aborted: wrong format of custom GPR rules\n')
|
|
|
356 split_rules, tmp_genes = do_rules(rules)
|
|
|
357 gene_in_rule = {}
|
|
|
358 for i in tmp_genes:
|
|
|
359 gene_in_rule[i] = 'ok'
|
|
|
360 return (ids, split_rules, gene_in_rule)
|
|
|
361
|
|
|
362 ############################ resolve_methods ##################################
|
|
|
363
|
|
|
364 def replace_gene_value(l, d):
|
|
|
365 tmp = []
|
|
|
366 err = []
|
|
|
367 while l:
|
|
|
368 if isinstance(l[0], list):
|
|
|
369 tmp_rules, tmp_err = replace_gene_value(l[0], d)
|
|
|
370 tmp.append(tmp_rules)
|
|
|
371 err.extend(tmp_err)
|
|
|
372 else:
|
|
|
373 value = replace_gene(l[0],d)
|
|
|
374 tmp.append(value)
|
|
|
375 if value == None:
|
|
|
376 err.append(l[0])
|
|
|
377 l = l[1:]
|
|
|
378 return (tmp, err)
|
|
|
379
|
|
|
380 def replace_gene(l, d):
|
|
|
381 if l =='and' or l == 'or':
|
|
|
382 return l
|
|
|
383 else:
|
|
|
384 value = d.get(l, None)
|
|
|
385 if not(value == None or isinstance(value, (int, float))):
|
|
|
386 sys.exit('Execution aborted: ' + value + ' value not valid\n')
|
|
|
387 return value
|
|
|
388
|
|
|
389 def compute(val1, op, val2, cn):
|
|
|
390 if val1 != None and val2 != None:
|
|
|
391 if op == 'and':
|
|
|
392 return min(val1, val2)
|
|
|
393 else:
|
|
|
394 return val1 + val2
|
|
|
395 elif op == 'and':
|
|
|
396 if cn is True:
|
|
|
397 if val1 != None:
|
|
|
398 return val1
|
|
|
399 elif val2 != None:
|
|
|
400 return val2
|
|
|
401 else:
|
|
|
402 return None
|
|
|
403 else:
|
|
|
404 return None
|
|
|
405 else:
|
|
|
406 if val1 != None:
|
|
|
407 return val1
|
|
|
408 elif val2 != None:
|
|
|
409 return val2
|
|
|
410 else:
|
|
|
411 return None
|
|
|
412
|
|
|
413 def control(ris, l, cn):
|
|
|
414 if len(l) == 1:
|
|
|
415 if isinstance(l[0], (float, int)) or l[0] == None:
|
|
|
416 return l[0]
|
|
|
417 elif isinstance(l[0], list):
|
|
|
418 return control(None, l[0], cn)
|
|
|
419 else:
|
|
|
420 return False
|
|
|
421 elif len(l) > 2:
|
|
|
422 return control_list(ris, l, cn)
|
|
|
423 else:
|
|
|
424 return False
|
|
|
425
|
|
|
426 def control_list(ris, l, cn):
|
|
|
427 while l:
|
|
|
428 if len(l) == 1:
|
|
|
429 return False
|
|
|
430 elif (isinstance(l[0], (float, int)) or
|
|
|
431 l[0] == None) and l[1] in ['and', 'or']:
|
|
|
432 if isinstance(l[2], (float, int)) or l[2] == None:
|
|
|
433 ris = compute(l[0], l[1], l[2], cn)
|
|
|
434 elif isinstance(l[2], list):
|
|
|
435 tmp = control(None, l[2], cn)
|
|
|
436 if tmp is False:
|
|
|
437 return False
|
|
|
438 else:
|
|
|
439 ris = compute(l[0], l[1], tmp, cn)
|
|
|
440 else:
|
|
|
441 return False
|
|
|
442 l = l[3:]
|
|
|
443 elif l[0] in ['and', 'or']:
|
|
|
444 if isinstance(l[1], (float, int)) or l[1] == None:
|
|
|
445 ris = compute(ris, l[0], l[1], cn)
|
|
|
446 elif isinstance(l[1], list):
|
|
|
447 tmp = control(None,l[1], cn)
|
|
|
448 if tmp is False:
|
|
|
449 return False
|
|
|
450 else:
|
|
|
451 ris = compute(ris, l[0], tmp, cn)
|
|
|
452 else:
|
|
|
453 return False
|
|
|
454 l = l[2:]
|
|
|
455 elif isinstance(l[0], list) and l[1] in ['and', 'or']:
|
|
|
456 if isinstance(l[2], (float, int)) or l[2] == None:
|
|
|
457 tmp = control(None, l[0], cn)
|
|
|
458 if tmp is False:
|
|
|
459 return False
|
|
|
460 else:
|
|
|
461 ris = compute(tmp, l[1], l[2], cn)
|
|
|
462 elif isinstance(l[2], list):
|
|
|
463 tmp = control(None, l[0], cn)
|
|
|
464 tmp2 = control(None, l[2], cn)
|
|
|
465 if tmp is False or tmp2 is False:
|
|
|
466 return False
|
|
|
467 else:
|
|
|
468 ris = compute(tmp, l[1], tmp2, cn)
|
|
|
469 else:
|
|
|
470 return False
|
|
|
471 l = l[3:]
|
|
|
472 else:
|
|
|
473 return False
|
|
|
474 return ris
|
|
|
475
|
|
|
476 ############################ gene #############################################
|
|
|
477
|
|
|
478 def data_gene(gene, type_gene, name, gene_custom):
|
|
|
479 args = process_args(sys.argv)
|
|
|
480 for i in range(len(gene)):
|
|
|
481 tmp = gene.iloc[i, 0]
|
|
|
482 if tmp.startswith(' ') or tmp.endswith(' '):
|
|
|
483 gene.iloc[i, 0] = (tmp.lstrip()).rstrip()
|
|
|
484 gene_dup = [item for item, count in
|
|
|
485 collections.Counter(gene[gene.columns[0]]).items() if count > 1]
|
|
|
486 pat_dup = [item for item, count in
|
|
|
487 collections.Counter(list(gene.columns)).items() if count > 1]
|
|
|
488 if gene_dup:
|
|
|
489 if gene_custom == None:
|
|
|
490 if args.rules_selector == 'HMRcore':
|
|
|
491 gene_in_rule = pk.load(open(args.tool_dir +
|
|
|
492 '/local/HMRcore_genes.p', 'rb'))
|
|
|
493 elif args.rules_selector == 'Recon':
|
|
|
494 gene_in_rule = pk.load(open(args.tool_dir +
|
|
|
495 '/local/Recon_genes.p', 'rb'))
|
|
|
496 gene_in_rule = gene_in_rule.get(type_gene)
|
|
|
497 else:
|
|
|
498 gene_in_rule = gene_custom
|
|
|
499 tmp = []
|
|
|
500 for i in gene_dup:
|
|
|
501 if gene_in_rule.get(i) == 'ok':
|
|
|
502 tmp.append(i)
|
|
|
503 if tmp:
|
|
|
504 sys.exit('Execution aborted because gene ID '
|
|
|
505 + str(tmp) + ' in ' + name + ' is duplicated\n')
|
|
|
506 if pat_dup:
|
|
|
507 sys.exit('Execution aborted: duplicated label\n'
|
|
|
508 + str(pat_dup) + 'in ' + name + '\n')
|
|
|
509 return (gene.set_index(gene.columns[0])).to_dict()
|
|
|
510
|
|
|
511 ############################ resolve ##########################################
|
|
|
512
|
|
|
513 def resolve(genes, rules, ids, resolve_none, name):
|
|
|
514 resolve_rules = {}
|
|
|
515 not_found = []
|
|
|
516 flag = False
|
|
|
517 for key, value in genes.items():
|
|
|
518 tmp_resolve = []
|
|
|
519 for i in range(len(rules)):
|
|
|
520 tmp = rules[i]
|
|
|
521 if tmp:
|
|
|
522 tmp, err = replace_gene_value(tmp, value)
|
|
|
523 if err:
|
|
|
524 not_found.extend(err)
|
|
|
525 ris = control(None, tmp, resolve_none)
|
|
|
526 if ris is False or ris == None:
|
|
|
527 tmp_resolve.append(None)
|
|
|
528 else:
|
|
|
529 tmp_resolve.append(ris)
|
|
|
530 flag = True
|
|
|
531 else:
|
|
|
532 tmp_resolve.append(None)
|
|
|
533 resolve_rules[key] = tmp_resolve
|
|
|
534 if flag is False:
|
|
|
535 sys.exit('Execution aborted: no computable score' +
|
|
|
536 ' (due to missing gene values) for class '
|
|
|
537 + name + ', the class has been disregarded\n')
|
|
|
538 return (resolve_rules, list(set(not_found)))
|
|
|
539
|
|
|
540 ################################# clustering ##################################
|
|
|
541
|
|
|
542 def f_cluster(resolve_rules):
|
|
|
543 os.makedirs('cluster_out')
|
|
|
544 args = process_args(sys.argv)
|
|
13
|
545 k_min = args.k_min
|
|
|
546 k_max = args.k_max
|
|
|
547 if k_min > k_max:
|
|
|
548 warning('k range boundaries inverted.\n')
|
|
|
549 tmp = k_min
|
|
|
550 k_min = k_max
|
|
|
551 k_max = tmp
|
|
|
552 else:
|
|
|
553 warning('k range correct.\n')
|
|
0
|
554 cluster_data = pd.DataFrame.from_dict(resolve_rules, orient = 'index')
|
|
|
555 for i in cluster_data.columns:
|
|
|
556 tmp = cluster_data[i][0]
|
|
|
557 if tmp == None:
|
|
|
558 cluster_data = cluster_data.drop(columns=[i])
|
|
|
559 distorsion = []
|
|
13
|
560 for i in range(k_min, k_max+1):
|
|
0
|
561 tmp_kmeans = KMeans(n_clusters = i,
|
|
|
562 n_init = 100,
|
|
|
563 max_iter = 300,
|
|
|
564 random_state = 0).fit(cluster_data)
|
|
|
565 distorsion.append(tmp_kmeans.inertia_)
|
|
|
566 predict = tmp_kmeans.predict(cluster_data)
|
|
|
567 predict = [x+1 for x in predict]
|
|
7
|
568 classe = (pd.DataFrame(list(zip(cluster_data.index, predict)))).astype(str)
|
|
0
|
569 dest = 'cluster_out/K=' + str(i) + '_' + args.name+'.tsv'
|
|
|
570 classe.to_csv(dest, sep = '\t', index = False,
|
|
|
571 header = ['Patient_ID', 'Class'])
|
|
|
572 plt.figure(0)
|
|
13
|
573 plt.plot(range(k_min, k_max+1), distorsion, marker = 'o')
|
|
0
|
574 plt.xlabel('Number of cluster')
|
|
|
575 plt.ylabel('Distorsion')
|
|
|
576 plt.savefig(args.elbow, dpi = 240, format = 'pdf')
|
|
|
577 if args.cond_hier == 'yes':
|
|
|
578 import scipy.cluster.hierarchy as hier
|
|
|
579 lin = hier.linkage(cluster_data, args.linkage)
|
|
|
580 plt.figure(1)
|
|
|
581 plt.figure(figsize=(10, 5))
|
|
|
582 hier.dendrogram(lin, leaf_font_size = 2, labels = cluster_data.index)
|
|
|
583 plt.savefig(args.dendro, dpi = 480, format = 'pdf')
|
|
|
584 return None
|
|
|
585
|
|
|
586 ################################# main ########################################
|
|
|
587
|
|
|
588 def main():
|
|
|
589 args = process_args(sys.argv)
|
|
|
590 if args.rules_selector == 'HMRcore':
|
|
|
591 recon = pk.load(open(args.tool_dir + '/local/HMRcore_rules.p', 'rb'))
|
|
|
592 elif args.rules_selector == 'Recon':
|
|
|
593 recon = pk.load(open(args.tool_dir + '/local/Recon_rules.p', 'rb'))
|
|
|
594 elif args.rules_selector == 'Custom':
|
|
|
595 ids, rules, gene_in_rule = make_recon(args.custom)
|
|
|
596 resolve_none = check_bool(args.none)
|
|
|
597 dataset = read_dataset(args.data, args.name)
|
|
|
598 dataset.iloc[:, 0] = (dataset.iloc[:, 0]).astype(str)
|
|
|
599 type_gene = gene_type(dataset.iloc[0, 0], args.name)
|
|
|
600 if args.rules_selector != 'Custom':
|
|
|
601 genes = data_gene(dataset, type_gene, args.name, None)
|
|
|
602 ids, rules = load_id_rules(recon.get(type_gene))
|
|
|
603 elif args.rules_selector == 'Custom':
|
|
|
604 genes = data_gene(dataset, type_gene, args.name, gene_in_rule)
|
|
|
605 resolve_rules, err = resolve(genes, rules, ids, resolve_none, args.name)
|
|
|
606 if err:
|
|
|
607 warning('WARNING: gene\n' + str(err) + '\nnot found in class '
|
|
|
608 + args.name + ', the expression level for this gene ' +
|
|
|
609 'will be considered NaN\n')
|
|
|
610 f_cluster(resolve_rules)
|
|
|
611 warning('Execution succeeded')
|
|
|
612 return None
|
|
|
613
|
|
|
614 ###############################################################################
|
|
|
615
|
|
|
616 if __name__ == "__main__":
|
|
6
|
617 main()
|
