|
542
|
1 from __future__ import division
|
|
|
2 import csv
|
|
|
3 from enum import Enum
|
|
|
4 import re
|
|
|
5 import sys
|
|
|
6 import numpy as np
|
|
|
7 import pandas as pd
|
|
|
8 import itertools as it
|
|
|
9 import scipy.stats as st
|
|
|
10 import lxml.etree as ET
|
|
|
11 import math
|
|
|
12 try:
|
|
|
13 from .utils import general_utils as utils
|
|
|
14 except:
|
|
|
15 import utils.general_utils as utils
|
|
|
16 from PIL import Image
|
|
|
17 import os
|
|
|
18 import copy
|
|
|
19 import argparse
|
|
|
20 import pyvips
|
|
|
21 from PIL import Image
|
|
|
22 from typing import Tuple, Union, Optional, List, Dict
|
|
|
23 import matplotlib.pyplot as plt
|
|
|
24
|
|
|
25 ERRORS = []
|
|
|
26 ########################## argparse ##########################################
|
|
|
27 ARGS :argparse.Namespace
|
|
|
28 def process_args(args:List[str] = None) -> argparse.Namespace:
|
|
|
29 """
|
|
|
30 Interfaces the script of a module with its frontend, making the user's choices for various parameters available as values in code.
|
|
|
31
|
|
|
32 Args:
|
|
|
33 args : Always obtained (in file) from sys.argv
|
|
|
34
|
|
|
35 Returns:
|
|
|
36 Namespace : An object containing the parsed arguments
|
|
|
37 """
|
|
|
38 parser = argparse.ArgumentParser(
|
|
|
39 usage = "%(prog)s [options]",
|
|
|
40 description = "process some value's genes to create a comparison's map.")
|
|
|
41
|
|
|
42 #General:
|
|
|
43 parser.add_argument(
|
|
|
44 '-td', '--tool_dir',
|
|
|
45 type = str,
|
|
|
46 default = os.path.dirname(os.path.abspath(__file__)),
|
|
|
47 help = 'your tool directory (default: auto-detected package location)')
|
|
|
48
|
|
|
49 parser.add_argument('-on', '--control', type = str)
|
|
|
50 parser.add_argument('-ol', '--out_log', help = "Output log")
|
|
|
51
|
|
|
52 #Computation details:
|
|
|
53 parser.add_argument(
|
|
|
54 '-co', '--comparison',
|
|
|
55 type = str,
|
|
|
56 default = 'manyvsmany',
|
|
|
57 choices = ['manyvsmany', 'onevsrest', 'onevsmany'])
|
|
|
58
|
|
|
59 parser.add_argument(
|
|
|
60 '-te' ,'--test',
|
|
|
61 type = str,
|
|
|
62 default = 'ks',
|
|
|
63 choices = ['ks', 'ttest_p', 'ttest_ind', 'wilcoxon', 'mw'],
|
|
|
64 help = 'Statistical test to use (default: %(default)s)')
|
|
|
65
|
|
|
66 parser.add_argument(
|
|
|
67 '-pv' ,'--pValue',
|
|
|
68 type = float,
|
|
|
69 default = 0.1,
|
|
|
70 help = 'P-Value threshold (default: %(default)s)')
|
|
|
71
|
|
|
72 parser.add_argument(
|
|
|
73 '-adj' ,'--adjusted',
|
|
|
74 type = utils.Bool("adjusted"), default = False,
|
|
|
75 help = 'Apply the FDR (Benjamini-Hochberg) correction (default: %(default)s)')
|
|
|
76
|
|
|
77 parser.add_argument(
|
|
|
78 '-fc', '--fChange',
|
|
|
79 type = float,
|
|
|
80 default = 1.5,
|
|
|
81 help = 'Fold-Change threshold (default: %(default)s)')
|
|
|
82
|
|
|
83 parser.add_argument(
|
|
|
84 '-op', '--option',
|
|
|
85 type = str,
|
|
|
86 choices = ['datasets', 'dataset_class'],
|
|
|
87 help='dataset or dataset and class')
|
|
|
88
|
|
|
89 parser.add_argument(
|
|
|
90 '-idf', '--input_data_fluxes',
|
|
|
91 type = str,
|
|
|
92 help = 'input dataset fluxes')
|
|
|
93
|
|
|
94 parser.add_argument(
|
|
|
95 '-icf', '--input_class_fluxes',
|
|
|
96 type = str,
|
|
|
97 help = 'sample group specification fluxes')
|
|
|
98
|
|
|
99 parser.add_argument(
|
|
|
100 '-idsf', '--input_datas_fluxes',
|
|
|
101 type = str,
|
|
|
102 nargs = '+',
|
|
|
103 help = 'input datasets fluxes')
|
|
|
104
|
|
|
105 parser.add_argument(
|
|
|
106 '-naf', '--names_fluxes',
|
|
|
107 type = str,
|
|
|
108 nargs = '+',
|
|
|
109 help = 'input names fluxes')
|
|
|
110
|
|
|
111 #Output:
|
|
|
112 parser.add_argument(
|
|
|
113 "-gs", "--generate_svg",
|
|
|
114 type = utils.Bool("generate_svg"), default = True,
|
|
|
115 help = "choose whether to generate svg")
|
|
|
116
|
|
|
117 parser.add_argument(
|
|
|
118 "-gp", "--generate_pdf",
|
|
|
119 type = utils.Bool("generate_pdf"), default = True,
|
|
|
120 help = "choose whether to generate pdf")
|
|
|
121
|
|
|
122 parser.add_argument(
|
|
|
123 '-cm', '--custom_map',
|
|
|
124 type = str,
|
|
|
125 help='custom map to use')
|
|
|
126
|
|
|
127 parser.add_argument(
|
|
|
128 '-mc', '--choice_map',
|
|
|
129 type = utils.Model, default = utils.Model.HMRcore,
|
|
|
130 choices = [utils.Model.HMRcore, utils.Model.ENGRO2, utils.Model.Custom])
|
|
|
131
|
|
|
132 parser.add_argument(
|
|
|
133 '-colorm', '--color_map',
|
|
|
134 type = str,
|
|
|
135 choices = ["jet", "viridis"])
|
|
|
136
|
|
|
137 parser.add_argument(
|
|
|
138 '-idop', '--output_path',
|
|
|
139 type = str,
|
|
|
140 default='result',
|
|
|
141 help = 'output path for maps')
|
|
|
142
|
|
|
143 args :argparse.Namespace = parser.parse_args(args)
|
|
|
144 args.net = True # TODO SICCOME I FLUSSI POSSONO ESSERE ANCHE NEGATIVI SONO SEMPRE CONSIDERATI NETTI
|
|
|
145
|
|
|
146 return args
|
|
|
147
|
|
|
148 ############################ dataset input ####################################
|
|
|
149 def read_dataset(data :str, name :str) -> pd.DataFrame:
|
|
|
150 """
|
|
|
151 Tries to read the dataset from its path (data) as a tsv and turns it into a DataFrame.
|
|
|
152
|
|
|
153 Args:
|
|
|
154 data : filepath of a dataset (from frontend input params or literals upon calling)
|
|
|
155 name : name associated with the dataset (from frontend input params or literals upon calling)
|
|
|
156
|
|
|
157 Returns:
|
|
|
158 pd.DataFrame : dataset in a runtime operable shape
|
|
|
159
|
|
|
160 Raises:
|
|
|
161 sys.exit : if there's no data (pd.errors.EmptyDataError) or if the dataset has less than 2 columns
|
|
|
162 """
|
|
|
163 try:
|
|
|
164 dataset = pd.read_csv(data, sep = '\t', header = 0, engine='python')
|
|
|
165 except pd.errors.EmptyDataError:
|
|
|
166 sys.exit('Execution aborted: wrong format of ' + name + '\n')
|
|
|
167 if len(dataset.columns) < 2:
|
|
|
168 sys.exit('Execution aborted: wrong format of ' + name + '\n')
|
|
|
169 return dataset
|
|
|
170
|
|
|
171 ############################ dataset name #####################################
|
|
|
172 def name_dataset(name_data :str, count :int) -> str:
|
|
|
173 """
|
|
|
174 Produces a unique name for a dataset based on what was provided by the user. The default name for any dataset is "Dataset", thus if the user didn't change it this function appends f"_{count}" to make it unique.
|
|
|
175
|
|
|
176 Args:
|
|
|
177 name_data : name associated with the dataset (from frontend input params)
|
|
|
178 count : counter from 1 to make these names unique (external)
|
|
|
179
|
|
|
180 Returns:
|
|
|
181 str : the name made unique
|
|
|
182 """
|
|
|
183 if str(name_data) == 'Dataset':
|
|
|
184 return str(name_data) + '_' + str(count)
|
|
|
185 else:
|
|
|
186 return str(name_data)
|
|
|
187
|
|
|
188 ############################ map_methods ######################################
|
|
|
189 FoldChange = Union[float, int, str] # Union[float, Literal[0, "-INF", "INF"]]
|
|
|
190 def fold_change(avg1 :float, avg2 :float) -> FoldChange:
|
|
|
191 """
|
|
|
192 Calculates the fold change between two gene expression values.
|
|
|
193
|
|
|
194 Args:
|
|
|
195 avg1 : average expression value from one dataset avg2 : average expression value from the other dataset
|
|
|
196
|
|
|
197 Returns:
|
|
|
198 FoldChange :
|
|
|
199 0 : when both input values are 0
|
|
|
200 "-INF" : when avg1 is 0
|
|
|
201 "INF" : when avg2 is 0
|
|
|
202 float : for any other combination of values
|
|
|
203 """
|
|
|
204 if avg1 == 0 and avg2 == 0:
|
|
|
205 return 0
|
|
|
206 elif avg1 == 0:
|
|
|
207 return '-INF'
|
|
|
208 elif avg2 == 0:
|
|
|
209 return 'INF'
|
|
|
210 else: # (threshold_F_C - 1) / (abs(threshold_F_C) + 1) con threshold_F_C > 1
|
|
|
211 return (avg1 - avg2) / (abs(avg1) + abs(avg2))
|
|
|
212
|
|
|
213 def getElementById(reactionId :str, metabMap :ET.ElementTree) -> utils.Result[ET.Element, utils.Result.ResultErr]:
|
|
|
214 """
|
|
|
215 Finds any element in the given map with the given ID. ID uniqueness in an svg file is recommended but
|
|
|
216 not enforced, if more than one element with the exact ID is found only the first will be returned.
|
|
|
217
|
|
|
218 Args:
|
|
|
219 reactionId (str): exact ID of the requested element.
|
|
|
220 metabMap (ET.ElementTree): metabolic map containing the element.
|
|
|
221
|
|
|
222 Returns:
|
|
|
223 utils.Result[ET.Element, ResultErr]: result of the search, either the first match found or a ResultErr.
|
|
|
224 """
|
|
|
225 return utils.Result.Ok(
|
|
|
226 f"//*[@id=\"{reactionId}\"]").map(
|
|
|
227 lambda xPath : metabMap.xpath(xPath)[0]).mapErr(
|
|
|
228 lambda _ : utils.Result.ResultErr(f"No elements with ID \"{reactionId}\" found in map"))
|
|
|
229 # ^^^ we shamelessly ignore the contents of the IndexError, it offers nothing to the user.
|
|
|
230
|
|
|
231 def styleMapElement(element :ET.Element, styleStr :str) -> None:
|
|
|
232 currentStyles :str = element.get("style", "")
|
|
|
233 if re.search(r";stroke:[^;]+;stroke-width:[^;]+;stroke-dasharray:[^;]+$", currentStyles):
|
|
|
234 currentStyles = ';'.join(currentStyles.split(';')[:-3])
|
|
|
235
|
|
|
236 element.set("style", currentStyles + styleStr)
|
|
|
237
|
|
|
238 class ReactionDirection(Enum):
|
|
|
239 Unknown = ""
|
|
|
240 Direct = "_F"
|
|
|
241 Inverse = "_B"
|
|
|
242
|
|
|
243 @classmethod
|
|
|
244 def fromDir(cls, s :str) -> "ReactionDirection":
|
|
|
245 # vvv as long as there's so few variants I actually condone the if spam:
|
|
|
246 if s == ReactionDirection.Direct.value: return ReactionDirection.Direct
|
|
|
247 if s == ReactionDirection.Inverse.value: return ReactionDirection.Inverse
|
|
|
248 return ReactionDirection.Unknown
|
|
|
249
|
|
|
250 @classmethod
|
|
|
251 def fromReactionId(cls, reactionId :str) -> "ReactionDirection":
|
|
|
252 return ReactionDirection.fromDir(reactionId[-2:])
|
|
|
253
|
|
|
254 def getArrowBodyElementId(reactionId :str) -> str:
|
|
|
255 if reactionId.endswith("_RV"): reactionId = reactionId[:-3] #TODO: standardize _RV
|
|
|
256 elif ReactionDirection.fromReactionId(reactionId) is not ReactionDirection.Unknown: reactionId = reactionId[:-2]
|
|
|
257 return f"R_{reactionId}"
|
|
|
258
|
|
|
259 def getArrowHeadElementId(reactionId :str) -> Tuple[str, str]:
|
|
|
260 """
|
|
|
261 We attempt extracting the direction information from the provided reaction ID, if unsuccessful we provide the IDs of both directions.
|
|
|
262
|
|
|
263 Args:
|
|
|
264 reactionId : the provided reaction ID.
|
|
|
265
|
|
|
266 Returns:
|
|
|
267 Tuple[str, str]: either a single str ID for the correct arrow head followed by an empty string or both options to try.
|
|
|
268 """
|
|
|
269 if reactionId.endswith("_RV"): reactionId = reactionId[:-3] #TODO: standardize _RV
|
|
|
270 elif ReactionDirection.fromReactionId(reactionId) is not ReactionDirection.Unknown: return reactionId[:-3:-1] + reactionId[:-2], ""
|
|
|
271 return f"F_{reactionId}", f"B_{reactionId}"
|
|
|
272
|
|
|
273 class ArrowColor(Enum):
|
|
|
274 """
|
|
|
275 Encodes possible arrow colors based on their meaning in the enrichment process.
|
|
|
276 """
|
|
|
277 Invalid = "#BEBEBE" # gray, fold-change under treshold or not significant p-value
|
|
|
278 Transparent = "#ffffff00" # transparent, to make some arrow segments disappear
|
|
|
279 UpRegulated = "#ecac68" # red, up-regulated reaction
|
|
|
280 DownRegulated = "#6495ed" # blue, down-regulated reaction
|
|
|
281
|
|
|
282 UpRegulatedInv = "#FF0000"
|
|
|
283 # ^^^ different shade of red (actually orange), up-regulated net value for a reversible reaction with
|
|
|
284 # conflicting enrichment in the two directions.
|
|
|
285
|
|
|
286 DownRegulatedInv = "#0000FF"
|
|
|
287 # ^^^ different shade of blue (actually purple), down-regulated net value for a reversible reaction with
|
|
|
288 # conflicting enrichment in the two directions.
|
|
|
289
|
|
|
290 @classmethod
|
|
|
291 def fromFoldChangeSign(cls, foldChange :float, *, useAltColor = False) -> "ArrowColor":
|
|
|
292 colors = (cls.DownRegulated, cls.DownRegulatedInv) if foldChange < 0 else (cls.UpRegulated, cls.UpRegulatedInv)
|
|
|
293 return colors[useAltColor]
|
|
|
294
|
|
|
295 def __str__(self) -> str: return self.value
|
|
|
296
|
|
|
297 class Arrow:
|
|
|
298 """
|
|
|
299 Models the properties of a reaction arrow that change based on enrichment.
|
|
|
300 """
|
|
|
301 MIN_W = 2
|
|
|
302 MAX_W = 12
|
|
|
303
|
|
|
304 def __init__(self, width :int, col: ArrowColor, *, isDashed = False) -> None:
|
|
|
305 """
|
|
|
306 (Private) Initializes an instance of Arrow.
|
|
|
307
|
|
|
308 Args:
|
|
|
309 width : width of the arrow, ideally to be kept within Arrow.MIN_W and Arrow.MAX_W (not enforced).
|
|
|
310 col : color of the arrow.
|
|
|
311 isDashed : whether the arrow should be dashed, meaning the associated pValue resulted not significant.
|
|
|
312
|
|
|
313 Returns:
|
|
|
314 None : practically, a Arrow instance.
|
|
|
315 """
|
|
|
316 self.w = width
|
|
|
317 self.col = col
|
|
|
318 self.dash = isDashed
|
|
|
319
|
|
|
320 def applyTo(self, reactionId :str, metabMap :ET.ElementTree, styleStr :str) -> None:
|
|
|
321 if getElementById(reactionId, metabMap).map(lambda el : styleMapElement(el, styleStr)).isErr:
|
|
|
322 ERRORS.append(reactionId)
|
|
|
323
|
|
|
324 def styleReactionElements(self, metabMap :ET.ElementTree, reactionId :str, *, mindReactionDir = True) -> None:
|
|
|
325 if not mindReactionDir:
|
|
|
326 return self.applyTo(getArrowBodyElementId(reactionId), metabMap, self.toStyleStr())
|
|
|
327
|
|
|
328 # Now we style the arrow head(s):
|
|
|
329 idOpt1, idOpt2 = getArrowHeadElementId(reactionId)
|
|
|
330 self.applyTo(idOpt1, metabMap, self.toStyleStr(downSizedForTips = True))
|
|
|
331 if idOpt2: self.applyTo(idOpt2, metabMap, self.toStyleStr(downSizedForTips = True))
|
|
|
332
|
|
|
333 def styleReactionElementsMeanMedian(self, metabMap :ET.ElementTree, reactionId :str, isNegative:bool) -> None:
|
|
|
334
|
|
|
335 self.applyTo(getArrowBodyElementId(reactionId), metabMap, self.toStyleStr())
|
|
|
336 idOpt1, idOpt2 = getArrowHeadElementId(reactionId)
|
|
|
337
|
|
|
338 if(isNegative):
|
|
|
339 self.applyTo(idOpt2, metabMap, self.toStyleStr(downSizedForTips = True))
|
|
|
340 self.col = ArrowColor.Transparent
|
|
|
341 self.applyTo(idOpt1, metabMap, self.toStyleStr(downSizedForTips = True)) #trasp
|
|
|
342 else:
|
|
|
343 self.applyTo(idOpt1, metabMap, self.toStyleStr(downSizedForTips = True))
|
|
|
344 self.col = ArrowColor.Transparent
|
|
|
345 self.applyTo(idOpt2, metabMap, self.toStyleStr(downSizedForTips = True)) #trasp
|
|
|
346
|
|
|
347
|
|
|
348
|
|
|
349 def getMapReactionId(self, reactionId :str, mindReactionDir :bool) -> str:
|
|
|
350 """
|
|
|
351 Computes the reaction ID as encoded in the map for a given reaction ID from the dataset.
|
|
|
352
|
|
|
353 Args:
|
|
|
354 reactionId: the reaction ID, as encoded in the dataset.
|
|
|
355 mindReactionDir: if True forward (F_) and backward (B_) directions will be encoded in the result.
|
|
|
356
|
|
|
357 Returns:
|
|
|
358 str : the ID of an arrow's body or tips in the map.
|
|
|
359 """
|
|
|
360 # we assume the reactionIds also don't encode reaction dir if they don't mind it when styling the map.
|
|
|
361 if not mindReactionDir: return "R_" + reactionId
|
|
|
362
|
|
|
363 #TODO: this is clearly something we need to make consistent in fluxes
|
|
|
364 return (reactionId[:-3:-1] + reactionId[:-2]) if reactionId[:-2] in ["_F", "_B"] else f"F_{reactionId}" # "Pyr_F" --> "F_Pyr"
|
|
|
365
|
|
|
366 def toStyleStr(self, *, downSizedForTips = False) -> str:
|
|
|
367 """
|
|
|
368 Collapses the styles of this Arrow into a str, ready to be applied as part of the "style" property on an svg element.
|
|
|
369
|
|
|
370 Returns:
|
|
|
371 str : the styles string.
|
|
|
372 """
|
|
|
373 width = self.w
|
|
|
374 if downSizedForTips: width *= 0.8
|
|
|
375 return f";stroke:{self.col};stroke-width:{width};stroke-dasharray:{'5,5' if self.dash else 'none'}"
|
|
|
376
|
|
|
377 # vvv These constants could be inside the class itself a static properties, but python
|
|
|
378 # was built by brainless organisms so here we are!
|
|
|
379 INVALID_ARROW = Arrow(Arrow.MIN_W, ArrowColor.Invalid)
|
|
|
380 INSIGNIFICANT_ARROW = Arrow(Arrow.MIN_W, ArrowColor.Invalid, isDashed = True)
|
|
|
381
|
|
|
382 def applyFluxesEnrichmentToMap(fluxesEnrichmentRes :Dict[str, Union[Tuple[float, FoldChange], Tuple[float, FoldChange, float, float]]], metabMap :ET.ElementTree, maxNumericZScore :float) -> None:
|
|
|
383 """
|
|
|
384 Applies fluxes enrichment results to the provided metabolic map.
|
|
|
385
|
|
|
386 Args:
|
|
|
387 fluxesEnrichmentRes : fluxes enrichment results.
|
|
|
388 metabMap : the metabolic map to edit.
|
|
|
389 maxNumericZScore : biggest finite z-score value found.
|
|
|
390
|
|
|
391 Side effects:
|
|
|
392 metabMap : mut
|
|
|
393
|
|
|
394 Returns:
|
|
|
395 None
|
|
|
396 """
|
|
|
397 for reactionId, values in fluxesEnrichmentRes.items():
|
|
|
398 pValue = values[0]
|
|
|
399 foldChange = values[1]
|
|
|
400 z_score = values[2]
|
|
|
401
|
|
|
402 if math.isnan(pValue) or (isinstance(foldChange, float) and math.isnan(foldChange)):
|
|
|
403 continue
|
|
|
404
|
|
|
405 if isinstance(foldChange, str): foldChange = float(foldChange)
|
|
|
406 if pValue > ARGS.pValue: # pValue above tresh: dashed arrow
|
|
|
407 INSIGNIFICANT_ARROW.styleReactionElements(metabMap, reactionId)
|
|
|
408 INSIGNIFICANT_ARROW.styleReactionElements(metabMap, reactionId, mindReactionDir = False)
|
|
|
409
|
|
|
410 continue
|
|
|
411
|
|
|
412 if abs(foldChange) < (ARGS.fChange - 1) / (abs(ARGS.fChange) + 1):
|
|
|
413 INVALID_ARROW.styleReactionElements(metabMap, reactionId)
|
|
|
414 INVALID_ARROW.styleReactionElements(metabMap, reactionId, mindReactionDir = False)
|
|
|
415
|
|
|
416 continue
|
|
|
417
|
|
|
418 width = Arrow.MAX_W
|
|
|
419 if not math.isinf(z_score):
|
|
|
420 try:
|
|
|
421 width = min(
|
|
|
422 max(abs(z_score * Arrow.MAX_W) / maxNumericZScore, Arrow.MIN_W),
|
|
|
423 Arrow.MAX_W)
|
|
|
424
|
|
|
425 except ZeroDivisionError: pass
|
|
|
426 # TODO CHECK RV
|
|
|
427 #if not reactionId.endswith("_RV"): # RV stands for reversible reactions
|
|
|
428 # Arrow(width, ArrowColor.fromFoldChangeSign(foldChange)).styleReactionElements(metabMap, reactionId)
|
|
|
429 # continue
|
|
|
430
|
|
|
431 #reactionId = reactionId[:-3] # Remove "_RV"
|
|
|
432
|
|
|
433 inversionScore = (values[3] < 0) + (values[4] < 0) # Compacts the signs of averages into 1 easy to check score
|
|
|
434 if inversionScore == 2: foldChange *= -1
|
|
|
435 # ^^^ Style the inverse direction with the opposite sign netValue
|
|
|
436
|
|
|
437 # If the score is 1 (opposite signs) we use alternative colors vvv
|
|
|
438 arrow = Arrow(width, ArrowColor.fromFoldChangeSign(foldChange, useAltColor = inversionScore == 1))
|
|
|
439
|
|
|
440 # vvv These 2 if statements can both be true and can both happen
|
|
|
441 if ARGS.net: # style arrow head(s):
|
|
|
442 arrow.styleReactionElements(metabMap, reactionId + ("_B" if inversionScore == 2 else "_F"))
|
|
|
443 arrow.applyTo(("F_" if inversionScore == 2 else "B_") + reactionId, metabMap, f";stroke:{ArrowColor.Transparent};stroke-width:0;stroke-dasharray:None")
|
|
|
444
|
|
|
445 arrow.styleReactionElements(metabMap, reactionId, mindReactionDir = False)
|
|
|
446
|
|
|
447
|
|
|
448 ############################ split class ######################################
|
|
|
449 def split_class(classes :pd.DataFrame, resolve_rules :Dict[str, List[float]]) -> Dict[str, List[List[float]]]:
|
|
|
450 """
|
|
|
451 Generates a :dict that groups together data from a :DataFrame based on classes the data is related to.
|
|
|
452
|
|
|
453 Args:
|
|
|
454 classes : a :DataFrame of only string values, containing class information (rows) and keys to query the resolve_rules :dict
|
|
|
455 resolve_rules : a :dict containing :float data
|
|
|
456
|
|
|
457 Returns:
|
|
|
458 dict : the dict with data grouped by class
|
|
|
459
|
|
|
460 Side effects:
|
|
|
461 classes : mut
|
|
|
462 """
|
|
|
463 class_pat :Dict[str, List[List[float]]] = {}
|
|
|
464 for i in range(len(classes)):
|
|
|
465 classe :str = classes.iloc[i, 1]
|
|
|
466 if pd.isnull(classe): continue
|
|
|
467
|
|
|
468 l :List[List[float]] = []
|
|
|
469 for j in range(i, len(classes)):
|
|
|
470 if classes.iloc[j, 1] == classe:
|
|
|
471 pat_id :str = classes.iloc[j, 0]
|
|
|
472 tmp = resolve_rules.get(pat_id, None)
|
|
|
473 if tmp != None:
|
|
|
474 l.append(tmp)
|
|
|
475 classes.iloc[j, 1] = None
|
|
|
476
|
|
|
477 if l:
|
|
|
478 class_pat[classe] = list(map(list, zip(*l)))
|
|
|
479 continue
|
|
|
480
|
|
|
481 utils.logWarning(
|
|
|
482 f"Warning: no sample found in class \"{classe}\", the class has been disregarded", ARGS.out_log)
|
|
|
483
|
|
|
484 return class_pat
|
|
|
485
|
|
|
486 ############################ conversion ##############################################
|
|
|
487 #conversion from svg to png
|
|
|
488 def svg_to_png_with_background(svg_path :utils.FilePath, png_path :utils.FilePath, dpi :int = 72, scale :int = 1, size :Optional[float] = None) -> None:
|
|
|
489 """
|
|
|
490 Internal utility to convert an SVG to PNG (forced opaque) to aid in PDF conversion.
|
|
|
491
|
|
|
492 Args:
|
|
|
493 svg_path : path to SVG file
|
|
|
494 png_path : path for new PNG file
|
|
|
495 dpi : dots per inch of the generated PNG
|
|
|
496 scale : scaling factor for the generated PNG, computed internally when a size is provided
|
|
|
497 size : final effective width of the generated PNG
|
|
|
498
|
|
|
499 Returns:
|
|
|
500 None
|
|
|
501 """
|
|
|
502 if size:
|
|
|
503 image = pyvips.Image.new_from_file(svg_path.show(), dpi=dpi, scale=1)
|
|
|
504 scale = size / image.width
|
|
|
505 image = image.resize(scale)
|
|
|
506 else:
|
|
|
507 image = pyvips.Image.new_from_file(svg_path.show(), dpi=dpi, scale=scale)
|
|
|
508
|
|
|
509 white_background = pyvips.Image.black(image.width, image.height).new_from_image([255, 255, 255])
|
|
|
510 white_background = white_background.affine([scale, 0, 0, scale])
|
|
|
511
|
|
|
512 if white_background.bands != image.bands:
|
|
|
513 white_background = white_background.extract_band(0)
|
|
|
514
|
|
|
515 composite_image = white_background.composite2(image, 'over')
|
|
|
516 composite_image.write_to_file(png_path.show())
|
|
|
517
|
|
|
518 #funzione unica, lascio fuori i file e li passo in input
|
|
|
519 #conversion from png to pdf
|
|
|
520 def convert_png_to_pdf(png_file :utils.FilePath, pdf_file :utils.FilePath) -> None:
|
|
|
521 """
|
|
|
522 Internal utility to convert a PNG to PDF to aid from SVG conversion.
|
|
|
523
|
|
|
524 Args:
|
|
|
525 png_file : path to PNG file
|
|
|
526 pdf_file : path to new PDF file
|
|
|
527
|
|
|
528 Returns:
|
|
|
529 None
|
|
|
530 """
|
|
|
531 image = Image.open(png_file.show())
|
|
|
532 image = image.convert("RGB")
|
|
|
533 image.save(pdf_file.show(), "PDF", resolution=100.0)
|
|
|
534
|
|
|
535 #function called to reduce redundancy in the code
|
|
|
536 def convert_to_pdf(file_svg :utils.FilePath, file_png :utils.FilePath, file_pdf :utils.FilePath) -> None:
|
|
|
537 """
|
|
|
538 Converts the SVG map at the provided path to PDF.
|
|
|
539
|
|
|
540 Args:
|
|
|
541 file_svg : path to SVG file
|
|
|
542 file_png : path to PNG file
|
|
|
543 file_pdf : path to new PDF file
|
|
|
544
|
|
|
545 Returns:
|
|
|
546 None
|
|
|
547 """
|
|
|
548 svg_to_png_with_background(file_svg, file_png)
|
|
|
549 try:
|
|
|
550 convert_png_to_pdf(file_png, file_pdf)
|
|
|
551 print(f'PDF file {file_pdf.filePath} successfully generated.')
|
|
|
552
|
|
|
553 except Exception as e:
|
|
|
554 raise utils.DataErr(file_pdf.show(), f'Error generating PDF file: {e}')
|
|
|
555
|
|
|
556 ############################ map ##############################################
|
|
|
557 def buildOutputPath(dataset1Name :str, dataset2Name = "rest", *, details = "", ext :utils.FileFormat) -> utils.FilePath:
|
|
|
558 """
|
|
|
559 Builds a FilePath instance from the names of confronted datasets ready to point to a location in the
|
|
|
560 "result/" folder, used by this tool for output files in collections.
|
|
|
561
|
|
|
562 Args:
|
|
|
563 dataset1Name : _description_
|
|
|
564 dataset2Name : _description_. Defaults to "rest".
|
|
|
565 details : _description_
|
|
|
566 ext : _description_
|
|
|
567
|
|
|
568 Returns:
|
|
|
569 utils.FilePath : _description_
|
|
|
570 """
|
|
|
571 # This function returns a util data structure but is extremely specific to this module.
|
|
|
572 # RAS also uses collections as output and as such might benefit from a method like this, but I'd wait
|
|
|
573 # TODO: until a third tool with multiple outputs appears before porting this to utils.
|
|
|
574 return utils.FilePath(
|
|
|
575 f"{dataset1Name}_vs_{dataset2Name}" + (f" ({details})" if details else ""),
|
|
|
576 # ^^^ yes this string is built every time even if the form is the same for the same 2 datasets in
|
|
|
577 # all output files: I don't care, this was never the performance bottleneck of the tool and
|
|
|
578 # there is no other net gain in saving and re-using the built string.
|
|
|
579 ext,
|
|
|
580 prefix = ARGS.output_path)
|
|
|
581
|
|
|
582 FIELD_NOT_AVAILABLE = '/'
|
|
|
583 def writeToCsv(rows: List[list], fieldNames :List[str], outPath :utils.FilePath) -> None:
|
|
|
584 fieldsAmt = len(fieldNames)
|
|
|
585 with open(outPath.show(), "w", newline = "") as fd:
|
|
|
586 writer = csv.DictWriter(fd, fieldnames = fieldNames, delimiter = '\t')
|
|
|
587 writer.writeheader()
|
|
|
588
|
|
|
589 for row in rows:
|
|
|
590 sizeMismatch = fieldsAmt - len(row)
|
|
|
591 if sizeMismatch > 0: row.extend([FIELD_NOT_AVAILABLE] * sizeMismatch)
|
|
|
592 writer.writerow({ field : data for field, data in zip(fieldNames, row) })
|
|
|
593
|
|
|
594 OldEnrichedScores = Dict[str, List[Union[float, FoldChange]]] #TODO: try to use Tuple whenever possible
|
|
|
595 def writeTabularResult(enrichedScores : OldEnrichedScores, outPath :utils.FilePath) -> None:
|
|
|
596 fieldNames = ["ids", "P_Value", "fold change", "z-score"]
|
|
|
597 fieldNames.extend(["average_1", "average_2"])
|
|
|
598
|
|
|
599 writeToCsv([ [reactId] + values for reactId, values in enrichedScores.items() ], fieldNames, outPath)
|
|
|
600
|
|
|
601 def temp_thingsInCommon(tmp :Dict[str, List[Union[float, FoldChange]]], core_map :ET.ElementTree, max_z_score :float, dataset1Name :str, dataset2Name = "rest") -> None:
|
|
|
602 # this function compiles the things always in common between comparison modes after enrichment.
|
|
|
603 # TODO: organize, name better.
|
|
|
604 writeTabularResult(tmp, buildOutputPath(dataset1Name, dataset2Name, details = "Tabular Result", ext = utils.FileFormat.TSV))
|
|
|
605 for reactId, enrichData in tmp.items(): tmp[reactId] = tuple(enrichData)
|
|
|
606 applyFluxesEnrichmentToMap(tmp, core_map, max_z_score)
|
|
|
607
|
|
|
608 def computePValue(dataset1Data: List[float], dataset2Data: List[float]) -> Tuple[float, float]:
|
|
|
609 """
|
|
|
610 Computes the statistical significance score (P-value) of the comparison between coherent data
|
|
|
611 from two datasets. The data is supposed to, in both datasets:
|
|
|
612 - be related to the same reaction ID;
|
|
|
613 - be ordered by sample, such that the item at position i in both lists is related to the
|
|
|
614 same sample or cell line.
|
|
|
615
|
|
|
616 Args:
|
|
|
617 dataset1Data : data from the 1st dataset.
|
|
|
618 dataset2Data : data from the 2nd dataset.
|
|
|
619
|
|
|
620 Returns:
|
|
|
621 tuple: (P-value, Z-score)
|
|
|
622 - P-value from the selected test on the provided data.
|
|
|
623 - Z-score of the difference between means of the two datasets.
|
|
|
624 """
|
|
|
625
|
|
|
626 match ARGS.test:
|
|
|
627 case "ks":
|
|
|
628 # Perform Kolmogorov-Smirnov test
|
|
|
629 _, p_value = st.ks_2samp(dataset1Data, dataset2Data)
|
|
|
630 case "ttest_p":
|
|
|
631 # Datasets should have same size
|
|
|
632 if len(dataset1Data) != len(dataset2Data):
|
|
|
633 raise ValueError("Datasets must have the same size for paired t-test.")
|
|
|
634 # Perform t-test for paired samples
|
|
|
635 _, p_value = st.ttest_rel(dataset1Data, dataset2Data)
|
|
|
636 case "ttest_ind":
|
|
|
637 # Perform t-test for independent samples
|
|
|
638 _, p_value = st.ttest_ind(dataset1Data, dataset2Data)
|
|
|
639 case "wilcoxon":
|
|
|
640 # Datasets should have same size
|
|
|
641 if len(dataset1Data) != len(dataset2Data):
|
|
|
642 raise ValueError("Datasets must have the same size for Wilcoxon signed-rank test.")
|
|
|
643 # Perform Wilcoxon signed-rank test
|
|
|
644 np.random.seed(42) # Ensure reproducibility since zsplit method is used
|
|
|
645 _, p_value = st.wilcoxon(dataset1Data, dataset2Data, zero_method="zsplit")
|
|
|
646 case "mw":
|
|
|
647 # Perform Mann-Whitney U test
|
|
|
648 _, p_value = st.mannwhitneyu(dataset1Data, dataset2Data)
|
|
|
649
|
|
|
650 # Calculate means and standard deviations
|
|
|
651 mean1 = np.nanmean(dataset1Data)
|
|
|
652 mean2 = np.nanmean(dataset2Data)
|
|
|
653 std1 = np.nanstd(dataset1Data, ddof=1)
|
|
|
654 std2 = np.nanstd(dataset2Data, ddof=1)
|
|
|
655
|
|
|
656 n1 = len(dataset1Data)
|
|
|
657 n2 = len(dataset2Data)
|
|
|
658
|
|
|
659 # Calculate Z-score
|
|
|
660 z_score = (mean1 - mean2) / np.sqrt((std1**2 / n1) + (std2**2 / n2))
|
|
|
661
|
|
|
662 return p_value, z_score
|
|
|
663
|
|
|
664 def compareDatasetPair(dataset1Data :List[List[float]], dataset2Data :List[List[float]], ids :List[str]) -> Tuple[Dict[str, List[Union[float, FoldChange]]], float]:
|
|
|
665 #TODO: the following code still suffers from "dumbvarnames-osis"
|
|
|
666 comparisonResult :Dict[str, List[Union[float, FoldChange]]] = {}
|
|
|
667 count = 0
|
|
|
668 max_z_score = 0
|
|
|
669 for l1, l2 in zip(dataset1Data, dataset2Data):
|
|
|
670 reactId = ids[count]
|
|
|
671 count += 1
|
|
|
672 if not reactId: continue # we skip ids that have already been processed
|
|
|
673
|
|
|
674 try:
|
|
|
675 p_value, z_score = computePValue(l1, l2)
|
|
|
676 avg1 = sum(l1) / len(l1)
|
|
|
677 avg2 = sum(l2) / len(l2)
|
|
|
678 f_c = fold_change(avg1, avg2)
|
|
|
679 if np.isfinite(z_score) and max_z_score < abs(z_score): max_z_score = abs(z_score)
|
|
|
680
|
|
|
681 comparisonResult[reactId] = [float(p_value), f_c, z_score, avg1, avg2]
|
|
|
682 except (TypeError, ZeroDivisionError): continue
|
|
|
683
|
|
|
684 # Apply multiple testing correction if set by the user
|
|
|
685 if ARGS.adjusted:
|
|
|
686
|
|
|
687 # Retrieve the p-values from the comparisonResult dictionary, they have to be different from NaN
|
|
|
688 validPValues = [(reactId, result[0]) for reactId, result in comparisonResult.items() if not np.isnan(result[0])]
|
|
|
689
|
|
|
690 if not validPValues:
|
|
|
691 return comparisonResult, max_z_score
|
|
|
692
|
|
|
693 # Unpack the valid p-values
|
|
|
694 reactIds, pValues = zip(*validPValues)
|
|
|
695 # Adjust the p-values using the Benjamini-Hochberg method
|
|
|
696 adjustedPValues = st.false_discovery_control(pValues)
|
|
|
697 # Update the comparisonResult dictionary with the adjusted p-values
|
|
|
698 for reactId , adjustedPValue in zip(reactIds, adjustedPValues):
|
|
|
699 comparisonResult[reactId][0] = adjustedPValue
|
|
|
700
|
|
|
701 return comparisonResult, max_z_score
|
|
|
702
|
|
|
703 def computeEnrichment(class_pat :Dict[str, List[List[float]]], ids :List[str]) -> List[Tuple[str, str, dict, float]]:
|
|
|
704 """
|
|
|
705 Compares clustered data based on a given comparison mode and applies enrichment-based styling on the
|
|
|
706 provided metabolic map.
|
|
|
707
|
|
|
708 Args:
|
|
|
709 class_pat : the clustered data.
|
|
|
710 ids : ids for data association.
|
|
|
711
|
|
|
712
|
|
|
713 Returns:
|
|
|
714 List[Tuple[str, str, dict, float]]: List of tuples with pairs of dataset names, comparison dictionary, and max z-score.
|
|
|
715
|
|
|
716 Raises:
|
|
|
717 sys.exit : if there are less than 2 classes for comparison
|
|
|
718
|
|
|
719 """
|
|
|
720 class_pat = { k.strip() : v for k, v in class_pat.items() }
|
|
|
721 #TODO: simplfy this stuff vvv and stop using sys.exit (raise the correct utils error)
|
|
|
722 if (not class_pat) or (len(class_pat.keys()) < 2): sys.exit('Execution aborted: classes provided for comparisons are less than two\n')
|
|
|
723
|
|
|
724 enrichment_results = []
|
|
|
725
|
|
|
726
|
|
|
727 if ARGS.comparison == "manyvsmany":
|
|
|
728 for i, j in it.combinations(class_pat.keys(), 2):
|
|
|
729 comparisonDict, max_z_score = compareDatasetPair(class_pat.get(i), class_pat.get(j), ids)
|
|
|
730 enrichment_results.append((i, j, comparisonDict, max_z_score))
|
|
|
731
|
|
|
732 elif ARGS.comparison == "onevsrest":
|
|
|
733 for single_cluster in class_pat.keys():
|
|
|
734 rest = [item for k, v in class_pat.items() if k != single_cluster for item in v]
|
|
|
735
|
|
|
736 comparisonDict, max_z_score = compareDatasetPair(class_pat.get(single_cluster), rest, ids)
|
|
|
737 enrichment_results.append((single_cluster, "rest", comparisonDict, max_z_score))
|
|
|
738
|
|
|
739 #elif ARGS.comparison == "onevsmany":
|
|
|
740 # controlItems = class_pat.get(ARGS.control)
|
|
|
741 # for otherDataset in class_pat.keys():
|
|
|
742 # if otherDataset == ARGS.control:
|
|
|
743 # continue
|
|
|
744 # comparisonDict, max_z_score = compareDatasetPair(controlItems, class_pat.get(otherDataset), ids)
|
|
|
745 # enrichment_results.append((ARGS.control, otherDataset, comparisonDict, max_z_score))
|
|
|
746 elif ARGS.comparison == "onevsmany":
|
|
|
747 controlItems = class_pat.get(ARGS.control)
|
|
|
748 for otherDataset in class_pat.keys():
|
|
|
749 if otherDataset == ARGS.control:
|
|
|
750 continue
|
|
|
751 comparisonDict, max_z_score = compareDatasetPair(class_pat.get(otherDataset),controlItems, ids)
|
|
|
752 enrichment_results.append(( otherDataset,ARGS.control, comparisonDict, max_z_score))
|
|
|
753
|
|
|
754 return enrichment_results
|
|
|
755
|
|
|
756 def createOutputMaps(dataset1Name :str, dataset2Name :str, core_map :ET.ElementTree) -> None:
|
|
|
757 svgFilePath = buildOutputPath(dataset1Name, dataset2Name, details="SVG Map", ext=utils.FileFormat.SVG)
|
|
|
758 utils.writeSvg(svgFilePath, core_map)
|
|
|
759
|
|
|
760 if ARGS.generate_pdf:
|
|
|
761 pngPath = buildOutputPath(dataset1Name, dataset2Name, details="PNG Map", ext=utils.FileFormat.PNG)
|
|
|
762 pdfPath = buildOutputPath(dataset1Name, dataset2Name, details="PDF Map", ext=utils.FileFormat.PDF)
|
|
|
763 convert_to_pdf(svgFilePath, pngPath, pdfPath)
|
|
|
764
|
|
|
765 if not ARGS.generate_svg:
|
|
|
766 os.remove(svgFilePath.show())
|
|
|
767
|
|
|
768 ClassPat = Dict[str, List[List[float]]]
|
|
|
769 def getClassesAndIdsFromDatasets(datasetsPaths :List[str], datasetPath :str, classPath :str, names :List[str]) -> Tuple[List[str], ClassPat]:
|
|
|
770 # TODO: I suggest creating dicts with ids as keys instead of keeping class_pat and ids separate,
|
|
|
771 # for the sake of everyone's sanity.
|
|
|
772 class_pat :ClassPat = {}
|
|
|
773 if ARGS.option == 'datasets':
|
|
|
774 num = 1 #TODO: the dataset naming function could be a generator
|
|
|
775 for path, name in zip(datasetsPaths, names):
|
|
|
776 name = name_dataset(name, num)
|
|
|
777 resolve_rules_float, ids = getDatasetValues(path, name)
|
|
|
778 if resolve_rules_float != None:
|
|
|
779 class_pat[name] = list(map(list, zip(*resolve_rules_float.values())))
|
|
|
780
|
|
|
781 num += 1
|
|
|
782
|
|
|
783 elif ARGS.option == "dataset_class":
|
|
|
784 classes = read_dataset(classPath, "class")
|
|
|
785 classes = classes.astype(str)
|
|
|
786 resolve_rules_float, ids = getDatasetValues(datasetPath, "Dataset Class (not actual name)")
|
|
|
787 #check if classes have match on ids
|
|
|
788 if not all(classes.iloc[:, 0].isin(ids)):
|
|
|
789 utils.logWarning(
|
|
|
790 "No match between classes and sample IDs", ARGS.out_log)
|
|
|
791 if resolve_rules_float != None: class_pat = split_class(classes, resolve_rules_float)
|
|
|
792
|
|
|
793 return ids, class_pat
|
|
|
794 #^^^ TODO: this could be a match statement over an enum, make it happen future marea dev with python 3.12! (it's why I kept the ifs)
|
|
|
795
|
|
|
796 #TODO: create these damn args as FilePath objects
|
|
|
797 def getDatasetValues(datasetPath :str, datasetName :str) -> Tuple[ClassPat, List[str]]:
|
|
|
798 """
|
|
|
799 Opens the dataset at the given path and extracts the values (expected nullable numerics) and the IDs.
|
|
|
800
|
|
|
801 Args:
|
|
|
802 datasetPath : path to the dataset
|
|
|
803 datasetName (str): dataset name, used in error reporting
|
|
|
804
|
|
|
805 Returns:
|
|
|
806 Tuple[ClassPat, List[str]]: values and IDs extracted from the dataset
|
|
|
807 """
|
|
|
808 dataset = read_dataset(datasetPath, datasetName)
|
|
|
809
|
|
|
810 # Ensure the first column is treated as the reaction name
|
|
|
811 dataset = dataset.set_index(dataset.columns[0])
|
|
|
812
|
|
|
813 # Check if required reactions exist in the dataset
|
|
|
814 required_reactions = ['EX_lac__L_e', 'EX_glc__D_e', 'EX_gln__L_e', 'EX_glu__L_e']
|
|
|
815 missing_reactions = [reaction for reaction in required_reactions if reaction not in dataset.index]
|
|
|
816
|
|
|
817 if missing_reactions:
|
|
|
818 sys.exit(f'Execution aborted: Missing required reactions {missing_reactions} in {datasetName}\n')
|
|
|
819
|
|
|
820 # Calculate new rows using safe division
|
|
|
821 lact_glc = np.divide(
|
|
|
822 np.clip(dataset.loc['EX_lac__L_e'].to_numpy(), a_min=0, a_max=None),
|
|
|
823 np.clip(dataset.loc['EX_glc__D_e'].to_numpy(), a_min=None, a_max=0),
|
|
|
824 out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan), # Prepara un array con NaN come output di default
|
|
|
825 where=dataset.loc['EX_glc__D_e'].to_numpy() != 0 # Condizione per evitare la divisione per zero
|
|
|
826 )
|
|
|
827 lact_gln = np.divide(
|
|
|
828 np.clip(dataset.loc['EX_lac__L_e'].to_numpy(), a_min=0, a_max=None),
|
|
|
829 np.clip(dataset.loc['EX_gln__L_e'].to_numpy(), a_min=None, a_max=0),
|
|
|
830 out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan),
|
|
|
831 where=dataset.loc['EX_gln__L_e'].to_numpy() != 0
|
|
|
832 )
|
|
|
833 lact_o2 = np.divide(
|
|
|
834 np.clip(dataset.loc['EX_lac__L_e'].to_numpy(), a_min=0, a_max=None),
|
|
|
835 np.clip(dataset.loc['EX_o2_e'].to_numpy(), a_min=None, a_max=0),
|
|
|
836 out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan),
|
|
|
837 where=dataset.loc['EX_o2_e'].to_numpy() != 0
|
|
|
838 )
|
|
|
839 glu_gln = np.divide(
|
|
|
840 dataset.loc['EX_glu__L_e'].to_numpy(),
|
|
|
841 np.clip(dataset.loc['EX_gln__L_e'].to_numpy(), a_min=None, a_max=0),
|
|
|
842 out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan),
|
|
|
843 where=dataset.loc['EX_gln__L_e'].to_numpy() != 0
|
|
|
844 )
|
|
|
845
|
|
|
846
|
|
|
847 values = {'lact_glc': lact_glc, 'lact_gln': lact_gln, 'lact_o2': lact_o2, 'glu_gln': glu_gln}
|
|
|
848
|
|
|
849 # Sostituzione di inf e NaN con 0 se necessario
|
|
|
850 for key in values:
|
|
|
851 values[key] = np.nan_to_num(values[key], nan=0.0, posinf=0.0, neginf=0.0)
|
|
|
852
|
|
|
853 # Creazione delle nuove righe da aggiungere al dataset
|
|
|
854 new_rows = pd.DataFrame({
|
|
|
855 dataset.index.name: ['LactGlc', 'LactGln', 'LactO2', 'GluGln'],
|
|
|
856 **{col: [values['lact_glc'][i], values['lact_gln'][i], values['lact_o2'][i], values['glu_gln'][i]]
|
|
|
857 for i, col in enumerate(dataset.columns)}
|
|
|
858 })
|
|
|
859
|
|
|
860 #print(new_rows)
|
|
|
861
|
|
|
862 # Ritorna il dataset originale con le nuove righe
|
|
|
863 dataset.reset_index(inplace=True)
|
|
|
864 dataset = pd.concat([dataset, new_rows], ignore_index=True)
|
|
|
865
|
|
|
866 IDs = pd.Series.tolist(dataset.iloc[:, 0].astype(str))
|
|
|
867
|
|
|
868 dataset = dataset.drop(dataset.columns[0], axis = "columns").to_dict("list")
|
|
|
869 return { id : list(map(utils.Float("Dataset values, not an argument"), values)) for id, values in dataset.items() }, IDs
|
|
|
870
|
|
|
871 def rgb_to_hex(rgb):
|
|
|
872 """
|
|
|
873 Convert RGB values (0-1 range) to hexadecimal color format.
|
|
|
874
|
|
|
875 Args:
|
|
|
876 rgb (numpy.ndarray): An array of RGB color components (in the range [0, 1]).
|
|
|
877
|
|
|
878 Returns:
|
|
|
879 str: The color in hexadecimal format (e.g., '#ff0000' for red).
|
|
|
880 """
|
|
|
881 # Convert RGB values (0-1 range) to hexadecimal format
|
|
|
882 rgb = (np.array(rgb) * 255).astype(int)
|
|
|
883 return '#{:02x}{:02x}{:02x}'.format(rgb[0], rgb[1], rgb[2])
|
|
|
884
|
|
|
885 def save_colormap_image(min_value: float, max_value: float, path: utils.FilePath, colorMap:str="viridis"):
|
|
|
886 """
|
|
|
887 Create and save an image of the colormap showing the gradient and its range.
|
|
|
888
|
|
|
889 Args:
|
|
|
890 min_value (float): The minimum value of the colormap range.
|
|
|
891 max_value (float): The maximum value of the colormap range.
|
|
|
892 filename (str): The filename for saving the image.
|
|
|
893 """
|
|
|
894
|
|
|
895 # Create a colormap using matplotlib
|
|
|
896 cmap = plt.get_cmap(colorMap)
|
|
|
897
|
|
|
898 # Create a figure and axis
|
|
|
899 fig, ax = plt.subplots(figsize=(6, 1))
|
|
|
900 fig.subplots_adjust(bottom=0.5)
|
|
|
901
|
|
|
902 # Create a gradient image
|
|
|
903 gradient = np.linspace(0, 1, 256)
|
|
|
904 gradient = np.vstack((gradient, gradient))
|
|
|
905
|
|
|
906 # Add min and max value annotations
|
|
|
907 ax.text(0, 0.5, f'{np.round(min_value, 3)}', va='center', ha='right', transform=ax.transAxes, fontsize=12, color='black')
|
|
|
908 ax.text(1, 0.5, f'{np.round(max_value, 3)}', va='center', ha='left', transform=ax.transAxes, fontsize=12, color='black')
|
|
|
909
|
|
|
910
|
|
|
911 # Display the gradient image
|
|
|
912 ax.imshow(gradient, aspect='auto', cmap=cmap)
|
|
|
913 ax.set_axis_off()
|
|
|
914
|
|
|
915 # Save the image
|
|
|
916 plt.savefig(path.show(), bbox_inches='tight', pad_inches=0)
|
|
|
917 plt.close()
|
|
|
918 pass
|
|
|
919
|
|
|
920 def min_nonzero_abs(arr):
|
|
|
921 # Flatten the array and filter out zeros, then find the minimum of the remaining values
|
|
|
922 non_zero_elements = np.abs(arr)[np.abs(arr) > 0]
|
|
|
923 return np.min(non_zero_elements) if non_zero_elements.size > 0 else None
|
|
|
924
|
|
|
925 def computeEnrichmentMeanMedian(metabMap: ET.ElementTree, class_pat: Dict[str, List[List[float]]], ids: List[str], colormap:str) -> None:
|
|
|
926 """
|
|
|
927 Compute and visualize the metabolic map based on mean and median of the input fluxes.
|
|
|
928 The fluxes are normalised across classes/datasets and visualised using the given colormap.
|
|
|
929
|
|
|
930 Args:
|
|
|
931 metabMap (ET.ElementTree): An XML tree representing the metabolic map.
|
|
|
932 class_pat (Dict[str, List[List[float]]]): A dictionary where keys are class names and values are lists of enrichment values.
|
|
|
933 ids (List[str]): A list of reaction IDs to be used for coloring arrows.
|
|
|
934
|
|
|
935 Returns:
|
|
|
936 None
|
|
|
937 """
|
|
|
938 # Create copies only if they are needed
|
|
|
939 metabMap_mean = copy.deepcopy(metabMap)
|
|
|
940 metabMap_median = copy.deepcopy(metabMap)
|
|
|
941
|
|
|
942 # Compute medians and means
|
|
|
943 medians = {key: np.round(np.nanmedian(np.array(value), axis=1), 6) for key, value in class_pat.items()}
|
|
|
944 means = {key: np.round(np.nanmean(np.array(value), axis=1),6) for key, value in class_pat.items()}
|
|
|
945
|
|
|
946 # Normalize medians and means
|
|
|
947 max_flux_medians = max(np.max(np.abs(arr)) for arr in medians.values())
|
|
|
948 max_flux_means = max(np.max(np.abs(arr)) for arr in means.values())
|
|
|
949
|
|
|
950 min_flux_medians = min(min_nonzero_abs(arr) for arr in medians.values())
|
|
|
951 min_flux_means = min(min_nonzero_abs(arr) for arr in means.values())
|
|
|
952
|
|
|
953 medians = {key: median/max_flux_medians for key, median in medians.items()}
|
|
|
954 means = {key: mean/max_flux_means for key, mean in means.items()}
|
|
|
955
|
|
|
956 save_colormap_image(min_flux_medians, max_flux_medians, utils.FilePath("Color map median", ext=utils.FileFormat.PNG, prefix=ARGS.output_path), colormap)
|
|
|
957 save_colormap_image(min_flux_means, max_flux_means, utils.FilePath("Color map mean", ext=utils.FileFormat.PNG, prefix=ARGS.output_path), colormap)
|
|
|
958
|
|
|
959 cmap = plt.get_cmap(colormap)
|
|
|
960
|
|
|
961 min_width = 2.0 # Minimum arrow width
|
|
|
962 max_width = 15.0 # Maximum arrow width
|
|
|
963
|
|
|
964 for key in class_pat:
|
|
|
965 # Create color mappings for median and mean
|
|
|
966 colors_median = {
|
|
|
967 rxn_id: rgb_to_hex(cmap(abs(medians[key][i]))) if medians[key][i] != 0 else '#bebebe' #grey blocked
|
|
|
968 for i, rxn_id in enumerate(ids)
|
|
|
969 }
|
|
|
970
|
|
|
971 colors_mean = {
|
|
|
972 rxn_id: rgb_to_hex(cmap(abs(means[key][i]))) if means[key][i] != 0 else '#bebebe' #grey blocked
|
|
|
973 for i, rxn_id in enumerate(ids)
|
|
|
974 }
|
|
|
975
|
|
|
976 for i, rxn_id in enumerate(ids):
|
|
|
977 # Calculate arrow width for median
|
|
|
978 width_median = np.interp(abs(medians[key][i]), [0, 1], [min_width, max_width])
|
|
|
979 isNegative = medians[key][i] < 0
|
|
|
980 apply_arrow(metabMap_median, rxn_id, colors_median[rxn_id], isNegative, width_median)
|
|
|
981
|
|
|
982 # Calculate arrow width for mean
|
|
|
983 width_mean = np.interp(abs(means[key][i]), [0, 1], [min_width, max_width])
|
|
|
984 isNegative = means[key][i] < 0
|
|
|
985 apply_arrow(metabMap_mean, rxn_id, colors_mean[rxn_id], isNegative, width_mean)
|
|
|
986
|
|
|
987 # Save and convert the SVG files
|
|
|
988 save_and_convert(metabMap_mean, "mean", key)
|
|
|
989 save_and_convert(metabMap_median, "median", key)
|
|
|
990
|
|
|
991 def apply_arrow(metabMap, rxn_id, color, isNegative, width=5):
|
|
|
992 """
|
|
|
993 Apply an arrow to a specific reaction in the metabolic map with a given color.
|
|
|
994
|
|
|
995 Args:
|
|
|
996 metabMap (ET.ElementTree): An XML tree representing the metabolic map.
|
|
|
997 rxn_id (str): The ID of the reaction to which the arrow will be applied.
|
|
|
998 color (str): The color of the arrow in hexadecimal format.
|
|
|
999 isNegative (bool): A boolean indicating if the arrow represents a negative value.
|
|
|
1000 width (int): The width of the arrow.
|
|
|
1001
|
|
|
1002 Returns:
|
|
|
1003 None
|
|
|
1004 """
|
|
|
1005 arrow = Arrow(width=width, col=color)
|
|
|
1006 arrow.styleReactionElementsMeanMedian(metabMap, rxn_id, isNegative)
|
|
|
1007 pass
|
|
|
1008
|
|
|
1009 def save_and_convert(metabMap, map_type, key):
|
|
|
1010 """
|
|
|
1011 Save the metabolic map as an SVG file and optionally convert it to PNG and PDF formats.
|
|
|
1012
|
|
|
1013 Args:
|
|
|
1014 metabMap (ET.ElementTree): An XML tree representing the metabolic map.
|
|
|
1015 map_type (str): The type of map ('mean' or 'median').
|
|
|
1016 key (str): The key identifying the specific map.
|
|
|
1017
|
|
|
1018 Returns:
|
|
|
1019 None
|
|
|
1020 """
|
|
|
1021 svgFilePath = utils.FilePath(f"SVG Map {map_type} - {key}", ext=utils.FileFormat.SVG, prefix=ARGS.output_path)
|
|
|
1022 utils.writeSvg(svgFilePath, metabMap)
|
|
|
1023 if ARGS.generate_pdf:
|
|
|
1024 pngPath = utils.FilePath(f"PNG Map {map_type} - {key}", ext=utils.FileFormat.PNG, prefix=ARGS.output_path)
|
|
|
1025 pdfPath = utils.FilePath(f"PDF Map {map_type} - {key}", ext=utils.FileFormat.PDF, prefix=ARGS.output_path)
|
|
|
1026 convert_to_pdf(svgFilePath, pngPath, pdfPath)
|
|
|
1027 if not ARGS.generate_svg:
|
|
|
1028 os.remove(svgFilePath.show())
|
|
|
1029
|
|
|
1030 ############################ MAIN #############################################
|
|
|
1031 def main(args:List[str] = None) -> None:
|
|
|
1032 """
|
|
|
1033 Initializes everything and sets the program in motion based on the fronted input arguments.
|
|
|
1034
|
|
|
1035 Returns:
|
|
|
1036 None
|
|
|
1037
|
|
|
1038 Raises:
|
|
|
1039 sys.exit : if a user-provided custom map is in the wrong format (ET.XMLSyntaxError, ET.XMLSchemaParseError)
|
|
|
1040 """
|
|
|
1041
|
|
|
1042 global ARGS
|
|
|
1043 ARGS = process_args(args)
|
|
|
1044
|
|
|
1045 if ARGS.custom_map == 'None':
|
|
|
1046 ARGS.custom_map = None
|
|
|
1047
|
|
|
1048 if os.path.isdir(ARGS.output_path) == False: os.makedirs(ARGS.output_path)
|
|
|
1049
|
|
|
1050 core_map :ET.ElementTree = ARGS.choice_map.getMap(
|
|
|
1051 ARGS.tool_dir,
|
|
|
1052 utils.FilePath.fromStrPath(ARGS.custom_map) if ARGS.custom_map else None)
|
|
|
1053 # TODO: ^^^ ugly but fine for now, the argument is None if the model isn't custom because no file was given.
|
|
|
1054 # getMap will None-check the customPath and panic when the model IS custom but there's no file (good). A cleaner
|
|
|
1055 # solution can be derived from my comment in FilePath.fromStrPath
|
|
|
1056
|
|
|
1057 ids, class_pat = getClassesAndIdsFromDatasets(ARGS.input_datas_fluxes, ARGS.input_data_fluxes, ARGS.input_class_fluxes, ARGS.names_fluxes)
|
|
|
1058
|
|
|
1059 if(ARGS.choice_map == utils.Model.HMRcore):
|
|
|
1060 temp_map = utils.Model.HMRcore_no_legend
|
|
|
1061 computeEnrichmentMeanMedian(temp_map.getMap(ARGS.tool_dir), class_pat, ids, ARGS.color_map)
|
|
|
1062 elif(ARGS.choice_map == utils.Model.ENGRO2):
|
|
|
1063 temp_map = utils.Model.ENGRO2_no_legend
|
|
|
1064 computeEnrichmentMeanMedian(temp_map.getMap(ARGS.tool_dir), class_pat, ids, ARGS.color_map)
|
|
|
1065 else:
|
|
|
1066 computeEnrichmentMeanMedian(core_map, class_pat, ids, ARGS.color_map)
|
|
|
1067
|
|
|
1068
|
|
|
1069 enrichment_results = computeEnrichment(class_pat, ids)
|
|
|
1070 for i, j, comparisonDict, max_z_score in enrichment_results:
|
|
|
1071 map_copy = copy.deepcopy(core_map)
|
|
|
1072 temp_thingsInCommon(comparisonDict, map_copy, max_z_score, i, j)
|
|
|
1073 createOutputMaps(i, j, map_copy)
|
|
|
1074
|
|
|
1075 if not ERRORS: return
|
|
|
1076 utils.logWarning(
|
|
|
1077 f"The following reaction IDs were mentioned in the dataset but weren't found in the map: {ERRORS}",
|
|
|
1078 ARGS.out_log)
|
|
|
1079
|
|
|
1080 print('Execution succeded')
|
|
|
1081
|
|
|
1082 ###############################################################################
|
|
|
1083 if __name__ == "__main__":
|
|
|
1084 main()
|
|
|
1085
|
