|
410
|
1 import argparse
|
|
|
2 import utils.general_utils as utils
|
|
|
3 from typing import Optional, List
|
|
|
4 import os
|
|
|
5 import numpy as np
|
|
|
6 import pandas as pd
|
|
|
7 import cobra
|
|
|
8 import utils.CBS_backend as CBS_backend
|
|
|
9 from joblib import Parallel, delayed, cpu_count
|
|
|
10 from cobra.sampling import OptGPSampler
|
|
|
11 import sys
|
|
411
|
12 import utils.general_utils as utils
|
|
410
|
13
|
|
|
14
|
|
|
15 ################################# process args ###############################
|
|
|
16 def process_args(args :List[str] = None) -> argparse.Namespace:
|
|
|
17 """
|
|
|
18 Processes command-line arguments.
|
|
|
19
|
|
|
20 Args:
|
|
|
21 args (list): List of command-line arguments.
|
|
|
22
|
|
|
23 Returns:
|
|
|
24 Namespace: An object containing parsed arguments.
|
|
|
25 """
|
|
|
26 parser = argparse.ArgumentParser(usage = '%(prog)s [options]',
|
|
|
27 description = 'process some value\'s')
|
|
|
28
|
|
|
29 parser.add_argument("-mo", "--model_upload", type = str,
|
|
|
30 help = "path to input file with custom rules, if provided")
|
|
|
31
|
|
|
32 parser.add_argument('-ol', '--out_log',
|
|
|
33 help = "Output log")
|
|
|
34
|
|
|
35 parser.add_argument('-td', '--tool_dir',
|
|
|
36 type = str,
|
|
|
37 required = True,
|
|
|
38 help = 'your tool directory')
|
|
|
39
|
|
|
40 parser.add_argument('-in', '--input',
|
|
|
41 required = True,
|
|
|
42 type=str,
|
|
|
43 help = 'inputs bounds')
|
|
|
44
|
|
|
45 parser.add_argument('-ni', '--names',
|
|
|
46 required = True,
|
|
|
47 type=str,
|
|
|
48 help = 'cell names')
|
|
|
49
|
|
|
50 parser.add_argument('-a', '--algorithm',
|
|
|
51 type = str,
|
|
|
52 choices = ['OPTGP', 'CBS'],
|
|
|
53 required = True,
|
|
|
54 help = 'choose sampling algorithm')
|
|
|
55
|
|
|
56 parser.add_argument('-th', '--thinning',
|
|
|
57 type = int,
|
|
|
58 default= 100,
|
|
|
59 required=False,
|
|
|
60 help = 'choose thinning')
|
|
|
61
|
|
|
62 parser.add_argument('-ns', '--n_samples',
|
|
|
63 type = int,
|
|
|
64 required = True,
|
|
|
65 help = 'choose how many samples')
|
|
|
66
|
|
|
67 parser.add_argument('-sd', '--seed',
|
|
|
68 type = int,
|
|
|
69 required = True,
|
|
|
70 help = 'seed')
|
|
|
71
|
|
|
72 parser.add_argument('-nb', '--n_batches',
|
|
|
73 type = int,
|
|
|
74 required = True,
|
|
|
75 help = 'choose how many batches')
|
|
|
76
|
|
|
77 parser.add_argument('-ot', '--output_type',
|
|
|
78 type = str,
|
|
|
79 required = True,
|
|
|
80 help = 'output type')
|
|
|
81
|
|
|
82 parser.add_argument('-ota', '--output_type_analysis',
|
|
|
83 type = str,
|
|
|
84 required = False,
|
|
|
85 help = 'output type analysis')
|
|
|
86
|
|
|
87 parser.add_argument('-idop', '--output_path',
|
|
|
88 type = str,
|
|
|
89 default='flux_simulation',
|
|
|
90 help = 'output path for maps')
|
|
|
91
|
|
|
92 ARGS = parser.parse_args(args)
|
|
|
93 return ARGS
|
|
|
94
|
|
|
95 ########################### warning ###########################################
|
|
|
96 def warning(s :str) -> None:
|
|
|
97 """
|
|
|
98 Log a warning message to an output log file and print it to the console.
|
|
|
99
|
|
|
100 Args:
|
|
|
101 s (str): The warning message to be logged and printed.
|
|
|
102
|
|
|
103 Returns:
|
|
|
104 None
|
|
|
105 """
|
|
|
106 with open(ARGS.out_log, 'a') as log:
|
|
|
107 log.write(s + "\n\n")
|
|
|
108 print(s)
|
|
|
109
|
|
|
110
|
|
|
111 def write_to_file(dataset: pd.DataFrame, name: str, keep_index:bool=False)->None:
|
|
|
112 dataset.index.name = 'Reactions'
|
|
|
113 dataset.to_csv(ARGS.output_path + "/" + name + ".csv", sep = '\t', index = keep_index)
|
|
|
114
|
|
|
115 ############################ dataset input ####################################
|
|
|
116 def read_dataset(data :str, name :str) -> pd.DataFrame:
|
|
|
117 """
|
|
|
118 Read a dataset from a CSV file and return it as a pandas DataFrame.
|
|
|
119
|
|
|
120 Args:
|
|
|
121 data (str): Path to the CSV file containing the dataset.
|
|
|
122 name (str): Name of the dataset, used in error messages.
|
|
|
123
|
|
|
124 Returns:
|
|
|
125 pandas.DataFrame: DataFrame containing the dataset.
|
|
|
126
|
|
|
127 Raises:
|
|
|
128 pd.errors.EmptyDataError: If the CSV file is empty.
|
|
|
129 sys.exit: If the CSV file has the wrong format, the execution is aborted.
|
|
|
130 """
|
|
|
131 try:
|
|
|
132 dataset = pd.read_csv(data, sep = '\t', header = 0, index_col=0, engine='python')
|
|
|
133 except pd.errors.EmptyDataError:
|
|
|
134 sys.exit('Execution aborted: wrong format of ' + name + '\n')
|
|
|
135 if len(dataset.columns) < 2:
|
|
|
136 sys.exit('Execution aborted: wrong format of ' + name + '\n')
|
|
|
137 return dataset
|
|
|
138
|
|
|
139
|
|
|
140
|
|
|
141 def OPTGP_sampler(model:cobra.Model, model_name:str, n_samples:int=1000, thinning:int=100, n_batches:int=1, seed:int=0)-> None:
|
|
|
142 """
|
|
|
143 Samples from the OPTGP (Optimal Global Perturbation) algorithm and saves the results to CSV files.
|
|
|
144
|
|
|
145 Args:
|
|
|
146 model (cobra.Model): The COBRA model to sample from.
|
|
|
147 model_name (str): The name of the model, used in naming output files.
|
|
|
148 n_samples (int, optional): Number of samples per batch. Default is 1000.
|
|
|
149 thinning (int, optional): Thinning parameter for the sampler. Default is 100.
|
|
|
150 n_batches (int, optional): Number of batches to run. Default is 1.
|
|
|
151 seed (int, optional): Random seed for reproducibility. Default is 0.
|
|
|
152
|
|
|
153 Returns:
|
|
|
154 None
|
|
|
155 """
|
|
|
156
|
|
|
157 for i in range(0, n_batches):
|
|
|
158 optgp = OptGPSampler(model, thinning, seed)
|
|
|
159 samples = optgp.sample(n_samples)
|
|
|
160 samples.to_csv(ARGS.output_path + "/" + model_name + '_'+ str(i)+'_OPTGP.csv', index=False)
|
|
|
161 seed+=1
|
|
|
162 samplesTotal = pd.DataFrame()
|
|
|
163 for i in range(0, n_batches):
|
|
|
164 samples_batch = pd.read_csv(ARGS.output_path + "/" + model_name + '_'+ str(i)+'_OPTGP.csv')
|
|
|
165 samplesTotal = pd.concat([samplesTotal, samples_batch], ignore_index = True)
|
|
|
166
|
|
|
167 write_to_file(samplesTotal.T, model_name, True)
|
|
|
168
|
|
|
169 for i in range(0, n_batches):
|
|
|
170 os.remove(ARGS.output_path + "/" + model_name + '_'+ str(i)+'_OPTGP.csv')
|
|
|
171 pass
|
|
|
172
|
|
|
173
|
|
|
174 def CBS_sampler(model:cobra.Model, model_name:str, n_samples:int=1000, n_batches:int=1, seed:int=0)-> None:
|
|
|
175 """
|
|
|
176 Samples using the CBS (Constraint-based Sampling) algorithm and saves the results to CSV files.
|
|
|
177
|
|
|
178 Args:
|
|
|
179 model (cobra.Model): The COBRA model to sample from.
|
|
|
180 model_name (str): The name of the model, used in naming output files.
|
|
|
181 n_samples (int, optional): Number of samples per batch. Default is 1000.
|
|
|
182 n_batches (int, optional): Number of batches to run. Default is 1.
|
|
|
183 seed (int, optional): Random seed for reproducibility. Default is 0.
|
|
|
184
|
|
|
185 Returns:
|
|
|
186 None
|
|
|
187 """
|
|
|
188
|
|
|
189 df_FVA = cobra.flux_analysis.flux_variability_analysis(model,fraction_of_optimum=0).round(6)
|
|
|
190
|
|
|
191 df_coefficients = CBS_backend.randomObjectiveFunction(model, n_samples*n_batches, df_FVA, seed=seed)
|
|
|
192
|
|
|
193 for i in range(0, n_batches):
|
|
|
194 samples = pd.DataFrame(columns =[reaction.id for reaction in model.reactions], index = range(n_samples))
|
|
|
195 try:
|
|
|
196 CBS_backend.randomObjectiveFunctionSampling(model, n_samples, df_coefficients.iloc[:,i*n_samples:(i+1)*n_samples], samples)
|
|
|
197 except Exception as e:
|
|
|
198 utils.logWarning(
|
|
|
199 "Warning: GLPK solver has failed for " + model_name + ". Trying with COBRA interface. Error:" + str(e),
|
|
|
200 ARGS.out_log)
|
|
|
201 CBS_backend.randomObjectiveFunctionSampling_cobrapy(model, n_samples, df_coefficients.iloc[:,i*n_samples:(i+1)*n_samples],
|
|
|
202 samples)
|
|
|
203 utils.logWarning(ARGS.output_path + "/" + model_name + '_'+ str(i)+'_CBS.csv', ARGS.out_log)
|
|
|
204 samples.to_csv(ARGS.output_path + "/" + model_name + '_'+ str(i)+'_CBS.csv', index=False)
|
|
|
205
|
|
|
206 samplesTotal = pd.DataFrame()
|
|
|
207 for i in range(0, n_batches):
|
|
|
208 samples_batch = pd.read_csv(ARGS.output_path + "/" + model_name + '_'+ str(i)+'_CBS.csv')
|
|
|
209 samplesTotal = pd.concat([samplesTotal, samples_batch], ignore_index = True)
|
|
|
210
|
|
|
211 write_to_file(samplesTotal.T, model_name, True)
|
|
|
212
|
|
|
213 for i in range(0, n_batches):
|
|
|
214 os.remove(ARGS.output_path + "/" + model_name + '_'+ str(i)+'_CBS.csv')
|
|
|
215 pass
|
|
|
216
|
|
|
217
|
|
|
218 def model_sampler(model_input_original:cobra.Model, bounds_path:str, cell_name:str)-> List[pd.DataFrame]:
|
|
|
219 """
|
|
|
220 Prepares the model with bounds from the dataset and performs sampling and analysis based on the selected algorithm.
|
|
|
221
|
|
|
222 Args:
|
|
|
223 model_input_original (cobra.Model): The original COBRA model.
|
|
|
224 bounds_path (str): Path to the CSV file containing the bounds dataset.
|
|
|
225 cell_name (str): Name of the cell, used to generate filenames for output.
|
|
|
226
|
|
|
227 Returns:
|
|
|
228 List[pd.DataFrame]: A list of DataFrames containing statistics and analysis results.
|
|
|
229 """
|
|
|
230
|
|
|
231 model_input = model_input_original.copy()
|
|
|
232 bounds_df = read_dataset(bounds_path, "bounds dataset")
|
|
|
233 for rxn_index, row in bounds_df.iterrows():
|
|
|
234 model_input.reactions.get_by_id(rxn_index).lower_bound = row.lower_bound
|
|
|
235 model_input.reactions.get_by_id(rxn_index).upper_bound = row.upper_bound
|
|
|
236
|
|
|
237
|
|
|
238 if ARGS.algorithm == 'OPTGP':
|
|
|
239 OPTGP_sampler(model_input, cell_name, ARGS.n_samples, ARGS.thinning, ARGS.n_batches, ARGS.seed)
|
|
|
240
|
|
|
241 elif ARGS.algorithm == 'CBS':
|
|
|
242 CBS_sampler(model_input, cell_name, ARGS.n_samples, ARGS.n_batches, ARGS.seed)
|
|
|
243
|
|
|
244 df_mean, df_median, df_quantiles = fluxes_statistics(cell_name, ARGS.output_types)
|
|
|
245
|
|
|
246 if("fluxes" not in ARGS.output_types):
|
|
|
247 os.remove(ARGS.output_path + "/" + cell_name + '.csv')
|
|
|
248
|
|
|
249 returnList = []
|
|
|
250 returnList.append(df_mean)
|
|
|
251 returnList.append(df_median)
|
|
|
252 returnList.append(df_quantiles)
|
|
|
253
|
|
|
254 df_pFBA, df_FVA, df_sensitivity = fluxes_analysis(model_input, cell_name, ARGS.output_type_analysis)
|
|
|
255
|
|
|
256 if("pFBA" in ARGS.output_type_analysis):
|
|
|
257 returnList.append(df_pFBA)
|
|
|
258 if("FVA" in ARGS.output_type_analysis):
|
|
|
259 returnList.append(df_FVA)
|
|
|
260 if("sensitivity" in ARGS.output_type_analysis):
|
|
|
261 returnList.append(df_sensitivity)
|
|
|
262
|
|
|
263 return returnList
|
|
|
264
|
|
|
265 def fluxes_statistics(model_name: str, output_types:List)-> List[pd.DataFrame]:
|
|
|
266 """
|
|
|
267 Computes statistics (mean, median, quantiles) for the fluxes.
|
|
|
268
|
|
|
269 Args:
|
|
|
270 model_name (str): Name of the model, used in filename for input.
|
|
|
271 output_types (List[str]): Types of statistics to compute (mean, median, quantiles).
|
|
|
272
|
|
|
273 Returns:
|
|
|
274 List[pd.DataFrame]: List of DataFrames containing mean, median, and quantiles statistics.
|
|
|
275 """
|
|
|
276
|
|
|
277 df_mean = pd.DataFrame()
|
|
|
278 df_median= pd.DataFrame()
|
|
|
279 df_quantiles= pd.DataFrame()
|
|
|
280
|
|
|
281 df_samples = pd.read_csv(ARGS.output_path + "/" + model_name + '.csv', sep = '\t', index_col = 0).T
|
|
|
282 df_samples = df_samples.round(8)
|
|
|
283
|
|
|
284 for output_type in output_types:
|
|
|
285 if(output_type == "mean"):
|
|
|
286 df_mean = df_samples.mean()
|
|
|
287 df_mean = df_mean.to_frame().T
|
|
|
288 df_mean = df_mean.reset_index(drop=True)
|
|
|
289 df_mean.index = [model_name]
|
|
|
290 elif(output_type == "median"):
|
|
|
291 df_median = df_samples.median()
|
|
|
292 df_median = df_median.to_frame().T
|
|
|
293 df_median = df_median.reset_index(drop=True)
|
|
|
294 df_median.index = [model_name]
|
|
|
295 elif(output_type == "quantiles"):
|
|
|
296 newRow = []
|
|
|
297 cols = []
|
|
|
298 for rxn in df_samples.columns:
|
|
|
299 quantiles = df_samples[rxn].quantile([0.25, 0.50, 0.75])
|
|
|
300 newRow.append(quantiles[0.25])
|
|
|
301 cols.append(rxn + "_q1")
|
|
|
302 newRow.append(quantiles[0.5])
|
|
|
303 cols.append(rxn + "_q2")
|
|
|
304 newRow.append(quantiles[0.75])
|
|
|
305 cols.append(rxn + "_q3")
|
|
|
306 df_quantiles = pd.DataFrame(columns=cols)
|
|
|
307 df_quantiles.loc[0] = newRow
|
|
|
308 df_quantiles = df_quantiles.reset_index(drop=True)
|
|
|
309 df_quantiles.index = [model_name]
|
|
|
310
|
|
|
311 return df_mean, df_median, df_quantiles
|
|
|
312
|
|
|
313 def fluxes_analysis(model:cobra.Model, model_name:str, output_types:List)-> List[pd.DataFrame]:
|
|
|
314 """
|
|
|
315 Performs flux analysis including pFBA, FVA, and sensitivity analysis.
|
|
|
316
|
|
|
317 Args:
|
|
|
318 model (cobra.Model): The COBRA model to analyze.
|
|
|
319 model_name (str): Name of the model, used in filenames for output.
|
|
|
320 output_types (List[str]): Types of analysis to perform (pFBA, FVA, sensitivity).
|
|
|
321
|
|
|
322 Returns:
|
|
|
323 List[pd.DataFrame]: List of DataFrames containing pFBA, FVA, and sensitivity analysis results.
|
|
|
324 """
|
|
|
325
|
|
|
326 df_pFBA = pd.DataFrame()
|
|
|
327 df_FVA= pd.DataFrame()
|
|
|
328 df_sensitivity= pd.DataFrame()
|
|
|
329
|
|
|
330 for output_type in output_types:
|
|
|
331 if(output_type == "pFBA"):
|
|
|
332 model.objective = "Biomass"
|
|
|
333 solution = cobra.flux_analysis.pfba(model)
|
|
|
334 fluxes = solution.fluxes
|
|
|
335 df_pFBA.loc[0,[rxn._id for rxn in model.reactions]] = fluxes.tolist()
|
|
|
336 df_pFBA = df_pFBA.reset_index(drop=True)
|
|
|
337 df_pFBA.index = [model_name]
|
|
|
338 df_pFBA = df_pFBA.astype(float).round(6)
|
|
|
339 elif(output_type == "FVA"):
|
|
|
340 fva = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0, processes=1).round(8)
|
|
|
341 columns = []
|
|
|
342 for rxn in fva.index.to_list():
|
|
|
343 columns.append(rxn + "_min")
|
|
|
344 columns.append(rxn + "_max")
|
|
|
345 df_FVA= pd.DataFrame(columns = columns)
|
|
|
346 for index_rxn, row in fva.iterrows():
|
|
|
347 df_FVA.loc[0, index_rxn+ "_min"] = fva.loc[index_rxn, "minimum"]
|
|
|
348 df_FVA.loc[0, index_rxn+ "_max"] = fva.loc[index_rxn, "maximum"]
|
|
|
349 df_FVA = df_FVA.reset_index(drop=True)
|
|
|
350 df_FVA.index = [model_name]
|
|
|
351 df_FVA = df_FVA.astype(float).round(6)
|
|
|
352 elif(output_type == "sensitivity"):
|
|
|
353 model.objective = "Biomass"
|
|
|
354 solution_original = model.optimize().objective_value
|
|
|
355 reactions = model.reactions
|
|
|
356 single = cobra.flux_analysis.single_reaction_deletion(model)
|
|
|
357 newRow = []
|
|
|
358 df_sensitivity = pd.DataFrame(columns = [rxn.id for rxn in reactions], index = [model_name])
|
|
|
359 for rxn in reactions:
|
|
|
360 newRow.append(single.knockout[rxn.id].growth.values[0]/solution_original)
|
|
|
361 df_sensitivity.loc[model_name] = newRow
|
|
|
362 df_sensitivity = df_sensitivity.astype(float).round(6)
|
|
|
363 return df_pFBA, df_FVA, df_sensitivity
|
|
|
364
|
|
|
365 ############################# main ###########################################
|
|
|
366 def main(args :List[str] = None) -> None:
|
|
|
367 """
|
|
|
368 Initializes everything and sets the program in motion based on the fronted input arguments.
|
|
|
369
|
|
|
370 Returns:
|
|
|
371 None
|
|
|
372 """
|
|
|
373
|
|
|
374 num_processors = cpu_count()
|
|
|
375
|
|
|
376 global ARGS
|
|
|
377 ARGS = process_args(args)
|
|
|
378
|
|
|
379 if not os.path.exists(ARGS.output_path):
|
|
|
380 os.makedirs(ARGS.output_path)
|
|
|
381
|
|
|
382 #model_type :utils.Model = ARGS.model_selector
|
|
|
383 #if model_type is utils.Model.Custom:
|
|
|
384 # model = model_type.getCOBRAmodel(customPath = utils.FilePath.fromStrPath(ARGS.model), customExtension = utils.FilePath.fromStrPath(ARGS.model_name).ext)
|
|
|
385 #else:
|
|
|
386 # model = model_type.getCOBRAmodel(toolDir=ARGS.tool_dir)
|
|
|
387
|
|
|
388 model = utils.build_cobra_model_from_csv(ARGS.model_upload)
|
|
|
389
|
|
|
390 validation = utils.validate_model(model)
|
|
|
391
|
|
|
392 print("\n=== VALIDAZIONE MODELLO ===")
|
|
|
393 for key, value in validation.items():
|
|
|
394 print(f"{key}: {value}")
|
|
|
395
|
|
|
396 #Set solver verbosity to 1 to see warning and error messages only.
|
|
|
397 model.solver.configuration.verbosity = 1
|
|
|
398
|
|
|
399 ARGS.bounds = ARGS.input.split(",")
|
|
|
400 ARGS.bounds_name = ARGS.names.split(",")
|
|
|
401 ARGS.output_types = ARGS.output_type.split(",")
|
|
|
402 ARGS.output_type_analysis = ARGS.output_type_analysis.split(",")
|
|
|
403
|
|
|
404
|
|
|
405 results = Parallel(n_jobs=num_processors)(delayed(model_sampler)(model, bounds_path, cell_name) for bounds_path, cell_name in zip(ARGS.bounds, ARGS.bounds_name))
|
|
|
406
|
|
|
407 all_mean = pd.concat([result[0] for result in results], ignore_index=False)
|
|
|
408 all_median = pd.concat([result[1] for result in results], ignore_index=False)
|
|
|
409 all_quantiles = pd.concat([result[2] for result in results], ignore_index=False)
|
|
|
410
|
|
|
411 if("mean" in ARGS.output_types):
|
|
|
412 all_mean = all_mean.fillna(0.0)
|
|
|
413 all_mean = all_mean.sort_index()
|
|
|
414 write_to_file(all_mean.T, "mean", True)
|
|
|
415
|
|
|
416 if("median" in ARGS.output_types):
|
|
|
417 all_median = all_median.fillna(0.0)
|
|
|
418 all_median = all_median.sort_index()
|
|
|
419 write_to_file(all_median.T, "median", True)
|
|
|
420
|
|
|
421 if("quantiles" in ARGS.output_types):
|
|
|
422 all_quantiles = all_quantiles.fillna(0.0)
|
|
|
423 all_quantiles = all_quantiles.sort_index()
|
|
|
424 write_to_file(all_quantiles.T, "quantiles", True)
|
|
|
425
|
|
|
426 index_result = 3
|
|
|
427 if("pFBA" in ARGS.output_type_analysis):
|
|
|
428 all_pFBA = pd.concat([result[index_result] for result in results], ignore_index=False)
|
|
|
429 all_pFBA = all_pFBA.sort_index()
|
|
|
430 write_to_file(all_pFBA.T, "pFBA", True)
|
|
|
431 index_result+=1
|
|
|
432 if("FVA" in ARGS.output_type_analysis):
|
|
|
433 all_FVA= pd.concat([result[index_result] for result in results], ignore_index=False)
|
|
|
434 all_FVA = all_FVA.sort_index()
|
|
|
435 write_to_file(all_FVA.T, "FVA", True)
|
|
|
436 index_result+=1
|
|
|
437 if("sensitivity" in ARGS.output_type_analysis):
|
|
|
438 all_sensitivity = pd.concat([result[index_result] for result in results], ignore_index=False)
|
|
|
439 all_sensitivity = all_sensitivity.sort_index()
|
|
|
440 write_to_file(all_sensitivity.T, "sensitivity", True)
|
|
|
441
|
|
|
442 pass
|
|
|
443
|
|
|
444 ##############################################################################
|
|
|
445 if __name__ == "__main__":
|
|
|
446 main() |
