-
Notifications
You must be signed in to change notification settings - Fork 0
/
celltype-analysis_topo.py
435 lines (365 loc) · 18.5 KB
/
celltype-analysis_topo.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
# analyse 10x nested x-val data by cell type:
import numpy as np
import pandas as pd
from config import Config
from matplotlib import pyplot as plt
from DBM_toolbox.data_manipulation.data_utils import pickle_objects, unpickle_objects
from sklearn.metrics import roc_curve, auc
import seaborn as sns
import glob
from pathlib import Path
###############################################################################
### 1 Tables with performances per cell type
#measures = ['EIGENVECTOR','BETWEENNESS','CLOSENESS', 'PAGERANK','AVNEIGHBOUR', 'HARMONIC', 'INFORMATION', 'CONSTRAINT', 'ECCENTRICITY',
# 'SUBGRAPH', 'APPROXCURRBET','CLIQUENO','SQUARECLUSTERING','DEGREECENT' ]
#,'BETWEENNESS','CLOSENESS', 'PAGERANK','AVNEIGHBOUR', 'HARMONIC', 'INFORMATION', 'CONSTRAINT', 'ECCENTRICITY',
# 'SUBGRAPH', 'APPROXCURRBET','CLIQUENO','SQUARECLUSTERING','DEGREECENT'
#disc = ['DISCRETIZED']
#base_omics = ['RPPA'] #'RPPA','RNA',
#combinations = base_omics + measures
#name_str = [x[0:3] for x in combinations]
#comb_name = ''.join([y for y in name_str])
#results_folder = f'{comb_name}'
#Path(f"{results_folder}/").mkdir(parents=True, exist_ok=True)
############################
######
def generate_folder_comb(res_filename:str, additional_string=None):
filename = res_filename
fi1 = filename.split('_')
fi2 = fi1[0]
x = len(fi2)
y = [fi2[i:i+3] for i in range(0,x,3)]
base_omics_dict = dict()
measures_dict = dict()
disc_dict = dict()
base_omics_dict['RNA'] = 'RNA'
base_omics_dict['RPP'] = 'RPPA'
base_omics_dict['DNA'] = 'DNA'
measures_dict['EIG'] = 'EIGENVECTOR'
measures_dict['BET'] = 'BETWEENNESS'
measures_dict['CLO'] = 'CLOSENESS'
measures_dict['PAG'] = 'PAGERANK'
measures_dict['AVN'] = 'AVNEIGHBOUR'
measures_dict['HAR'] = 'HARMONIC'
measures_dict['INF'] = 'INFORMATION'
measures_dict['CON'] = 'CONSTRAINT'
measures_dict['ECC'] = 'ECCENTRICITY'
measures_dict['SUB'] = 'SUBGRAPH'
measures_dict['APP'] = 'APPROXCURRBET'
measures_dict['CLI'] = 'CLIQUENO'
measures_dict['SQU'] = 'SQUARECLUSTERING'
measures_dict['DEG'] = 'DEGREECENT'
disc_dict['DIS'] = 'DISCRETIZED'
base_omics_list = []
measures_list = []
disc_list = []
for xi in y:
if xi in base_omics_dict:
base_omics_list.append(base_omics_dict[xi])
elif xi in measures_dict:
measures_list.append(measures_dict[xi])
elif xi in disc_dict:
disc_list.append(disc_dict[xi])
else:
pass
combinations = base_omics_list + measures_list + disc_list
if additional_string is not None:
results_folder = fi1[0] + additional_string
else:
results_folder = fi1[0]
Path(f"{results_folder}/").mkdir(parents=True, exist_ok=True)
return combinations, results_folder
def final_balanced_acc_plot(final_results, comb_name, results_folder):
drugs_list = list(final_results.keys())
tumors_list = [
"NERVOUS",
"SKIN",
"LUNG",
"BREAST",
"LARGE_INTESTINE",
]
exp = [x + '_N' for x in drugs_list]
df_all = pd.DataFrame(index=tumors_list, columns=drugs_list + exp)
df_pos = pd.DataFrame(index=tumors_list, columns=drugs_list + exp)
df_neg = pd.DataFrame(index=tumors_list, columns=drugs_list + exp)
df_prec = pd.DataFrame(index=tumors_list, columns=drugs_list + exp)
df_recall = pd.DataFrame(index=tumors_list, columns=drugs_list + exp)
df_ba = pd.DataFrame(index=tumors_list, columns=drugs_list + exp)
for drug in drugs_list:
for tumor in tumors_list:
print(f'drug: {drug}, tumor: {tumor}')
res = final_results[drug].loc[:, [drug, 'pred2_RFC']]
res = res[res.index.str.contains(tumor)]
N_pos = sum(res.loc[:, drug] == 1)
N_neg = sum(res.loc[:, drug] == 0)
N_all = N_pos + N_neg
if N_all >= 10:
res['sum'] = res.sum(axis=1)
res['tp'] = res['sum'] >= 1.5
res['tn'] = res['sum'] <= 0.5
res['fp'] = (res['sum'] > 0.5) & (res['sum'] <= 1)
res['fn'] = (res['sum'] > 1) & (res['sum'] <= 1.5)
res['pp'] = res['pred2_RFC'] > 0.5
res['pn'] = res['pred2_RFC'] < 0.5
df_pos.loc[tumor, drug +'_N'] = N_pos
df_neg.loc[tumor, drug +'_N'] = N_neg
df_all.loc[tumor, drug +'_N'] = N_all
if N_pos > 0:
df_pos.loc[tumor, drug] = res.loc[:, 'tp'].sum().sum() / N_pos
if N_neg > 0:
df_neg.loc[tumor, drug] = res.loc[:, 'tn'].sum().sum() / N_neg
if res.loc[:, 'pp'].astype(int).sum() > 0:
df_prec.loc[tumor, drug] = res.loc[:, 'tp'].sum().sum().astype(int) / (res.loc[:, 'pp'].astype(int)).sum().sum()
df_recall.loc[tumor, drug] = res.loc[:, 'tp'].sum().sum().astype(int) / (res.loc[:, 'tp'].astype(int) + res.loc[:, 'fn'].astype(int)).sum().sum()
df_all.loc[tumor, drug] = res.loc[:, ['tp', 'tn']].sum().sum().astype(int) / N_all
tpr = (res.loc[:, 'tp'].sum().sum().astype(int) / res.loc[:, ['tp', 'fn']].sum().sum().astype(int))
tnr = (res.loc[:, 'tn'].sum().sum().astype(int) / res.loc[:, ['tn', 'fp']].sum().sum().astype(int))
df_ba.loc[tumor, drug] = (tpr + tnr) / 2
df_pos = df_pos.dropna(how='all')
df_neg = df_neg.dropna(how='all')
df_all = df_all.dropna(how='all')
df_prec = df_prec.dropna(how='all')
df_recall = df_recall.dropna(how='all')
df_ba = df_ba.dropna(how='all')
df_ba = df_ba.dropna(axis=1, how='all')
df_pos.to_csv(f'{results_folder}/Topo_{comb_name}_celltype_results_pos.csv')
df_neg.to_csv(f'{results_folder}/Topo_{comb_name}_celltype_results_neg.csv')
df_all.to_csv(f'{results_folder}/Topo_{comb_name}_celltype_results_all.csv')
df_prec.to_csv(f'{results_folder}/Topo_{comb_name}_celltype_results_prec.csv')
df_recall.to_csv(f'{results_folder}/Topo_{comb_name}_celltype_results_recall.csv')
df_ba.to_csv(f'{results_folder}/Topo_{comb_name}_celltype_results_bal-accuracy.csv')
df_ba = df_ba.apply(pd.to_numeric, errors='coerce')
fig, ax = plt.subplots(figsize=(30, 15))
cmap = sns.cm.rocket_r
sns.set(font_scale=1.4)
sns.heatmap(df_ba, cmap=cmap, annot=True, fmt='.2f', linewidths=.1, ax=ax, annot_kws={
'fontsize': 16,
'fontweight': 'bold',
'fontfamily': 'serif'
})
#cmap=sns.color_palette("rocket", as_cmap=True)
plt.xlabel('')
plt.savefig(f'{results_folder}/{comb_name}_balanced_accuracy')
plt.close()
#####################################################################
##################| ROC Curves
def plot_roc_curves(final_results, comb_name, results_folder):
hfont = {'fontname':'Times New Roman'}
targets_names = list(final_results.keys())
fprs = dict()
tprs = dict()
roc_aucs = dict()
for target_name in targets_names:
fpr = dict()
tpr = dict()
roc_auc = dict()
global_models = [x for x in final_results[target_name].columns if x.startswith('pred2')]
plt.subplots(figsize=(8, 8))
for global_model in global_models:
y_score = final_results[target_name].loc[:, global_model]
y_score = y_score[y_score.isna() == False]
y_true = final_results[target_name].loc[:, target_name]
y_true = y_true.loc[y_score.index]
fpr[global_model], tpr[global_model], _ = roc_curve(y_true, y_score)
roc_auc[global_model] = auc(fpr[global_model], tpr[global_model])
plt.plot(fpr[global_model], tpr[global_model], linewidth=3, label=f"AUC: {round(roc_auc[global_model], 3)}")
plt.plot([0, 1], [0, 1], color="black", linestyle="--")
plt.xlim([0.0, 1.05])
plt.ylim([0.0, 1.05])
plt.xlabel("False Positive Rate", fontsize=20, **hfont)
plt.ylabel("True Positive Rate", fontsize=20, **hfont)
plt.title(f"{target_name.split('_')[0]}",fontsize=20, **hfont)
plt.legend(loc="lower right", fontsize=20)
plt.savefig(f'{results_folder}/{comb_name}_{target_name}')
plt.close()
fprs[target_name] = fpr
tprs[target_name] = tpr
roc_aucs[target_name] = roc_auc
roc_aucs_df = pd.DataFrame.from_dict(roc_aucs).T
roc_aucs_df.to_csv(f'{results_folder}/{comb_name}_AUCS.csv')
###################################################################
### clustergrams
def plot_clustergram(features_file, comb_name, results_folder):
#'f_test_topo2_features_2023-04-26-11-20-18-361678.pkl'
fi = unpickle_objects(features_file)
fi_names = fi['Lapatinib_ActArea'][0]['RFC'].index
res_fi = pd.DataFrame(0, columns=fi_names, index=fi.keys())
res_sd = pd.DataFrame(0, columns=fi_names, index=fi.keys())
for target_name, target_dict in fi.items():
these_means = pd.Series()
for loop_number, loop_dict in target_dict.items():
these_means = pd.concat([these_means, loop_dict['RFC']], axis=1)
this_mean = these_means.mean(axis=1)
this_sd = these_means.std(axis=1)
res_fi.loc[target_name, :] = this_mean
res_sd.loc[target_name, :] = this_sd
colnames = res_fi.columns
colnames = ['+'.join(x.split('_')[1:]) for x in colnames]
indexnames = res_fi.index
indexnames = [x.split('_')[0] for x in indexnames]
res_fi.columns = colnames
res_fi.index = indexnames
res_sd.columns = colnames
res_sd.index = indexnames
res_fi.to_csv(f'{results_folder}/{comb_name}_RFC_contributions_table.csv')
res_sd.to_csv(f'{results_folder}/{comb_name}_RFC_contributions_sd_table.csv')
fig, ax = plt.subplots(figsize=(25, 11))
cmap = sns.cm.rocket_r
sns.heatmap(res_fi, cmap = cmap, annot=False, fmt='.2f', ax=ax)
plt.savefig(f'{results_folder}/{comb_name}_FINAL_RFC_contributions')
plt.close()
sns.clustermap(res_fi, cmap=cmap, xticklabels=True, figsize=(25, 11))
#sns.color_palette("rocket", as_cmap=True)
plt.savefig(f'{results_folder}/{comb_name}_RFC_contributions_cluster')
plt.close()
#
# fig, axs = plt.subplots(nrows=len(indexnames), figsize=(10, 50), sharex='all')
#
# for i, (mean_row, std_row) in enumerate(zip(res_fi.iterrows(), res_sd.iterrows())):
# index = mean_row[0]
# means = mean_row[1]
# stds = std_row[1]
#
# axs[i].bar(means.index, means, yerr=stds, capsize=5)
# # axs[i].set_xlabel(colnames)
# axs[i].set_xticklabels(axs[i].get_xticklabels(), rotation=90)
# axs[i].set_ylabel(f'{index}')
#
# fig.tight_layout()
# fig.savefig("FINAL_contributions_bar_plots.png")
# plt.close()
fig, axs = plt.subplots(nrows=23, figsize=(20, 50), sharex=True, sharey=True, gridspec_kw={"bottom": 0.2})
for i, (mean_row, std_row) in enumerate(zip(res_fi.iterrows(), res_sd.iterrows())):
index = mean_row[0]
data = mean_row[1]
errors = std_row[1]
ax = sns.barplot(x=data.index, y=data.values, ax=axs[i], capsize=.2, ci="sd", color="grey")
ax.errorbar(x=data.index, y=data.values, yerr=errors.values, fmt='none', ecolor='black', capsize=2)
axs[i].set_xlabel("")
axs[i].set_ylabel(index)
axs[i].tick_params(axis='x', rotation=90)
axs[0].xaxis.set_ticks_position('top')
axs[0].xaxis.set_label_position('top')
fig.savefig(f'{results_folder}/{comb_name}_contributions_bar_plots.svg')
plt.close()
#############################################################################
def extract_contributions(basemodels_file, datafile, comb_name, results_folder, combinations):
base_models = unpickle_objects(basemodels_file)
dataset = unpickle_objects(datafile)
#measures = ['EIGENVECTOR','BETWEENNESS','CLOSENESS', 'PAGERANK','AVNEIGHBOUR', 'HARMONIC', 'INFORMATION', 'CONSTRAINT', 'ECCENTRICITY',
# 'SUBGRAPH', 'APPROXCURRBET','CLIQUENO' ]
#disc = ['DISCRETIZED']
#base_omics = ['RPPA','RNA', 'DNA' ]
omics_list = combinations
features_names_dict = {}
for omic in omics_list:
features_names_dict[omic] = list(dataset.to_pandas(omic=omic).columns)
all_features_names = list(dataset.extract(omics_list=omics_list).dataframe.columns)
targets_names = base_models.keys()
fi_dict = {}
for target_name, target_dict in base_models.items():
rr = pd.DataFrame(columns=all_features_names)
idx = 0
for loop1_number, loop1_dict in target_dict.items():
for loop2_number, loop2_dict in loop1_dict.items():
for omic_name, omic_dict in loop2_dict.items():
for algo_name, algo_data in omic_dict.items():
if not algo_name in ['Ridge', 'KNN', 'SVC']:
if algo_name in ['SVM', 'Logistic', 'EN']:
algo_data = algo_data[0]
print(f'{target_name}/{loop1_number}/{loop2_number}/{omic_name}/{algo_name}')
sorted_indices = np.argsort(-algo_data)
ranks = np.empty_like(sorted_indices)
ranks[sorted_indices] = np.arange(len(algo_data))
rr.loc[idx, features_names_dict[omic_name]] = ranks
idx += 1
rr.loc['average', :] = rr.mean(axis=0)
fi_dict[target_name] = rr
for this_key in fi_dict.keys():
fi_dict[this_key].to_csv(f'{results_folder}/fi_{this_key}_{comb_name}.csv')
fi_dict = {}
# for this_key in fi_dict.keys():
# fi_dict[this_key].to_csv(f'{results_folder}/fi_{this_key}_{comb_name}.csv')
#############################################################
def extract_features(datafile, comb_name, results_folder, combinations):
file_list = glob.glob(f'{results_folder}/fi_*')
#file_list.remove('FI_plots')
#dataset = unpickle_objects('Topo_integration_2023-05-01-06-25-23-229169.pkl')
dataset = unpickle_objects(datafile)
omics_list = combinations
algos_list = ['SVC', 'RFC', 'Logistic', 'EN', 'ET', 'XGB', 'Ada']
omics_biglist = dataset.omic
all_col_meds = pd.DataFrame(index=dataset.dataframe.columns)
for file in file_list:
df = pd.read_csv(file)
drugname = file.split('_')[1]
for omic in omics_list:
fig, axs = plt.subplots(nrows=8, figsize=(8, 90))
these_features = omics_biglist[omics_biglist == omic].index
idxs = [x for x in df.columns if x in these_features]
this_df = df.loc[:, idxs]
this_df = this_df.iloc[:-1, :]
# col_means = this_df.mean()
col_meds = this_df.median()
# all_col_meds = pd.concat([all_col_meds, col_meds], axis=1)
# sorted_cols = col_means.sort_values()
sorted_cols = col_meds.sort_values()
sorted_cols = sorted_cols.index[:30]
sorted_df = this_df[sorted_cols].dropna()
ax = axs[0]
colnames = sorted_df.columns
colnames = [x.split('_ENS')[0].split('_Cautio')[0].split('_nmiR')[0].split('/isoval')[0].split('/tauro')[0] for x in colnames]
sorted_df.columns = colnames
sns.boxplot(data=sorted_df + 1, ax=ax, color='white', orient='h')
sns.set_context('paper')
for line in ax.lines:
if line.get_linestyle() == '-':
line.set_color('black')
ax.set_title('All Algorithms')
ax.set_xticklabels(ax.get_xticks(), rotation=90)
ax.set_ylabel('rank')
for i, algo in enumerate(algos_list):
idxs = list(range(i, this_df.shape[0], 7))
algo_df = this_df.iloc[idxs, :]
col_meds = algo_df.median()
all_col_meds = pd.concat([all_col_meds, col_meds], axis=1)
sorted_cols = col_meds.sort_values()
sorted_cols = sorted_cols.index[:30]
sorted_df = algo_df[sorted_cols].dropna()
ax = axs[i+1]
colnames = sorted_df.columns
colnames = [x.split('_ENS')[0].split('_Cautio')[0].split('_nmiR')[0].split('/isoval')[0].split('/tauro')[0] for x in colnames]
sorted_df.columns = colnames
sns.boxplot(data=sorted_df + 1, ax=ax, color='white', orient='h')
sorted_df.to_csv(f"{results_folder}/{comb_name}_FI_med_{drugname}_{omic}_{algo}")
sns.set_context('paper')
for line in ax.lines:
if line.get_linestyle() == '-':
line.set_color('black')
ax.set_title(algo)
ax.set_xticklabels(ax.get_xticks(), rotation=90)
ax.set_ylabel('rank')
plt.suptitle(f'Distribution of features ranks: {drugname} / {omic}')
plt.tight_layout()
plt.savefig(f'{results_folder}/{comb_name}_FI_med_{drugname}_{omic}.tif', format='tif', dpi=300)
plt.close()
all_col_meds.to_csv(f"{results_folder}/{comb_name}_all_col_meds.csv")
########################################################
### features, models, dataset files
final_results = unpickle_objects('RPPSQU_partial_final_results_2023-05-21-14-01-00-898299.pkl') #EIGBETCLOPAGAVNHARINFCONECCSUBAPPCLISQUDEG_final_results_2023-05-08-05-33-00-213462
config = Config("testall/config_test_topo_integrated.yaml")
base_models = 'RPPSQU_partial_base_models_2023-05-21-14-01-00-903747.pkl' #EIGBETCLOPAGAVNHARINFCONECCSUBAPPCLISQUDEG_base_models_2023-05-08-05-33-00-225061
#dataset = 'Topo_integration_2023-05-06-19-07-06-236917.pkl'
dataset = 'Topo_PARTIAL_2023-05-12-13-32-14-031004.pkl'
datafile = dataset
features_file = 'RPPSQU_partial_features_2023-05-21-14-01-00-900328.pkl' #EIGBETCLOPAGAVNHARINFCONECCSUBAPPCLISQUDEG_features_2023-05-08-05-33-00-220327
basemodels_file = base_models
##### Plots
comb, results_folder = generate_folder_comb(features_file, 'PARTIAL')
comb_name = results_folder
combinations = comb
final_balanced_acc_plot(final_results, comb_name, results_folder)
plot_roc_curves(final_results, comb_name, results_folder)
plot_clustergram(features_file, comb_name, results_folder)
extract_contributions(basemodels_file, datafile, comb_name, results_folder, combinations)
extract_features(datafile, comb_name, results_folder, combinations)