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