Mercurial > repos > bimib > marea

comparison marea_cluster.py @ 38:4e1b466935cd draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bimib
date Mon, 25 Nov 2019 12:02:28 -0500
parents
children
comparison
equal deleted inserted replaced
37:2495c7772ca8 38:4e1b466935cd
1 # -*- coding: utf-8 -*-
2 """
3 Created on Mon Jun 3 19:51:00 2019
4 @author: Narger
5 """
6
7 import sys
8 import argparse
9 import os
10 from sklearn.datasets import make_blobs
11 from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
12 from sklearn.metrics import silhouette_samples, silhouette_score, davies_bouldin_score, cluster
13 import matplotlib
14 matplotlib.use('agg')
15 import matplotlib.pyplot as plt
16 import scipy.cluster.hierarchy as shc
17 import matplotlib.cm as cm
18 import numpy as np
19 import pandas as pd
20
21 ################################# process args ###############################
22
23 def process_args(args):
24 parser = argparse.ArgumentParser(usage = '%(prog)s [options]',
25 description = 'process some value\'s' +
26 ' genes to create class.')
27
28 parser.add_argument('-ol', '--out_log',
29 help = "Output log")
30
31 parser.add_argument('-in', '--input',
32 type = str,
33 help = 'input dataset')
34
35 parser.add_argument('-cy', '--cluster_type',
36 type = str,
37 choices = ['kmeans', 'dbscan', 'hierarchy'],
38 default = 'kmeans',
39 help = 'choose clustering algorythm')
40
41 parser.add_argument('-k1', '--k_min',
42 type = int,
43 default = 2,
44 help = 'choose minimun cluster number to be generated')
45
46 parser.add_argument('-k2', '--k_max',
47 type = int,
48 default = 7,
49 help = 'choose maximum cluster number to be generated')
50
51 parser.add_argument('-el', '--elbow',
52 type = str,
53 default = 'false',
54 choices = ['true', 'false'],
55 help = 'choose if you want to generate an elbow plot for kmeans')
56
57 parser.add_argument('-si', '--silhouette',
58 type = str,
59 default = 'false',
60 choices = ['true', 'false'],
61 help = 'choose if you want silhouette plots')
62
63 parser.add_argument('-td', '--tool_dir',
64 type = str,
65 required = True,
66 help = 'your tool directory')
67
68 parser.add_argument('-ms', '--min_samples',
69 type = float,
70 help = 'min samples for dbscan (optional)')
71
72 parser.add_argument('-ep', '--eps',
73 type = float,
74 help = 'eps for dbscan (optional)')
75
76 parser.add_argument('-bc', '--best_cluster',
77 type = str,
78 help = 'output of best cluster tsv')
79
80
81
82 args = parser.parse_args()
83 return args
84
85 ########################### warning ###########################################
86
87 def warning(s):
88 args = process_args(sys.argv)
89 with open(args.out_log, 'a') as log:
90 log.write(s + "\n\n")
91 print(s)
92
93 ########################## read dataset ######################################
94
95 def read_dataset(dataset):
96 try:
97 dataset = pd.read_csv(dataset, sep = '\t', header = 0)
98 except pd.errors.EmptyDataError:
99 sys.exit('Execution aborted: wrong format of dataset\n')
100 if len(dataset.columns) < 2:
101 sys.exit('Execution aborted: wrong format of dataset\n')
102 return dataset
103
104 ############################ rewrite_input ###################################
105
106 def rewrite_input(dataset):
107 #Riscrivo il dataset come dizionario di liste,
108 #non come dizionario di dizionari
109
110 dataset.pop('Reactions', None)
111
112 for key, val in dataset.items():
113 l = []
114 for i in val:
115 if i == 'None':
116 l.append(None)
117 else:
118 l.append(float(i))
119
120 dataset[key] = l
121
122 return dataset
123
124 ############################## write to csv ##################################
125
126 def write_to_csv (dataset, labels, name):
127 #labels = predict
128 predict = [x+1 for x in labels]
129
130 classe = (pd.DataFrame(list(zip(dataset.index, predict)))).astype(str)
131
132 dest = name
133 classe.to_csv(dest, sep = '\t', index = False,
134 header = ['Patient_ID', 'Class'])
135
136 ########################### trova il massimo in lista ########################
137 def max_index (lista):
138 best = -1
139 best_index = 0
140 for i in range(len(lista)):
141 if lista[i] > best:
142 best = lista [i]
143 best_index = i
144
145 return best_index
146
147 ################################ kmeans #####################################
148
149 def kmeans (k_min, k_max, dataset, elbow, silhouette, best_cluster):
150 if not os.path.exists('clustering'):
151 os.makedirs('clustering')
152
153
154 if elbow == 'true':
155 elbow = True
156 else:
157 elbow = False
158
159 if silhouette == 'true':
160 silhouette = True
161 else:
162 silhouette = False
163
164 range_n_clusters = [i for i in range(k_min, k_max+1)]
165 distortions = []
166 scores = []
167 all_labels = []
168
169 clusterer = KMeans(n_clusters=1, random_state=10)
170 distortions.append(clusterer.fit(dataset).inertia_)
171
172
173 for n_clusters in range_n_clusters:
174 clusterer = KMeans(n_clusters=n_clusters, random_state=10)
175 cluster_labels = clusterer.fit_predict(dataset)
176
177 all_labels.append(cluster_labels)
178 if n_clusters == 1:
179 silhouette_avg = 0
180 else:
181 silhouette_avg = silhouette_score(dataset, cluster_labels)
182 scores.append(silhouette_avg)
183 distortions.append(clusterer.fit(dataset).inertia_)
184
185 best = max_index(scores) + k_min
186
187 for i in range(len(all_labels)):
188 prefix = ''
189 if (i + k_min == best):
190 prefix = '_BEST'
191
192 write_to_csv(dataset, all_labels[i], 'clustering/kmeans_with_' + str(i + k_min) + prefix + '_clusters.tsv')
193
194
195 if (prefix == '_BEST'):
196 labels = all_labels[i]
197 predict = [x+1 for x in labels]
198 classe = (pd.DataFrame(list(zip(dataset.index, predict)))).astype(str)
199 classe.to_csv(best_cluster, sep = '\t', index = False, header = ['Patient_ID', 'Class'])
200
201
202
203
204 if silhouette:
205 silihouette_draw(dataset, all_labels[i], i + k_min, 'clustering/silhouette_with_' + str(i + k_min) + prefix + '_clusters.png')
206
207
208 if elbow:
209 elbow_plot(distortions, k_min,k_max)
210
211
212
213
214
215 ############################## elbow_plot ####################################
216
217 def elbow_plot (distortions, k_min, k_max):
218 plt.figure(0)
219 x = list(range(k_min, k_max + 1))
220 x.insert(0, 1)
221 plt.plot(x, distortions, marker = 'o')
222 plt.xlabel('Number of clusters (k)')
223 plt.ylabel('Distortion')
224 s = 'clustering/elbow_plot.png'
225 fig = plt.gcf()
226 fig.set_size_inches(18.5, 10.5, forward = True)
227 fig.savefig(s, dpi=100)
228
229
230 ############################## silhouette plot ###############################
231 def silihouette_draw(dataset, labels, n_clusters, path):
232 if n_clusters == 1:
233 return None
234
235 silhouette_avg = silhouette_score(dataset, labels)
236 warning("For n_clusters = " + str(n_clusters) +
237 " The average silhouette_score is: " + str(silhouette_avg))
238
239 plt.close('all')
240 # Create a subplot with 1 row and 2 columns
241 fig, (ax1) = plt.subplots(1, 1)
242
243 fig.set_size_inches(18, 7)
244
245 # The 1st subplot is the silhouette plot
246 # The silhouette coefficient can range from -1, 1 but in this example all
247 # lie within [-0.1, 1]
248 ax1.set_xlim([-1, 1])
249 # The (n_clusters+1)*10 is for inserting blank space between silhouette
250 # plots of individual clusters, to demarcate them clearly.
251 ax1.set_ylim([0, len(dataset) + (n_clusters + 1) * 10])
252
253 # Compute the silhouette scores for each sample
254 sample_silhouette_values = silhouette_samples(dataset, labels)
255
256 y_lower = 10
257 for i in range(n_clusters):
258 # Aggregate the silhouette scores for samples belonging to
259 # cluster i, and sort them
260 ith_cluster_silhouette_values = \
261 sample_silhouette_values[labels == i]
262
263 ith_cluster_silhouette_values.sort()
264
265 size_cluster_i = ith_cluster_silhouette_values.shape[0]
266 y_upper = y_lower + size_cluster_i
267
268 color = cm.nipy_spectral(float(i) / n_clusters)
269 ax1.fill_betweenx(np.arange(y_lower, y_upper),
270 0, ith_cluster_silhouette_values,
271 facecolor=color, edgecolor=color, alpha=0.7)
272
273 # Label the silhouette plots with their cluster numbers at the middle
274 ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))
275
276 # Compute the new y_lower for next plot
277 y_lower = y_upper + 10 # 10 for the 0 samples
278
279 ax1.set_title("The silhouette plot for the various clusters.")
280 ax1.set_xlabel("The silhouette coefficient values")
281 ax1.set_ylabel("Cluster label")
282
283 # The vertical line for average silhouette score of all the values
284 ax1.axvline(x=silhouette_avg, color="red", linestyle="--")
285
286 ax1.set_yticks([]) # Clear the yaxis labels / ticks
287 ax1.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])
288
289
290 plt.suptitle(("Silhouette analysis for clustering on sample data "
291 "with n_clusters = " + str(n_clusters) + "\nAverage silhouette_score = " + str(silhouette_avg)), fontsize=12, fontweight='bold')
292
293
294 plt.savefig(path, bbox_inches='tight')
295
296 ######################## dbscan ##############################################
297
298 def dbscan(dataset, eps, min_samples, best_cluster):
299 if not os.path.exists('clustering'):
300 os.makedirs('clustering')
301
302 if eps is not None:
303 clusterer = DBSCAN(eps = eps, min_samples = min_samples)
304 else:
305 clusterer = DBSCAN()
306
307 clustering = clusterer.fit(dataset)
308
309 core_samples_mask = np.zeros_like(clustering.labels_, dtype=bool)
310 core_samples_mask[clustering.core_sample_indices_] = True
311 labels = clustering.labels_
312
313 # Number of clusters in labels, ignoring noise if present.
314 n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
315
316
317 labels = labels
318 predict = [x+1 for x in labels]
319 classe = (pd.DataFrame(list(zip(dataset.index, predict)))).astype(str)
320 classe.to_csv(best_cluster, sep = '\t', index = False, header = ['Patient_ID', 'Class'])
321
322
323 ########################## hierachical #######################################
324
325 def hierachical_agglomerative(dataset, k_min, k_max, best_cluster, silhouette):
326
327 if not os.path.exists('clustering'):
328 os.makedirs('clustering')
329
330 plt.figure(figsize=(10, 7))
331 plt.title("Classes Dendogram")
332 shc.dendrogram(shc.linkage(dataset, method='ward'), labels=dataset.index.values.tolist())
333 fig = plt.gcf()
334 fig.savefig('clustering/dendogram.png', dpi=200)
335
336 range_n_clusters = [i for i in range(k_min, k_max+1)]
337
338 scores = []
339 labels = []
340
341 for n_clusters in range_n_clusters:
342 cluster = AgglomerativeClustering(n_clusters=n_clusters, affinity='euclidean', linkage='ward')
343 cluster.fit_predict(dataset)
344 cluster_labels = cluster.labels_
345 labels.append(cluster_labels)
346 write_to_csv(dataset, cluster_labels, 'clustering/hierarchical_with_' + str(n_clusters) + '_clusters.tsv')
347
348 best = max_index(scores) + k_min
349
350 for i in range(len(labels)):
351 prefix = ''
352 if (i + k_min == best):
353 prefix = '_BEST'
354 if silhouette == 'true':
355 silihouette_draw(dataset, labels[i], i + k_min, 'clustering/silhouette_with_' + str(i + k_min) + prefix + '_clusters.png')
356
357 for i in range(len(labels)):
358 if (i + k_min == best):
359 labels = labels[i]
360 predict = [x+1 for x in labels]
361 classe = (pd.DataFrame(list(zip(dataset.index, predict)))).astype(str)
362 classe.to_csv(best_cluster, sep = '\t', index = False, header = ['Patient_ID', 'Class'])
363
364
365 ############################# main ###########################################
366
367
368 def main():
369 if not os.path.exists('clustering'):
370 os.makedirs('clustering')
371
372 args = process_args(sys.argv)
373
374 #Data read
375
376 X = read_dataset(args.input)
377 X = pd.DataFrame.to_dict(X, orient='list')
378 X = rewrite_input(X)
379 X = pd.DataFrame.from_dict(X, orient = 'index')
380
381 for i in X.columns:
382 tmp = X[i][0]
383 if tmp == None:
384 X = X.drop(columns=[i])
385
386
387 if args.cluster_type == 'kmeans':
388 kmeans(args.k_min, args.k_max, X, args.elbow, args.silhouette, args.best_cluster)
389
390 if args.cluster_type == 'dbscan':
391 dbscan(X, args.eps, args.min_samples, args.best_cluster)
392
393 if args.cluster_type == 'hierarchy':
394 hierachical_agglomerative(X, args.k_min, args.k_max, args.best_cluster, args.silhouette)
395
396 ##############################################################################
397
398 if __name__ == "__main__":
399 main()