|
489
|
1 """
|
|
|
2 Flux sampling and analysis utilities for COBRA models.
|
|
|
3
|
|
|
4 This script supports two modes:
|
|
|
5 - Mode 1 (model_and_bounds=True): load a base model and apply bounds from
|
|
|
6 separate files before sampling.
|
|
|
7 - Mode 2 (model_and_bounds=False): load complete models and sample directly.
|
|
|
8
|
|
|
9 Sampling algorithms supported: OPTGP and CBS. Outputs include flux samples
|
|
|
10 and optional analyses (pFBA, FVA, sensitivity), saved as tabular files.
|
|
|
11 """
|
|
|
12
|
|
4
|
13 import argparse
|
|
|
14 import utils.general_utils as utils
|
|
489
|
15 from typing import List
|
|
4
|
16 import os
|
|
489
|
17 import pandas as pd
|
|
4
|
18 import numpy as np
|
|
|
19 import cobra
|
|
|
20 import utils.CBS_backend as CBS_backend
|
|
|
21 from joblib import Parallel, delayed, cpu_count
|
|
|
22 from cobra.sampling import OptGPSampler
|
|
|
23 import sys
|
|
489
|
24 import utils.model_utils as model_utils
|
|
|
25
|
|
4
|
26
|
|
|
27 ################################# process args ###############################
|
|
489
|
28 def process_args(args: List[str] = None) -> argparse.Namespace:
|
|
4
|
29 """
|
|
|
30 Processes command-line arguments.
|
|
489
|
31
|
|
4
|
32 Args:
|
|
|
33 args (list): List of command-line arguments.
|
|
489
|
34
|
|
4
|
35 Returns:
|
|
|
36 Namespace: An object containing parsed arguments.
|
|
|
37 """
|
|
489
|
38 parser = argparse.ArgumentParser(usage='%(prog)s [options]',
|
|
|
39 description='process some value\'s')
|
|
|
40
|
|
|
41 parser.add_argument("-mo", "--model_upload", type=str,
|
|
|
42 help="path to input file with custom rules, if provided")
|
|
|
43
|
|
|
44 parser.add_argument("-mab", "--model_and_bounds", type=str,
|
|
|
45 choices=['True', 'False'],
|
|
|
46 required=True,
|
|
|
47 help="upload mode: True for model+bounds, False for complete models")
|
|
|
48
|
|
|
49 parser.add_argument("-ens", "--sampling_enabled", type=str,
|
|
|
50 choices=['true', 'false'],
|
|
|
51 required=True,
|
|
|
52 help="enable sampling: 'true' for sampling, 'false' for no sampling")
|
|
|
53
|
|
|
54 parser.add_argument('-ol', '--out_log',
|
|
|
55 help="Output log")
|
|
4
|
56
|
|
|
57 parser.add_argument('-td', '--tool_dir',
|
|
489
|
58 type=str,
|
|
|
59 required=True,
|
|
|
60 help='your tool directory')
|
|
4
|
61
|
|
|
62 parser.add_argument('-in', '--input',
|
|
489
|
63 required=True,
|
|
4
|
64 type=str,
|
|
489
|
65 help='input bounds files or complete model files')
|
|
4
|
66
|
|
489
|
67 parser.add_argument('-ni', '--name',
|
|
|
68 required=True,
|
|
|
69 type=str,
|
|
|
70 help='cell names')
|
|
4
|
71
|
|
|
72 parser.add_argument('-a', '--algorithm',
|
|
489
|
73 type=str,
|
|
|
74 choices=['OPTGP', 'CBS'],
|
|
|
75 required=True,
|
|
|
76 help='choose sampling algorithm')
|
|
4
|
77
|
|
489
|
78 parser.add_argument('-th', '--thinning',
|
|
|
79 type=int,
|
|
|
80 default=100,
|
|
|
81 required=True,
|
|
|
82 help='choose thinning')
|
|
4
|
83
|
|
489
|
84 parser.add_argument('-ns', '--n_samples',
|
|
|
85 type=int,
|
|
|
86 required=True,
|
|
|
87 help='choose how many samples (set to 0 for optimization only)')
|
|
|
88
|
|
|
89 parser.add_argument('-sd', '--seed',
|
|
|
90 type=int,
|
|
|
91 required=True,
|
|
|
92 help='seed for random number generation')
|
|
4
|
93
|
|
489
|
94 parser.add_argument('-nb', '--n_batches',
|
|
|
95 type=int,
|
|
|
96 required=True,
|
|
|
97 help='choose how many batches')
|
|
4
|
98
|
|
489
|
99 parser.add_argument('-opt', '--perc_opt',
|
|
|
100 type=float,
|
|
|
101 default=0.9,
|
|
|
102 required=False,
|
|
|
103 help='choose the fraction of optimality for FVA (0-1)')
|
|
4
|
104
|
|
489
|
105 parser.add_argument('-ot', '--output_type',
|
|
|
106 type=str,
|
|
|
107 required=True,
|
|
|
108 help='output type for sampling results')
|
|
4
|
109
|
|
489
|
110 parser.add_argument('-ota', '--output_type_analysis',
|
|
|
111 type=str,
|
|
|
112 required=False,
|
|
|
113 help='output type analysis (optimization methods)')
|
|
515
|
114
|
|
489
|
115 parser.add_argument('-idop', '--output_path',
|
|
|
116 type=str,
|
|
515
|
117 default='flux_simulation/',
|
|
|
118 help = 'output path for fluxes')
|
|
|
119
|
|
|
120 parser.add_argument('-otm', '--out_mean',
|
|
|
121 type = str,
|
|
517
|
122 required=False,
|
|
515
|
123 help = 'output of mean of fluxes')
|
|
147
|
124
|
|
515
|
125 parser.add_argument('-otmd', '--out_median',
|
|
|
126 type = str,
|
|
517
|
127 required=False,
|
|
515
|
128 help = 'output of median of fluxes')
|
|
|
129
|
|
|
130 parser.add_argument('-otq', '--out_quantiles',
|
|
|
131 type = str,
|
|
517
|
132 required=False,
|
|
515
|
133 help = 'output of quantiles of fluxes')
|
|
|
134
|
|
|
135 parser.add_argument('-otf', '--out_fluxes',
|
|
|
136 type = str,
|
|
517
|
137 required=False,
|
|
515
|
138 help = 'output of fluxes')
|
|
|
139 parser.add_argument('-otfva', '--out_fva',
|
|
|
140 type = str,
|
|
517
|
141 required=False,
|
|
515
|
142 help = 'output of FVA results')
|
|
|
143 parser.add_argument('-otp', '--out_pfba',
|
|
|
144 type = str,
|
|
517
|
145 required=False,
|
|
515
|
146 help = 'output of pFBA results')
|
|
|
147 parser.add_argument('-ots', '--out_sensitivity',
|
|
|
148 type = str,
|
|
517
|
149 required=False,
|
|
515
|
150 help = 'output of sensitivity results')
|
|
147
|
151 ARGS = parser.parse_args(args)
|
|
4
|
152 return ARGS
|
|
|
153 ########################### warning ###########################################
|
|
|
154 def warning(s :str) -> None:
|
|
|
155 """
|
|
|
156 Log a warning message to an output log file and print it to the console.
|
|
|
157
|
|
|
158 Args:
|
|
|
159 s (str): The warning message to be logged and printed.
|
|
|
160
|
|
|
161 Returns:
|
|
|
162 None
|
|
|
163 """
|
|
|
164 with open(ARGS.out_log, 'a') as log:
|
|
|
165 log.write(s + "\n\n")
|
|
|
166 print(s)
|
|
|
167
|
|
|
168
|
|
519
|
169 def write_to_file(dataset: pd.DataFrame, path: str, keep_index:bool=False, name:str=None)->None:
|
|
489
|
170 """
|
|
515
|
171 Write a DataFrame to a TSV file under path with a given base name.
|
|
489
|
172
|
|
|
173 Args:
|
|
|
174 dataset: The DataFrame to write.
|
|
519
|
175 name: Base file name (without extension). If None, 'path' is treated as the full file path.
|
|
515
|
176 path: Directory path where the file will be saved.
|
|
489
|
177 keep_index: Whether to keep the DataFrame index in the file.
|
|
|
178
|
|
|
179 Returns:
|
|
|
180 None
|
|
|
181 """
|
|
4
|
182 dataset.index.name = 'Reactions'
|
|
519
|
183 if name:
|
|
|
184 dataset.to_csv(os.path.join(path, name + ".csv"), sep = '\t', index = keep_index)
|
|
|
185 else:
|
|
|
186 dataset.to_csv(path, sep = '\t', index = keep_index)
|
|
4
|
187
|
|
|
188 ############################ dataset input ####################################
|
|
|
189 def read_dataset(data :str, name :str) -> pd.DataFrame:
|
|
|
190 """
|
|
|
191 Read a dataset from a CSV file and return it as a pandas DataFrame.
|
|
|
192
|
|
|
193 Args:
|
|
|
194 data (str): Path to the CSV file containing the dataset.
|
|
|
195 name (str): Name of the dataset, used in error messages.
|
|
|
196
|
|
|
197 Returns:
|
|
|
198 pandas.DataFrame: DataFrame containing the dataset.
|
|
|
199
|
|
|
200 Raises:
|
|
|
201 pd.errors.EmptyDataError: If the CSV file is empty.
|
|
|
202 sys.exit: If the CSV file has the wrong format, the execution is aborted.
|
|
|
203 """
|
|
|
204 try:
|
|
|
205 dataset = pd.read_csv(data, sep = '\t', header = 0, index_col=0, engine='python')
|
|
|
206 except pd.errors.EmptyDataError:
|
|
|
207 sys.exit('Execution aborted: wrong format of ' + name + '\n')
|
|
|
208 if len(dataset.columns) < 2:
|
|
|
209 sys.exit('Execution aborted: wrong format of ' + name + '\n')
|
|
|
210 return dataset
|
|
|
211
|
|
|
212
|
|
|
213
|
|
489
|
214 def OPTGP_sampler(model: cobra.Model, model_name: str, n_samples: int = 1000, thinning: int = 100, n_batches: int = 1, seed: int = 0) -> None:
|
|
4
|
215 """
|
|
|
216 Samples from the OPTGP (Optimal Global Perturbation) algorithm and saves the results to CSV files.
|
|
489
|
217
|
|
4
|
218 Args:
|
|
|
219 model (cobra.Model): The COBRA model to sample from.
|
|
|
220 model_name (str): The name of the model, used in naming output files.
|
|
|
221 n_samples (int, optional): Number of samples per batch. Default is 1000.
|
|
|
222 thinning (int, optional): Thinning parameter for the sampler. Default is 100.
|
|
|
223 n_batches (int, optional): Number of batches to run. Default is 1.
|
|
|
224 seed (int, optional): Random seed for reproducibility. Default is 0.
|
|
489
|
225
|
|
4
|
226 Returns:
|
|
|
227 None
|
|
|
228 """
|
|
489
|
229 import numpy as np
|
|
|
230
|
|
|
231 # Get reaction IDs for consistent column ordering
|
|
|
232 reaction_ids = [rxn.id for rxn in model.reactions]
|
|
|
233
|
|
|
234 # Sample and save each batch as numpy file
|
|
|
235 for i in range(n_batches):
|
|
4
|
236 optgp = OptGPSampler(model, thinning, seed)
|
|
|
237 samples = optgp.sample(n_samples)
|
|
489
|
238
|
|
|
239 # Save as numpy array (more memory efficient)
|
|
|
240 batch_filename = f"{ARGS.output_path}/{model_name}_{i}_OPTGP.npy"
|
|
|
241 np.save(batch_filename, samples.to_numpy())
|
|
|
242
|
|
|
243 seed += 1
|
|
|
244
|
|
|
245 # Merge all batches into a single DataFrame
|
|
|
246 all_samples = []
|
|
|
247
|
|
|
248 for i in range(n_batches):
|
|
|
249 batch_filename = f"{ARGS.output_path}/{model_name}_{i}_OPTGP.npy"
|
|
|
250 batch_data = np.load(batch_filename, allow_pickle=True)
|
|
|
251 all_samples.append(batch_data)
|
|
|
252
|
|
|
253 # Concatenate all batches
|
|
|
254 samplesTotal_array = np.vstack(all_samples)
|
|
|
255
|
|
|
256 # Convert back to DataFrame with proper column names
|
|
|
257 samplesTotal = pd.DataFrame(samplesTotal_array, columns=reaction_ids)
|
|
|
258
|
|
|
259 # Save the final merged result as CSV
|
|
519
|
260 write_to_file(samplesTotal.T, ARGS.output_path, True, name=model_name)
|
|
489
|
261
|
|
|
262 # Clean up temporary numpy files
|
|
|
263 for i in range(n_batches):
|
|
|
264 batch_filename = f"{ARGS.output_path}/{model_name}_{i}_OPTGP.npy"
|
|
|
265 if os.path.exists(batch_filename):
|
|
|
266 os.remove(batch_filename)
|
|
4
|
267
|
|
|
268
|
|
489
|
269 def CBS_sampler(model: cobra.Model, model_name: str, n_samples: int = 1000, n_batches: int = 1, seed: int = 0) -> None:
|
|
4
|
270 """
|
|
|
271 Samples using the CBS (Constraint-based Sampling) algorithm and saves the results to CSV files.
|
|
489
|
272
|
|
4
|
273 Args:
|
|
|
274 model (cobra.Model): The COBRA model to sample from.
|
|
|
275 model_name (str): The name of the model, used in naming output files.
|
|
|
276 n_samples (int, optional): Number of samples per batch. Default is 1000.
|
|
|
277 n_batches (int, optional): Number of batches to run. Default is 1.
|
|
|
278 seed (int, optional): Random seed for reproducibility. Default is 0.
|
|
489
|
279
|
|
4
|
280 Returns:
|
|
|
281 None
|
|
|
282 """
|
|
489
|
283 import numpy as np
|
|
|
284
|
|
|
285 # Get reaction IDs for consistent column ordering
|
|
|
286 reaction_ids = [reaction.id for reaction in model.reactions]
|
|
|
287
|
|
|
288 # Perform FVA analysis once for all batches
|
|
|
289 df_FVA = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0).round(6)
|
|
|
290
|
|
|
291 # Generate random objective functions for all samples across all batches
|
|
|
292 df_coefficients = CBS_backend.randomObjectiveFunction(model, n_samples * n_batches, df_FVA, seed=seed)
|
|
4
|
293
|
|
489
|
294 # Sample and save each batch as numpy file
|
|
|
295 for i in range(n_batches):
|
|
|
296 samples = pd.DataFrame(columns=reaction_ids, index=range(n_samples))
|
|
|
297
|
|
4
|
298 try:
|
|
489
|
299 CBS_backend.randomObjectiveFunctionSampling(
|
|
|
300 model,
|
|
|
301 n_samples,
|
|
|
302 df_coefficients.iloc[:, i * n_samples:(i + 1) * n_samples],
|
|
|
303 samples
|
|
|
304 )
|
|
4
|
305 except Exception as e:
|
|
|
306 utils.logWarning(
|
|
489
|
307 f"Warning: GLPK solver has failed for {model_name}. Trying with COBRA interface. Error: {str(e)}",
|
|
|
308 ARGS.out_log
|
|
|
309 )
|
|
|
310 CBS_backend.randomObjectiveFunctionSampling_cobrapy(
|
|
|
311 model,
|
|
|
312 n_samples,
|
|
|
313 df_coefficients.iloc[:, i * n_samples:(i + 1) * n_samples],
|
|
|
314 samples
|
|
|
315 )
|
|
|
316
|
|
|
317 # Save as numpy array (more memory efficient)
|
|
|
318 batch_filename = f"{ARGS.output_path}/{model_name}_{i}_CBS.npy"
|
|
|
319 utils.logWarning(batch_filename, ARGS.out_log)
|
|
|
320 np.save(batch_filename, samples.to_numpy())
|
|
|
321
|
|
|
322 # Merge all batches into a single DataFrame
|
|
|
323 all_samples = []
|
|
|
324
|
|
|
325 for i in range(n_batches):
|
|
|
326 batch_filename = f"{ARGS.output_path}/{model_name}_{i}_CBS.npy"
|
|
|
327 batch_data = np.load(batch_filename, allow_pickle=True)
|
|
|
328 all_samples.append(batch_data)
|
|
|
329
|
|
|
330 # Concatenate all batches
|
|
|
331 samplesTotal_array = np.vstack(all_samples)
|
|
|
332
|
|
|
333 # Convert back to DataFrame with proper column namesq
|
|
|
334 samplesTotal = pd.DataFrame(samplesTotal_array, columns=reaction_ids)
|
|
|
335
|
|
|
336 # Save the final merged result as CSV
|
|
519
|
337 write_to_file(samplesTotal.T, ARGS.output_path, True, name=model_name)
|
|
489
|
338
|
|
|
339 # Clean up temporary numpy files
|
|
|
340 for i in range(n_batches):
|
|
|
341 batch_filename = f"{ARGS.output_path}/{model_name}_{i}_CBS.npy"
|
|
|
342 if os.path.exists(batch_filename):
|
|
|
343 os.remove(batch_filename)
|
|
4
|
344
|
|
|
345
|
|
489
|
346
|
|
|
347 def model_sampler_with_bounds(model_input_original: cobra.Model, bounds_path: str, cell_name: str) -> List[pd.DataFrame]:
|
|
4
|
348 """
|
|
489
|
349 MODE 1: Prepares the model with bounds from separate bounds file and performs sampling.
|
|
4
|
350
|
|
|
351 Args:
|
|
|
352 model_input_original (cobra.Model): The original COBRA model.
|
|
|
353 bounds_path (str): Path to the CSV file containing the bounds dataset.
|
|
|
354 cell_name (str): Name of the cell, used to generate filenames for output.
|
|
|
355
|
|
|
356 Returns:
|
|
|
357 List[pd.DataFrame]: A list of DataFrames containing statistics and analysis results.
|
|
|
358 """
|
|
|
359
|
|
|
360 model_input = model_input_original.copy()
|
|
|
361 bounds_df = read_dataset(bounds_path, "bounds dataset")
|
|
|
362
|
|
489
|
363 # Apply bounds to model
|
|
|
364 for rxn_index, row in bounds_df.iterrows():
|
|
|
365 try:
|
|
|
366 model_input.reactions.get_by_id(rxn_index).lower_bound = row.lower_bound
|
|
|
367 model_input.reactions.get_by_id(rxn_index).upper_bound = row.upper_bound
|
|
|
368 except KeyError:
|
|
|
369 warning(f"Warning: Reaction {rxn_index} not found in model. Skipping.")
|
|
4
|
370
|
|
489
|
371 return perform_sampling_and_analysis(model_input, cell_name)
|
|
|
372
|
|
4
|
373
|
|
489
|
374 def perform_sampling_and_analysis(model_input: cobra.Model, cell_name: str) -> List[pd.DataFrame]:
|
|
|
375 """
|
|
|
376 Common function to perform sampling and analysis on a prepared model.
|
|
4
|
377
|
|
489
|
378 Args:
|
|
|
379 model_input (cobra.Model): The prepared COBRA model with bounds applied.
|
|
|
380 cell_name (str): Name of the cell, used to generate filenames for output.
|
|
4
|
381
|
|
489
|
382 Returns:
|
|
|
383 List[pd.DataFrame]: A list of DataFrames containing statistics and analysis results.
|
|
|
384 """
|
|
4
|
385
|
|
|
386 returnList = []
|
|
489
|
387
|
|
|
388 if ARGS.sampling_enabled == "true":
|
|
|
389
|
|
|
390 if ARGS.algorithm == 'OPTGP':
|
|
|
391 OPTGP_sampler(model_input, cell_name, ARGS.n_samples, ARGS.thinning, ARGS.n_batches, ARGS.seed)
|
|
|
392 elif ARGS.algorithm == 'CBS':
|
|
|
393 CBS_sampler(model_input, cell_name, ARGS.n_samples, ARGS.n_batches, ARGS.seed)
|
|
|
394
|
|
|
395 df_mean, df_median, df_quantiles = fluxes_statistics(cell_name, ARGS.output_types)
|
|
|
396
|
|
|
397 if("fluxes" not in ARGS.output_types):
|
|
|
398 os.remove(ARGS.output_path + "/" + cell_name + '.csv')
|
|
|
399
|
|
|
400 returnList = [df_mean, df_median, df_quantiles]
|
|
4
|
401
|
|
210
|
402 df_pFBA, df_FVA, df_sensitivity = fluxes_analysis(model_input, cell_name, ARGS.output_type_analysis)
|
|
4
|
403
|
|
|
404 if("pFBA" in ARGS.output_type_analysis):
|
|
|
405 returnList.append(df_pFBA)
|
|
|
406 if("FVA" in ARGS.output_type_analysis):
|
|
|
407 returnList.append(df_FVA)
|
|
|
408 if("sensitivity" in ARGS.output_type_analysis):
|
|
|
409 returnList.append(df_sensitivity)
|
|
|
410
|
|
|
411 return returnList
|
|
|
412
|
|
|
413 def fluxes_statistics(model_name: str, output_types:List)-> List[pd.DataFrame]:
|
|
|
414 """
|
|
|
415 Computes statistics (mean, median, quantiles) for the fluxes.
|
|
|
416
|
|
|
417 Args:
|
|
|
418 model_name (str): Name of the model, used in filename for input.
|
|
|
419 output_types (List[str]): Types of statistics to compute (mean, median, quantiles).
|
|
|
420
|
|
|
421 Returns:
|
|
|
422 List[pd.DataFrame]: List of DataFrames containing mean, median, and quantiles statistics.
|
|
|
423 """
|
|
|
424
|
|
|
425 df_mean = pd.DataFrame()
|
|
|
426 df_median= pd.DataFrame()
|
|
|
427 df_quantiles= pd.DataFrame()
|
|
|
428
|
|
161
|
429 df_samples = pd.read_csv(ARGS.output_path + "/" + model_name + '.csv', sep = '\t', index_col = 0).T
|
|
4
|
430 df_samples = df_samples.round(8)
|
|
|
431
|
|
|
432 for output_type in output_types:
|
|
|
433 if(output_type == "mean"):
|
|
|
434 df_mean = df_samples.mean()
|
|
|
435 df_mean = df_mean.to_frame().T
|
|
|
436 df_mean = df_mean.reset_index(drop=True)
|
|
|
437 df_mean.index = [model_name]
|
|
|
438 elif(output_type == "median"):
|
|
|
439 df_median = df_samples.median()
|
|
|
440 df_median = df_median.to_frame().T
|
|
|
441 df_median = df_median.reset_index(drop=True)
|
|
|
442 df_median.index = [model_name]
|
|
|
443 elif(output_type == "quantiles"):
|
|
|
444 newRow = []
|
|
|
445 cols = []
|
|
|
446 for rxn in df_samples.columns:
|
|
|
447 quantiles = df_samples[rxn].quantile([0.25, 0.50, 0.75])
|
|
|
448 newRow.append(quantiles[0.25])
|
|
|
449 cols.append(rxn + "_q1")
|
|
|
450 newRow.append(quantiles[0.5])
|
|
|
451 cols.append(rxn + "_q2")
|
|
|
452 newRow.append(quantiles[0.75])
|
|
|
453 cols.append(rxn + "_q3")
|
|
|
454 df_quantiles = pd.DataFrame(columns=cols)
|
|
|
455 df_quantiles.loc[0] = newRow
|
|
|
456 df_quantiles = df_quantiles.reset_index(drop=True)
|
|
|
457 df_quantiles.index = [model_name]
|
|
|
458
|
|
|
459 return df_mean, df_median, df_quantiles
|
|
|
460
|
|
|
461 def fluxes_analysis(model:cobra.Model, model_name:str, output_types:List)-> List[pd.DataFrame]:
|
|
|
462 """
|
|
489
|
463 Performs flux analysis including pFBA, FVA, and sensitivity analysis. The objective function
|
|
|
464 is assumed to be already set in the model.
|
|
4
|
465
|
|
|
466 Args:
|
|
|
467 model (cobra.Model): The COBRA model to analyze.
|
|
|
468 model_name (str): Name of the model, used in filenames for output.
|
|
|
469 output_types (List[str]): Types of analysis to perform (pFBA, FVA, sensitivity).
|
|
|
470
|
|
|
471 Returns:
|
|
|
472 List[pd.DataFrame]: List of DataFrames containing pFBA, FVA, and sensitivity analysis results.
|
|
|
473 """
|
|
|
474
|
|
|
475 df_pFBA = pd.DataFrame()
|
|
|
476 df_FVA= pd.DataFrame()
|
|
|
477 df_sensitivity= pd.DataFrame()
|
|
|
478
|
|
|
479 for output_type in output_types:
|
|
|
480 if(output_type == "pFBA"):
|
|
|
481 solution = cobra.flux_analysis.pfba(model)
|
|
|
482 fluxes = solution.fluxes
|
|
489
|
483 df_pFBA.loc[0,[rxn.id for rxn in model.reactions]] = fluxes.tolist()
|
|
4
|
484 df_pFBA = df_pFBA.reset_index(drop=True)
|
|
|
485 df_pFBA.index = [model_name]
|
|
|
486 df_pFBA = df_pFBA.astype(float).round(6)
|
|
|
487 elif(output_type == "FVA"):
|
|
489
|
488 fva = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=ARGS.perc_opt, processes=1).round(8)
|
|
4
|
489 columns = []
|
|
|
490 for rxn in fva.index.to_list():
|
|
|
491 columns.append(rxn + "_min")
|
|
|
492 columns.append(rxn + "_max")
|
|
|
493 df_FVA= pd.DataFrame(columns = columns)
|
|
|
494 for index_rxn, row in fva.iterrows():
|
|
|
495 df_FVA.loc[0, index_rxn+ "_min"] = fva.loc[index_rxn, "minimum"]
|
|
|
496 df_FVA.loc[0, index_rxn+ "_max"] = fva.loc[index_rxn, "maximum"]
|
|
|
497 df_FVA = df_FVA.reset_index(drop=True)
|
|
|
498 df_FVA.index = [model_name]
|
|
|
499 df_FVA = df_FVA.astype(float).round(6)
|
|
|
500 elif(output_type == "sensitivity"):
|
|
|
501 solution_original = model.optimize().objective_value
|
|
|
502 reactions = model.reactions
|
|
|
503 single = cobra.flux_analysis.single_reaction_deletion(model)
|
|
|
504 newRow = []
|
|
|
505 df_sensitivity = pd.DataFrame(columns = [rxn.id for rxn in reactions], index = [model_name])
|
|
|
506 for rxn in reactions:
|
|
|
507 newRow.append(single.knockout[rxn.id].growth.values[0]/solution_original)
|
|
|
508 df_sensitivity.loc[model_name] = newRow
|
|
|
509 df_sensitivity = df_sensitivity.astype(float).round(6)
|
|
|
510 return df_pFBA, df_FVA, df_sensitivity
|
|
|
511
|
|
|
512 ############################# main ###########################################
|
|
489
|
513 def main(args: List[str] = None) -> None:
|
|
4
|
514 """
|
|
489
|
515 Initialize and run sampling/analysis based on the frontend input arguments.
|
|
4
|
516
|
|
|
517 Returns:
|
|
|
518 None
|
|
|
519 """
|
|
|
520
|
|
489
|
521 num_processors = max(1, cpu_count() - 1)
|
|
4
|
522
|
|
|
523 global ARGS
|
|
147
|
524 ARGS = process_args(args)
|
|
158
|
525
|
|
515
|
526 if not os.path.exists('flux_simulation'):
|
|
|
527 os.makedirs('flux_simulation')
|
|
489
|
528
|
|
|
529 # --- Normalize inputs (the tool may pass comma-separated --input and either --name or --names) ---
|
|
|
530 ARGS.input_files = ARGS.input.split(",") if ARGS.input else []
|
|
|
531 ARGS.file_names = ARGS.name.split(",")
|
|
|
532 # output types (required) -> list
|
|
|
533 ARGS.output_types = ARGS.output_type.split(",") if ARGS.output_type else []
|
|
|
534 # optional analysis output types -> list or empty
|
|
|
535 ARGS.output_type_analysis = ARGS.output_type_analysis.split(",") if ARGS.output_type_analysis else []
|
|
|
536
|
|
|
537 # Determine if sampling should be performed
|
|
|
538 if ARGS.sampling_enabled == "true":
|
|
|
539 perform_sampling = True
|
|
4
|
540 else:
|
|
489
|
541 perform_sampling = False
|
|
334
|
542
|
|
489
|
543 print("=== INPUT FILES ===")
|
|
|
544 print(f"{ARGS.input_files}")
|
|
|
545 print(f"{ARGS.file_names}")
|
|
|
546 print(f"{ARGS.output_type}")
|
|
|
547 print(f"{ARGS.output_types}")
|
|
|
548 print(f"{ARGS.output_type_analysis}")
|
|
|
549 print(f"Sampling enabled: {perform_sampling} (n_samples: {ARGS.n_samples})")
|
|
4
|
550
|
|
489
|
551 if ARGS.model_and_bounds == "True":
|
|
|
552 # MODE 1: Model + bounds (separate files)
|
|
|
553 print("=== MODE 1: Model + Bounds (separate files) ===")
|
|
|
554
|
|
|
555 # Load base model
|
|
|
556 if not ARGS.model_upload:
|
|
|
557 sys.exit("Error: model_upload is required for Mode 1")
|
|
4
|
558
|
|
489
|
559 base_model = model_utils.build_cobra_model_from_csv(ARGS.model_upload)
|
|
4
|
560
|
|
489
|
561 validation = model_utils.validate_model(base_model)
|
|
|
562
|
|
|
563 print("\n=== MODEL VALIDATION ===")
|
|
|
564 for key, value in validation.items():
|
|
|
565 print(f"{key}: {value}")
|
|
|
566
|
|
|
567 # Set solver verbosity to 1 to see warning and error messages only.
|
|
|
568 base_model.solver.configuration.verbosity = 1
|
|
4
|
569
|
|
489
|
570 # Process each bounds file with the base model
|
|
|
571 results = Parallel(n_jobs=num_processors)(
|
|
|
572 delayed(model_sampler_with_bounds)(base_model, bounds_file, cell_name)
|
|
|
573 for bounds_file, cell_name in zip(ARGS.input_files, ARGS.file_names)
|
|
|
574 )
|
|
4
|
575
|
|
489
|
576 else:
|
|
|
577 # MODE 2: Multiple complete models
|
|
|
578 print("=== MODE 2: Multiple complete models ===")
|
|
|
579
|
|
|
580 # Process each complete model file
|
|
|
581 results = Parallel(n_jobs=num_processors)(
|
|
|
582 delayed(perform_sampling_and_analysis)(model_utils.build_cobra_model_from_csv(model_file), cell_name)
|
|
|
583 for model_file, cell_name in zip(ARGS.input_files, ARGS.file_names)
|
|
|
584 )
|
|
|
585
|
|
|
586 # Handle sampling outputs (only if sampling was performed)
|
|
|
587 if perform_sampling:
|
|
|
588 print("=== PROCESSING SAMPLING RESULTS ===")
|
|
|
589
|
|
|
590 all_mean = pd.concat([result[0] for result in results], ignore_index=False)
|
|
|
591 all_median = pd.concat([result[1] for result in results], ignore_index=False)
|
|
|
592 all_quantiles = pd.concat([result[2] for result in results], ignore_index=False)
|
|
4
|
593
|
|
489
|
594 if "mean" in ARGS.output_types:
|
|
|
595 all_mean = all_mean.fillna(0.0)
|
|
|
596 all_mean = all_mean.sort_index()
|
|
519
|
597 write_to_file(all_mean.T, ARGS.out_mean, True)
|
|
489
|
598
|
|
|
599 if "median" in ARGS.output_types:
|
|
|
600 all_median = all_median.fillna(0.0)
|
|
|
601 all_median = all_median.sort_index()
|
|
519
|
602 write_to_file(all_median.T, ARGS.out_median, True)
|
|
489
|
603
|
|
|
604 if "quantiles" in ARGS.output_types:
|
|
|
605 all_quantiles = all_quantiles.fillna(0.0)
|
|
|
606 all_quantiles = all_quantiles.sort_index()
|
|
519
|
607 write_to_file(all_quantiles.T, ARGS.out_quantiles, True)
|
|
489
|
608 else:
|
|
|
609 print("=== SAMPLING SKIPPED (n_samples = 0 or sampling disabled) ===")
|
|
|
610
|
|
|
611 # Handle optimization analysis outputs (always available)
|
|
|
612 print("=== PROCESSING OPTIMIZATION RESULTS ===")
|
|
4
|
613
|
|
489
|
614 # Determine the starting index for optimization results
|
|
|
615 # If sampling was performed, optimization results start at index 3
|
|
|
616 # If no sampling, optimization results start at index 0
|
|
|
617 index_result = 3 if perform_sampling else 0
|
|
|
618
|
|
|
619 if "pFBA" in ARGS.output_type_analysis:
|
|
4
|
620 all_pFBA = pd.concat([result[index_result] for result in results], ignore_index=False)
|
|
|
621 all_pFBA = all_pFBA.sort_index()
|
|
519
|
622 write_to_file(all_pFBA.T, ARGS.out_pfba, True)
|
|
489
|
623 index_result += 1
|
|
|
624
|
|
|
625 if "FVA" in ARGS.output_type_analysis:
|
|
|
626 all_FVA = pd.concat([result[index_result] for result in results], ignore_index=False)
|
|
4
|
627 all_FVA = all_FVA.sort_index()
|
|
519
|
628 write_to_file(all_FVA.T, ARGS.out_fva, True)
|
|
489
|
629 index_result += 1
|
|
|
630
|
|
|
631 if "sensitivity" in ARGS.output_type_analysis:
|
|
4
|
632 all_sensitivity = pd.concat([result[index_result] for result in results], ignore_index=False)
|
|
|
633 all_sensitivity = all_sensitivity.sort_index()
|
|
519
|
634 write_to_file(all_sensitivity.T, ARGS.out_sensitivity, True)
|
|
4
|
635
|
|
489
|
636 return
|
|
4
|
637
|
|
|
638 ##############################################################################
|
|
|
639 if __name__ == "__main__":
|
|
|
640 main() |
