Mercurial > repos > bimib > cobraxy
diff COBRAxy/src/flux_to_map.py @ 539:2fb97466e404 draft
Uploaded
| author | francesco_lapi |
|---|---|
| date | Sat, 25 Oct 2025 14:55:13 +0000 |
| parents | |
| children | fcdbc81feb45 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/COBRAxy/src/flux_to_map.py Sat Oct 25 14:55:13 2025 +0000 @@ -0,0 +1,1082 @@ +from __future__ import division +import csv +from enum import Enum +import re +import sys +import numpy as np +import pandas as pd +import itertools as it +import scipy.stats as st +import lxml.etree as ET +import math +import utils.general_utils as utils +from PIL import Image +import os +import copy +import argparse +import pyvips +from PIL import Image +from typing import Tuple, Union, Optional, List, Dict +import matplotlib.pyplot as plt + +ERRORS = [] +########################## argparse ########################################## +ARGS :argparse.Namespace +def process_args(args:List[str] = None) -> argparse.Namespace: + """ + Interfaces the script of a module with its frontend, making the user's choices for various parameters available as values in code. + + Args: + args : Always obtained (in file) from sys.argv + + Returns: + Namespace : An object containing the parsed arguments + """ + parser = argparse.ArgumentParser( + usage = "%(prog)s [options]", + description = "process some value's genes to create a comparison's map.") + + #General: + parser.add_argument( + '-td', '--tool_dir', + type = str, + required = True, + help = 'your tool directory') + + parser.add_argument('-on', '--control', type = str) + parser.add_argument('-ol', '--out_log', help = "Output log") + + #Computation details: + parser.add_argument( + '-co', '--comparison', + type = str, + default = 'manyvsmany', + choices = ['manyvsmany', 'onevsrest', 'onevsmany']) + + parser.add_argument( + '-te' ,'--test', + type = str, + default = 'ks', + choices = ['ks', 'ttest_p', 'ttest_ind', 'wilcoxon', 'mw'], + help = 'Statistical test to use (default: %(default)s)') + + parser.add_argument( + '-pv' ,'--pValue', + type = float, + default = 0.1, + help = 'P-Value threshold (default: %(default)s)') + + parser.add_argument( + '-adj' ,'--adjusted', + type = utils.Bool("adjusted"), default = False, + help = 'Apply the FDR (Benjamini-Hochberg) correction (default: %(default)s)') + + parser.add_argument( + '-fc', '--fChange', + type = float, + default = 1.5, + help = 'Fold-Change threshold (default: %(default)s)') + + parser.add_argument( + '-op', '--option', + type = str, + choices = ['datasets', 'dataset_class'], + help='dataset or dataset and class') + + parser.add_argument( + '-idf', '--input_data_fluxes', + type = str, + help = 'input dataset fluxes') + + parser.add_argument( + '-icf', '--input_class_fluxes', + type = str, + help = 'sample group specification fluxes') + + parser.add_argument( + '-idsf', '--input_datas_fluxes', + type = str, + nargs = '+', + help = 'input datasets fluxes') + + parser.add_argument( + '-naf', '--names_fluxes', + type = str, + nargs = '+', + help = 'input names fluxes') + + #Output: + parser.add_argument( + "-gs", "--generate_svg", + type = utils.Bool("generate_svg"), default = True, + help = "choose whether to generate svg") + + parser.add_argument( + "-gp", "--generate_pdf", + type = utils.Bool("generate_pdf"), default = True, + help = "choose whether to generate pdf") + + parser.add_argument( + '-cm', '--custom_map', + type = str, + help='custom map to use') + + parser.add_argument( + '-mc', '--choice_map', + type = utils.Model, default = utils.Model.HMRcore, + choices = [utils.Model.HMRcore, utils.Model.ENGRO2, utils.Model.Custom]) + + parser.add_argument( + '-colorm', '--color_map', + type = str, + choices = ["jet", "viridis"]) + + parser.add_argument( + '-idop', '--output_path', + type = str, + default='result', + help = 'output path for maps') + + args :argparse.Namespace = parser.parse_args(args) + args.net = True # TODO SICCOME I FLUSSI POSSONO ESSERE ANCHE NEGATIVI SONO SEMPRE CONSIDERATI NETTI + + return args + +############################ dataset input #################################### +def read_dataset(data :str, name :str) -> pd.DataFrame: + """ + Tries to read the dataset from its path (data) as a tsv and turns it into a DataFrame. + + Args: + data : filepath of a dataset (from frontend input params or literals upon calling) + name : name associated with the dataset (from frontend input params or literals upon calling) + + Returns: + pd.DataFrame : dataset in a runtime operable shape + + Raises: + sys.exit : if there's no data (pd.errors.EmptyDataError) or if the dataset has less than 2 columns + """ + try: + dataset = pd.read_csv(data, sep = '\t', header = 0, engine='python') + except pd.errors.EmptyDataError: + sys.exit('Execution aborted: wrong format of ' + name + '\n') + if len(dataset.columns) < 2: + sys.exit('Execution aborted: wrong format of ' + name + '\n') + return dataset + +############################ dataset name ##################################### +def name_dataset(name_data :str, count :int) -> str: + """ + 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. + + Args: + name_data : name associated with the dataset (from frontend input params) + count : counter from 1 to make these names unique (external) + + Returns: + str : the name made unique + """ + if str(name_data) == 'Dataset': + return str(name_data) + '_' + str(count) + else: + return str(name_data) + +############################ map_methods ###################################### +FoldChange = Union[float, int, str] # Union[float, Literal[0, "-INF", "INF"]] +def fold_change(avg1 :float, avg2 :float) -> FoldChange: + """ + Calculates the fold change between two gene expression values. + + Args: + avg1 : average expression value from one dataset avg2 : average expression value from the other dataset + + Returns: + FoldChange : + 0 : when both input values are 0 + "-INF" : when avg1 is 0 + "INF" : when avg2 is 0 + float : for any other combination of values + """ + if avg1 == 0 and avg2 == 0: + return 0 + elif avg1 == 0: + return '-INF' + elif avg2 == 0: + return 'INF' + else: # (threshold_F_C - 1) / (abs(threshold_F_C) + 1) con threshold_F_C > 1 + return (avg1 - avg2) / (abs(avg1) + abs(avg2)) + +def getElementById(reactionId :str, metabMap :ET.ElementTree) -> utils.Result[ET.Element, utils.Result.ResultErr]: + """ + Finds any element in the given map with the given ID. ID uniqueness in an svg file is recommended but + not enforced, if more than one element with the exact ID is found only the first will be returned. + + Args: + reactionId (str): exact ID of the requested element. + metabMap (ET.ElementTree): metabolic map containing the element. + + Returns: + utils.Result[ET.Element, ResultErr]: result of the search, either the first match found or a ResultErr. + """ + return utils.Result.Ok( + f"//*[@id=\"{reactionId}\"]").map( + lambda xPath : metabMap.xpath(xPath)[0]).mapErr( + lambda _ : utils.Result.ResultErr(f"No elements with ID \"{reactionId}\" found in map")) + # ^^^ we shamelessly ignore the contents of the IndexError, it offers nothing to the user. + +def styleMapElement(element :ET.Element, styleStr :str) -> None: + currentStyles :str = element.get("style", "") + if re.search(r";stroke:[^;]+;stroke-width:[^;]+;stroke-dasharray:[^;]+$", currentStyles): + currentStyles = ';'.join(currentStyles.split(';')[:-3]) + + element.set("style", currentStyles + styleStr) + +class ReactionDirection(Enum): + Unknown = "" + Direct = "_F" + Inverse = "_B" + + @classmethod + def fromDir(cls, s :str) -> "ReactionDirection": + # vvv as long as there's so few variants I actually condone the if spam: + if s == ReactionDirection.Direct.value: return ReactionDirection.Direct + if s == ReactionDirection.Inverse.value: return ReactionDirection.Inverse + return ReactionDirection.Unknown + + @classmethod + def fromReactionId(cls, reactionId :str) -> "ReactionDirection": + return ReactionDirection.fromDir(reactionId[-2:]) + +def getArrowBodyElementId(reactionId :str) -> str: + if reactionId.endswith("_RV"): reactionId = reactionId[:-3] #TODO: standardize _RV + elif ReactionDirection.fromReactionId(reactionId) is not ReactionDirection.Unknown: reactionId = reactionId[:-2] + return f"R_{reactionId}" + +def getArrowHeadElementId(reactionId :str) -> Tuple[str, str]: + """ + We attempt extracting the direction information from the provided reaction ID, if unsuccessful we provide the IDs of both directions. + + Args: + reactionId : the provided reaction ID. + + Returns: + Tuple[str, str]: either a single str ID for the correct arrow head followed by an empty string or both options to try. + """ + if reactionId.endswith("_RV"): reactionId = reactionId[:-3] #TODO: standardize _RV + elif ReactionDirection.fromReactionId(reactionId) is not ReactionDirection.Unknown: return reactionId[:-3:-1] + reactionId[:-2], "" + return f"F_{reactionId}", f"B_{reactionId}" + +class ArrowColor(Enum): + """ + Encodes possible arrow colors based on their meaning in the enrichment process. + """ + Invalid = "#BEBEBE" # gray, fold-change under treshold or not significant p-value + Transparent = "#ffffff00" # transparent, to make some arrow segments disappear + UpRegulated = "#ecac68" # red, up-regulated reaction + DownRegulated = "#6495ed" # blue, down-regulated reaction + + UpRegulatedInv = "#FF0000" + # ^^^ different shade of red (actually orange), up-regulated net value for a reversible reaction with + # conflicting enrichment in the two directions. + + DownRegulatedInv = "#0000FF" + # ^^^ different shade of blue (actually purple), down-regulated net value for a reversible reaction with + # conflicting enrichment in the two directions. + + @classmethod + def fromFoldChangeSign(cls, foldChange :float, *, useAltColor = False) -> "ArrowColor": + colors = (cls.DownRegulated, cls.DownRegulatedInv) if foldChange < 0 else (cls.UpRegulated, cls.UpRegulatedInv) + return colors[useAltColor] + + def __str__(self) -> str: return self.value + +class Arrow: + """ + Models the properties of a reaction arrow that change based on enrichment. + """ + MIN_W = 2 + MAX_W = 12 + + def __init__(self, width :int, col: ArrowColor, *, isDashed = False) -> None: + """ + (Private) Initializes an instance of Arrow. + + Args: + width : width of the arrow, ideally to be kept within Arrow.MIN_W and Arrow.MAX_W (not enforced). + col : color of the arrow. + isDashed : whether the arrow should be dashed, meaning the associated pValue resulted not significant. + + Returns: + None : practically, a Arrow instance. + """ + self.w = width + self.col = col + self.dash = isDashed + + def applyTo(self, reactionId :str, metabMap :ET.ElementTree, styleStr :str) -> None: + if getElementById(reactionId, metabMap).map(lambda el : styleMapElement(el, styleStr)).isErr: + ERRORS.append(reactionId) + + def styleReactionElements(self, metabMap :ET.ElementTree, reactionId :str, *, mindReactionDir = True) -> None: + if not mindReactionDir: + return self.applyTo(getArrowBodyElementId(reactionId), metabMap, self.toStyleStr()) + + # Now we style the arrow head(s): + idOpt1, idOpt2 = getArrowHeadElementId(reactionId) + self.applyTo(idOpt1, metabMap, self.toStyleStr(downSizedForTips = True)) + if idOpt2: self.applyTo(idOpt2, metabMap, self.toStyleStr(downSizedForTips = True)) + + def styleReactionElementsMeanMedian(self, metabMap :ET.ElementTree, reactionId :str, isNegative:bool) -> None: + + self.applyTo(getArrowBodyElementId(reactionId), metabMap, self.toStyleStr()) + idOpt1, idOpt2 = getArrowHeadElementId(reactionId) + + if(isNegative): + self.applyTo(idOpt2, metabMap, self.toStyleStr(downSizedForTips = True)) + self.col = ArrowColor.Transparent + self.applyTo(idOpt1, metabMap, self.toStyleStr(downSizedForTips = True)) #trasp + else: + self.applyTo(idOpt1, metabMap, self.toStyleStr(downSizedForTips = True)) + self.col = ArrowColor.Transparent + self.applyTo(idOpt2, metabMap, self.toStyleStr(downSizedForTips = True)) #trasp + + + + def getMapReactionId(self, reactionId :str, mindReactionDir :bool) -> str: + """ + Computes the reaction ID as encoded in the map for a given reaction ID from the dataset. + + Args: + reactionId: the reaction ID, as encoded in the dataset. + mindReactionDir: if True forward (F_) and backward (B_) directions will be encoded in the result. + + Returns: + str : the ID of an arrow's body or tips in the map. + """ + # we assume the reactionIds also don't encode reaction dir if they don't mind it when styling the map. + if not mindReactionDir: return "R_" + reactionId + + #TODO: this is clearly something we need to make consistent in fluxes + return (reactionId[:-3:-1] + reactionId[:-2]) if reactionId[:-2] in ["_F", "_B"] else f"F_{reactionId}" # "Pyr_F" --> "F_Pyr" + + def toStyleStr(self, *, downSizedForTips = False) -> str: + """ + Collapses the styles of this Arrow into a str, ready to be applied as part of the "style" property on an svg element. + + Returns: + str : the styles string. + """ + width = self.w + if downSizedForTips: width *= 0.8 + return f";stroke:{self.col};stroke-width:{width};stroke-dasharray:{'5,5' if self.dash else 'none'}" + +# vvv These constants could be inside the class itself a static properties, but python +# was built by brainless organisms so here we are! +INVALID_ARROW = Arrow(Arrow.MIN_W, ArrowColor.Invalid) +INSIGNIFICANT_ARROW = Arrow(Arrow.MIN_W, ArrowColor.Invalid, isDashed = True) + +def applyFluxesEnrichmentToMap(fluxesEnrichmentRes :Dict[str, Union[Tuple[float, FoldChange], Tuple[float, FoldChange, float, float]]], metabMap :ET.ElementTree, maxNumericZScore :float) -> None: + """ + Applies fluxes enrichment results to the provided metabolic map. + + Args: + fluxesEnrichmentRes : fluxes enrichment results. + metabMap : the metabolic map to edit. + maxNumericZScore : biggest finite z-score value found. + + Side effects: + metabMap : mut + + Returns: + None + """ + for reactionId, values in fluxesEnrichmentRes.items(): + pValue = values[0] + foldChange = values[1] + z_score = values[2] + + if math.isnan(pValue) or (isinstance(foldChange, float) and math.isnan(foldChange)): + continue + + if isinstance(foldChange, str): foldChange = float(foldChange) + if pValue > ARGS.pValue: # pValue above tresh: dashed arrow + INSIGNIFICANT_ARROW.styleReactionElements(metabMap, reactionId) + INSIGNIFICANT_ARROW.styleReactionElements(metabMap, reactionId, mindReactionDir = False) + + continue + + if abs(foldChange) < (ARGS.fChange - 1) / (abs(ARGS.fChange) + 1): + INVALID_ARROW.styleReactionElements(metabMap, reactionId) + INVALID_ARROW.styleReactionElements(metabMap, reactionId, mindReactionDir = False) + + continue + + width = Arrow.MAX_W + if not math.isinf(z_score): + try: + width = min( + max(abs(z_score * Arrow.MAX_W) / maxNumericZScore, Arrow.MIN_W), + Arrow.MAX_W) + + except ZeroDivisionError: pass + # TODO CHECK RV + #if not reactionId.endswith("_RV"): # RV stands for reversible reactions + # Arrow(width, ArrowColor.fromFoldChangeSign(foldChange)).styleReactionElements(metabMap, reactionId) + # continue + + #reactionId = reactionId[:-3] # Remove "_RV" + + inversionScore = (values[3] < 0) + (values[4] < 0) # Compacts the signs of averages into 1 easy to check score + if inversionScore == 2: foldChange *= -1 + # ^^^ Style the inverse direction with the opposite sign netValue + + # If the score is 1 (opposite signs) we use alternative colors vvv + arrow = Arrow(width, ArrowColor.fromFoldChangeSign(foldChange, useAltColor = inversionScore == 1)) + + # vvv These 2 if statements can both be true and can both happen + if ARGS.net: # style arrow head(s): + arrow.styleReactionElements(metabMap, reactionId + ("_B" if inversionScore == 2 else "_F")) + arrow.applyTo(("F_" if inversionScore == 2 else "B_") + reactionId, metabMap, f";stroke:{ArrowColor.Transparent};stroke-width:0;stroke-dasharray:None") + + arrow.styleReactionElements(metabMap, reactionId, mindReactionDir = False) + + +############################ split class ###################################### +def split_class(classes :pd.DataFrame, resolve_rules :Dict[str, List[float]]) -> Dict[str, List[List[float]]]: + """ + Generates a :dict that groups together data from a :DataFrame based on classes the data is related to. + + Args: + classes : a :DataFrame of only string values, containing class information (rows) and keys to query the resolve_rules :dict + resolve_rules : a :dict containing :float data + + Returns: + dict : the dict with data grouped by class + + Side effects: + classes : mut + """ + class_pat :Dict[str, List[List[float]]] = {} + for i in range(len(classes)): + classe :str = classes.iloc[i, 1] + if pd.isnull(classe): continue + + l :List[List[float]] = [] + for j in range(i, len(classes)): + if classes.iloc[j, 1] == classe: + pat_id :str = classes.iloc[j, 0] + tmp = resolve_rules.get(pat_id, None) + if tmp != None: + l.append(tmp) + classes.iloc[j, 1] = None + + if l: + class_pat[classe] = list(map(list, zip(*l))) + continue + + utils.logWarning( + f"Warning: no sample found in class \"{classe}\", the class has been disregarded", ARGS.out_log) + + return class_pat + +############################ conversion ############################################## +#conversion from svg to png +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: + """ + Internal utility to convert an SVG to PNG (forced opaque) to aid in PDF conversion. + + Args: + svg_path : path to SVG file + png_path : path for new PNG file + dpi : dots per inch of the generated PNG + scale : scaling factor for the generated PNG, computed internally when a size is provided + size : final effective width of the generated PNG + + Returns: + None + """ + if size: + image = pyvips.Image.new_from_file(svg_path.show(), dpi=dpi, scale=1) + scale = size / image.width + image = image.resize(scale) + else: + image = pyvips.Image.new_from_file(svg_path.show(), dpi=dpi, scale=scale) + + white_background = pyvips.Image.black(image.width, image.height).new_from_image([255, 255, 255]) + white_background = white_background.affine([scale, 0, 0, scale]) + + if white_background.bands != image.bands: + white_background = white_background.extract_band(0) + + composite_image = white_background.composite2(image, 'over') + composite_image.write_to_file(png_path.show()) + +#funzione unica, lascio fuori i file e li passo in input +#conversion from png to pdf +def convert_png_to_pdf(png_file :utils.FilePath, pdf_file :utils.FilePath) -> None: + """ + Internal utility to convert a PNG to PDF to aid from SVG conversion. + + Args: + png_file : path to PNG file + pdf_file : path to new PDF file + + Returns: + None + """ + image = Image.open(png_file.show()) + image = image.convert("RGB") + image.save(pdf_file.show(), "PDF", resolution=100.0) + +#function called to reduce redundancy in the code +def convert_to_pdf(file_svg :utils.FilePath, file_png :utils.FilePath, file_pdf :utils.FilePath) -> None: + """ + Converts the SVG map at the provided path to PDF. + + Args: + file_svg : path to SVG file + file_png : path to PNG file + file_pdf : path to new PDF file + + Returns: + None + """ + svg_to_png_with_background(file_svg, file_png) + try: + convert_png_to_pdf(file_png, file_pdf) + print(f'PDF file {file_pdf.filePath} successfully generated.') + + except Exception as e: + raise utils.DataErr(file_pdf.show(), f'Error generating PDF file: {e}') + +############################ map ############################################## +def buildOutputPath(dataset1Name :str, dataset2Name = "rest", *, details = "", ext :utils.FileFormat) -> utils.FilePath: + """ + Builds a FilePath instance from the names of confronted datasets ready to point to a location in the + "result/" folder, used by this tool for output files in collections. + + Args: + dataset1Name : _description_ + dataset2Name : _description_. Defaults to "rest". + details : _description_ + ext : _description_ + + Returns: + utils.FilePath : _description_ + """ + # This function returns a util data structure but is extremely specific to this module. + # RAS also uses collections as output and as such might benefit from a method like this, but I'd wait + # TODO: until a third tool with multiple outputs appears before porting this to utils. + return utils.FilePath( + f"{dataset1Name}_vs_{dataset2Name}" + (f" ({details})" if details else ""), + # ^^^ yes this string is built every time even if the form is the same for the same 2 datasets in + # all output files: I don't care, this was never the performance bottleneck of the tool and + # there is no other net gain in saving and re-using the built string. + ext, + prefix = ARGS.output_path) + +FIELD_NOT_AVAILABLE = '/' +def writeToCsv(rows: List[list], fieldNames :List[str], outPath :utils.FilePath) -> None: + fieldsAmt = len(fieldNames) + with open(outPath.show(), "w", newline = "") as fd: + writer = csv.DictWriter(fd, fieldnames = fieldNames, delimiter = '\t') + writer.writeheader() + + for row in rows: + sizeMismatch = fieldsAmt - len(row) + if sizeMismatch > 0: row.extend([FIELD_NOT_AVAILABLE] * sizeMismatch) + writer.writerow({ field : data for field, data in zip(fieldNames, row) }) + +OldEnrichedScores = Dict[str, List[Union[float, FoldChange]]] #TODO: try to use Tuple whenever possible +def writeTabularResult(enrichedScores : OldEnrichedScores, outPath :utils.FilePath) -> None: + fieldNames = ["ids", "P_Value", "fold change", "z-score"] + fieldNames.extend(["average_1", "average_2"]) + + writeToCsv([ [reactId] + values for reactId, values in enrichedScores.items() ], fieldNames, outPath) + +def temp_thingsInCommon(tmp :Dict[str, List[Union[float, FoldChange]]], core_map :ET.ElementTree, max_z_score :float, dataset1Name :str, dataset2Name = "rest") -> None: + # this function compiles the things always in common between comparison modes after enrichment. + # TODO: organize, name better. + writeTabularResult(tmp, buildOutputPath(dataset1Name, dataset2Name, details = "Tabular Result", ext = utils.FileFormat.TSV)) + for reactId, enrichData in tmp.items(): tmp[reactId] = tuple(enrichData) + applyFluxesEnrichmentToMap(tmp, core_map, max_z_score) + +def computePValue(dataset1Data: List[float], dataset2Data: List[float]) -> Tuple[float, float]: + """ + Computes the statistical significance score (P-value) of the comparison between coherent data + from two datasets. The data is supposed to, in both datasets: + - be related to the same reaction ID; + - be ordered by sample, such that the item at position i in both lists is related to the + same sample or cell line. + + Args: + dataset1Data : data from the 1st dataset. + dataset2Data : data from the 2nd dataset. + + Returns: + tuple: (P-value, Z-score) + - P-value from the selected test on the provided data. + - Z-score of the difference between means of the two datasets. + """ + + match ARGS.test: + case "ks": + # Perform Kolmogorov-Smirnov test + _, p_value = st.ks_2samp(dataset1Data, dataset2Data) + case "ttest_p": + # Datasets should have same size + if len(dataset1Data) != len(dataset2Data): + raise ValueError("Datasets must have the same size for paired t-test.") + # Perform t-test for paired samples + _, p_value = st.ttest_rel(dataset1Data, dataset2Data) + case "ttest_ind": + # Perform t-test for independent samples + _, p_value = st.ttest_ind(dataset1Data, dataset2Data) + case "wilcoxon": + # Datasets should have same size + if len(dataset1Data) != len(dataset2Data): + raise ValueError("Datasets must have the same size for Wilcoxon signed-rank test.") + # Perform Wilcoxon signed-rank test + np.random.seed(42) # Ensure reproducibility since zsplit method is used + _, p_value = st.wilcoxon(dataset1Data, dataset2Data, zero_method="zsplit") + case "mw": + # Perform Mann-Whitney U test + _, p_value = st.mannwhitneyu(dataset1Data, dataset2Data) + + # Calculate means and standard deviations + mean1 = np.nanmean(dataset1Data) + mean2 = np.nanmean(dataset2Data) + std1 = np.nanstd(dataset1Data, ddof=1) + std2 = np.nanstd(dataset2Data, ddof=1) + + n1 = len(dataset1Data) + n2 = len(dataset2Data) + + # Calculate Z-score + z_score = (mean1 - mean2) / np.sqrt((std1**2 / n1) + (std2**2 / n2)) + + return p_value, z_score + +def compareDatasetPair(dataset1Data :List[List[float]], dataset2Data :List[List[float]], ids :List[str]) -> Tuple[Dict[str, List[Union[float, FoldChange]]], float]: + #TODO: the following code still suffers from "dumbvarnames-osis" + comparisonResult :Dict[str, List[Union[float, FoldChange]]] = {} + count = 0 + max_z_score = 0 + for l1, l2 in zip(dataset1Data, dataset2Data): + reactId = ids[count] + count += 1 + if not reactId: continue # we skip ids that have already been processed + + try: + p_value, z_score = computePValue(l1, l2) + avg1 = sum(l1) / len(l1) + avg2 = sum(l2) / len(l2) + f_c = fold_change(avg1, avg2) + if np.isfinite(z_score) and max_z_score < abs(z_score): max_z_score = abs(z_score) + + comparisonResult[reactId] = [float(p_value), f_c, z_score, avg1, avg2] + except (TypeError, ZeroDivisionError): continue + + # Apply multiple testing correction if set by the user + if ARGS.adjusted: + + # Retrieve the p-values from the comparisonResult dictionary, they have to be different from NaN + validPValues = [(reactId, result[0]) for reactId, result in comparisonResult.items() if not np.isnan(result[0])] + + if not validPValues: + return comparisonResult, max_z_score + + # Unpack the valid p-values + reactIds, pValues = zip(*validPValues) + # Adjust the p-values using the Benjamini-Hochberg method + adjustedPValues = st.false_discovery_control(pValues) + # Update the comparisonResult dictionary with the adjusted p-values + for reactId , adjustedPValue in zip(reactIds, adjustedPValues): + comparisonResult[reactId][0] = adjustedPValue + + return comparisonResult, max_z_score + +def computeEnrichment(class_pat :Dict[str, List[List[float]]], ids :List[str]) -> List[Tuple[str, str, dict, float]]: + """ + Compares clustered data based on a given comparison mode and applies enrichment-based styling on the + provided metabolic map. + + Args: + class_pat : the clustered data. + ids : ids for data association. + + + Returns: + List[Tuple[str, str, dict, float]]: List of tuples with pairs of dataset names, comparison dictionary, and max z-score. + + Raises: + sys.exit : if there are less than 2 classes for comparison + + """ + class_pat = { k.strip() : v for k, v in class_pat.items() } + #TODO: simplfy this stuff vvv and stop using sys.exit (raise the correct utils error) + if (not class_pat) or (len(class_pat.keys()) < 2): sys.exit('Execution aborted: classes provided for comparisons are less than two\n') + + enrichment_results = [] + + + if ARGS.comparison == "manyvsmany": + for i, j in it.combinations(class_pat.keys(), 2): + comparisonDict, max_z_score = compareDatasetPair(class_pat.get(i), class_pat.get(j), ids) + enrichment_results.append((i, j, comparisonDict, max_z_score)) + + elif ARGS.comparison == "onevsrest": + for single_cluster in class_pat.keys(): + rest = [item for k, v in class_pat.items() if k != single_cluster for item in v] + + comparisonDict, max_z_score = compareDatasetPair(class_pat.get(single_cluster), rest, ids) + enrichment_results.append((single_cluster, "rest", comparisonDict, max_z_score)) + + #elif ARGS.comparison == "onevsmany": + # controlItems = class_pat.get(ARGS.control) + # for otherDataset in class_pat.keys(): + # if otherDataset == ARGS.control: + # continue + # comparisonDict, max_z_score = compareDatasetPair(controlItems, class_pat.get(otherDataset), ids) + # enrichment_results.append((ARGS.control, otherDataset, comparisonDict, max_z_score)) + elif ARGS.comparison == "onevsmany": + controlItems = class_pat.get(ARGS.control) + for otherDataset in class_pat.keys(): + if otherDataset == ARGS.control: + continue + comparisonDict, max_z_score = compareDatasetPair(class_pat.get(otherDataset),controlItems, ids) + enrichment_results.append(( otherDataset,ARGS.control, comparisonDict, max_z_score)) + + return enrichment_results + +def createOutputMaps(dataset1Name :str, dataset2Name :str, core_map :ET.ElementTree) -> None: + svgFilePath = buildOutputPath(dataset1Name, dataset2Name, details="SVG Map", ext=utils.FileFormat.SVG) + utils.writeSvg(svgFilePath, core_map) + + if ARGS.generate_pdf: + pngPath = buildOutputPath(dataset1Name, dataset2Name, details="PNG Map", ext=utils.FileFormat.PNG) + pdfPath = buildOutputPath(dataset1Name, dataset2Name, details="PDF Map", ext=utils.FileFormat.PDF) + convert_to_pdf(svgFilePath, pngPath, pdfPath) + + if not ARGS.generate_svg: + os.remove(svgFilePath.show()) + +ClassPat = Dict[str, List[List[float]]] +def getClassesAndIdsFromDatasets(datasetsPaths :List[str], datasetPath :str, classPath :str, names :List[str]) -> Tuple[List[str], ClassPat]: + # TODO: I suggest creating dicts with ids as keys instead of keeping class_pat and ids separate, + # for the sake of everyone's sanity. + class_pat :ClassPat = {} + if ARGS.option == 'datasets': + num = 1 #TODO: the dataset naming function could be a generator + for path, name in zip(datasetsPaths, names): + name = name_dataset(name, num) + resolve_rules_float, ids = getDatasetValues(path, name) + if resolve_rules_float != None: + class_pat[name] = list(map(list, zip(*resolve_rules_float.values()))) + + num += 1 + + elif ARGS.option == "dataset_class": + classes = read_dataset(classPath, "class") + classes = classes.astype(str) + resolve_rules_float, ids = getDatasetValues(datasetPath, "Dataset Class (not actual name)") + #check if classes have match on ids + if not all(classes.iloc[:, 0].isin(ids)): + utils.logWarning( + "No match between classes and sample IDs", ARGS.out_log) + if resolve_rules_float != None: class_pat = split_class(classes, resolve_rules_float) + + return ids, class_pat + #^^^ 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) + +#TODO: create these damn args as FilePath objects +def getDatasetValues(datasetPath :str, datasetName :str) -> Tuple[ClassPat, List[str]]: + """ + Opens the dataset at the given path and extracts the values (expected nullable numerics) and the IDs. + + Args: + datasetPath : path to the dataset + datasetName (str): dataset name, used in error reporting + + Returns: + Tuple[ClassPat, List[str]]: values and IDs extracted from the dataset + """ + dataset = read_dataset(datasetPath, datasetName) + + # Ensure the first column is treated as the reaction name + dataset = dataset.set_index(dataset.columns[0]) + + # Check if required reactions exist in the dataset + required_reactions = ['EX_lac__L_e', 'EX_glc__D_e', 'EX_gln__L_e', 'EX_glu__L_e'] + missing_reactions = [reaction for reaction in required_reactions if reaction not in dataset.index] + + if missing_reactions: + sys.exit(f'Execution aborted: Missing required reactions {missing_reactions} in {datasetName}\n') + + # Calculate new rows using safe division + lact_glc = np.divide( + np.clip(dataset.loc['EX_lac__L_e'].to_numpy(), a_min=0, a_max=None), + np.clip(dataset.loc['EX_glc__D_e'].to_numpy(), a_min=None, a_max=0), + out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan), # Prepara un array con NaN come output di default + where=dataset.loc['EX_glc__D_e'].to_numpy() != 0 # Condizione per evitare la divisione per zero + ) + lact_gln = np.divide( + np.clip(dataset.loc['EX_lac__L_e'].to_numpy(), a_min=0, a_max=None), + np.clip(dataset.loc['EX_gln__L_e'].to_numpy(), a_min=None, a_max=0), + out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan), + where=dataset.loc['EX_gln__L_e'].to_numpy() != 0 + ) + lact_o2 = np.divide( + np.clip(dataset.loc['EX_lac__L_e'].to_numpy(), a_min=0, a_max=None), + np.clip(dataset.loc['EX_o2_e'].to_numpy(), a_min=None, a_max=0), + out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan), + where=dataset.loc['EX_o2_e'].to_numpy() != 0 + ) + glu_gln = np.divide( + dataset.loc['EX_glu__L_e'].to_numpy(), + np.clip(dataset.loc['EX_gln__L_e'].to_numpy(), a_min=None, a_max=0), + out=np.full_like(dataset.loc['EX_lac__L_e'].to_numpy(), np.nan), + where=dataset.loc['EX_gln__L_e'].to_numpy() != 0 + ) + + + values = {'lact_glc': lact_glc, 'lact_gln': lact_gln, 'lact_o2': lact_o2, 'glu_gln': glu_gln} + + # Sostituzione di inf e NaN con 0 se necessario + for key in values: + values[key] = np.nan_to_num(values[key], nan=0.0, posinf=0.0, neginf=0.0) + + # Creazione delle nuove righe da aggiungere al dataset + new_rows = pd.DataFrame({ + dataset.index.name: ['LactGlc', 'LactGln', 'LactO2', 'GluGln'], + **{col: [values['lact_glc'][i], values['lact_gln'][i], values['lact_o2'][i], values['glu_gln'][i]] + for i, col in enumerate(dataset.columns)} + }) + + #print(new_rows) + + # Ritorna il dataset originale con le nuove righe + dataset.reset_index(inplace=True) + dataset = pd.concat([dataset, new_rows], ignore_index=True) + + IDs = pd.Series.tolist(dataset.iloc[:, 0].astype(str)) + + dataset = dataset.drop(dataset.columns[0], axis = "columns").to_dict("list") + return { id : list(map(utils.Float("Dataset values, not an argument"), values)) for id, values in dataset.items() }, IDs + +def rgb_to_hex(rgb): + """ + Convert RGB values (0-1 range) to hexadecimal color format. + + Args: + rgb (numpy.ndarray): An array of RGB color components (in the range [0, 1]). + + Returns: + str: The color in hexadecimal format (e.g., '#ff0000' for red). + """ + # Convert RGB values (0-1 range) to hexadecimal format + rgb = (np.array(rgb) * 255).astype(int) + return '#{:02x}{:02x}{:02x}'.format(rgb[0], rgb[1], rgb[2]) + +def save_colormap_image(min_value: float, max_value: float, path: utils.FilePath, colorMap:str="viridis"): + """ + Create and save an image of the colormap showing the gradient and its range. + + Args: + min_value (float): The minimum value of the colormap range. + max_value (float): The maximum value of the colormap range. + filename (str): The filename for saving the image. + """ + + # Create a colormap using matplotlib + cmap = plt.get_cmap(colorMap) + + # Create a figure and axis + fig, ax = plt.subplots(figsize=(6, 1)) + fig.subplots_adjust(bottom=0.5) + + # Create a gradient image + gradient = np.linspace(0, 1, 256) + gradient = np.vstack((gradient, gradient)) + + # Add min and max value annotations + ax.text(0, 0.5, f'{np.round(min_value, 3)}', va='center', ha='right', transform=ax.transAxes, fontsize=12, color='black') + ax.text(1, 0.5, f'{np.round(max_value, 3)}', va='center', ha='left', transform=ax.transAxes, fontsize=12, color='black') + + + # Display the gradient image + ax.imshow(gradient, aspect='auto', cmap=cmap) + ax.set_axis_off() + + # Save the image + plt.savefig(path.show(), bbox_inches='tight', pad_inches=0) + plt.close() + pass + +def min_nonzero_abs(arr): + # Flatten the array and filter out zeros, then find the minimum of the remaining values + non_zero_elements = np.abs(arr)[np.abs(arr) > 0] + return np.min(non_zero_elements) if non_zero_elements.size > 0 else None + +def computeEnrichmentMeanMedian(metabMap: ET.ElementTree, class_pat: Dict[str, List[List[float]]], ids: List[str], colormap:str) -> None: + """ + Compute and visualize the metabolic map based on mean and median of the input fluxes. + The fluxes are normalised across classes/datasets and visualised using the given colormap. + + Args: + metabMap (ET.ElementTree): An XML tree representing the metabolic map. + class_pat (Dict[str, List[List[float]]]): A dictionary where keys are class names and values are lists of enrichment values. + ids (List[str]): A list of reaction IDs to be used for coloring arrows. + + Returns: + None + """ + # Create copies only if they are needed + metabMap_mean = copy.deepcopy(metabMap) + metabMap_median = copy.deepcopy(metabMap) + + # Compute medians and means + medians = {key: np.round(np.nanmedian(np.array(value), axis=1), 6) for key, value in class_pat.items()} + means = {key: np.round(np.nanmean(np.array(value), axis=1),6) for key, value in class_pat.items()} + + # Normalize medians and means + max_flux_medians = max(np.max(np.abs(arr)) for arr in medians.values()) + max_flux_means = max(np.max(np.abs(arr)) for arr in means.values()) + + min_flux_medians = min(min_nonzero_abs(arr) for arr in medians.values()) + min_flux_means = min(min_nonzero_abs(arr) for arr in means.values()) + + medians = {key: median/max_flux_medians for key, median in medians.items()} + means = {key: mean/max_flux_means for key, mean in means.items()} + + save_colormap_image(min_flux_medians, max_flux_medians, utils.FilePath("Color map median", ext=utils.FileFormat.PNG, prefix=ARGS.output_path), colormap) + save_colormap_image(min_flux_means, max_flux_means, utils.FilePath("Color map mean", ext=utils.FileFormat.PNG, prefix=ARGS.output_path), colormap) + + cmap = plt.get_cmap(colormap) + + min_width = 2.0 # Minimum arrow width + max_width = 15.0 # Maximum arrow width + + for key in class_pat: + # Create color mappings for median and mean + colors_median = { + rxn_id: rgb_to_hex(cmap(abs(medians[key][i]))) if medians[key][i] != 0 else '#bebebe' #grey blocked + for i, rxn_id in enumerate(ids) + } + + colors_mean = { + rxn_id: rgb_to_hex(cmap(abs(means[key][i]))) if means[key][i] != 0 else '#bebebe' #grey blocked + for i, rxn_id in enumerate(ids) + } + + for i, rxn_id in enumerate(ids): + # Calculate arrow width for median + width_median = np.interp(abs(medians[key][i]), [0, 1], [min_width, max_width]) + isNegative = medians[key][i] < 0 + apply_arrow(metabMap_median, rxn_id, colors_median[rxn_id], isNegative, width_median) + + # Calculate arrow width for mean + width_mean = np.interp(abs(means[key][i]), [0, 1], [min_width, max_width]) + isNegative = means[key][i] < 0 + apply_arrow(metabMap_mean, rxn_id, colors_mean[rxn_id], isNegative, width_mean) + + # Save and convert the SVG files + save_and_convert(metabMap_mean, "mean", key) + save_and_convert(metabMap_median, "median", key) + +def apply_arrow(metabMap, rxn_id, color, isNegative, width=5): + """ + Apply an arrow to a specific reaction in the metabolic map with a given color. + + Args: + metabMap (ET.ElementTree): An XML tree representing the metabolic map. + rxn_id (str): The ID of the reaction to which the arrow will be applied. + color (str): The color of the arrow in hexadecimal format. + isNegative (bool): A boolean indicating if the arrow represents a negative value. + width (int): The width of the arrow. + + Returns: + None + """ + arrow = Arrow(width=width, col=color) + arrow.styleReactionElementsMeanMedian(metabMap, rxn_id, isNegative) + pass + +def save_and_convert(metabMap, map_type, key): + """ + Save the metabolic map as an SVG file and optionally convert it to PNG and PDF formats. + + Args: + metabMap (ET.ElementTree): An XML tree representing the metabolic map. + map_type (str): The type of map ('mean' or 'median'). + key (str): The key identifying the specific map. + + Returns: + None + """ + svgFilePath = utils.FilePath(f"SVG Map {map_type} - {key}", ext=utils.FileFormat.SVG, prefix=ARGS.output_path) + utils.writeSvg(svgFilePath, metabMap) + if ARGS.generate_pdf: + pngPath = utils.FilePath(f"PNG Map {map_type} - {key}", ext=utils.FileFormat.PNG, prefix=ARGS.output_path) + pdfPath = utils.FilePath(f"PDF Map {map_type} - {key}", ext=utils.FileFormat.PDF, prefix=ARGS.output_path) + convert_to_pdf(svgFilePath, pngPath, pdfPath) + if not ARGS.generate_svg: + os.remove(svgFilePath.show()) + +############################ MAIN ############################################# +def main(args:List[str] = None) -> None: + """ + Initializes everything and sets the program in motion based on the fronted input arguments. + + Returns: + None + + Raises: + sys.exit : if a user-provided custom map is in the wrong format (ET.XMLSyntaxError, ET.XMLSchemaParseError) + """ + + global ARGS + ARGS = process_args(args) + + if ARGS.custom_map == 'None': + ARGS.custom_map = None + + if os.path.isdir(ARGS.output_path) == False: os.makedirs(ARGS.output_path) + + core_map :ET.ElementTree = ARGS.choice_map.getMap( + ARGS.tool_dir, + utils.FilePath.fromStrPath(ARGS.custom_map) if ARGS.custom_map else None) + # TODO: ^^^ ugly but fine for now, the argument is None if the model isn't custom because no file was given. + # getMap will None-check the customPath and panic when the model IS custom but there's no file (good). A cleaner + # solution can be derived from my comment in FilePath.fromStrPath + + ids, class_pat = getClassesAndIdsFromDatasets(ARGS.input_datas_fluxes, ARGS.input_data_fluxes, ARGS.input_class_fluxes, ARGS.names_fluxes) + + if(ARGS.choice_map == utils.Model.HMRcore): + temp_map = utils.Model.HMRcore_no_legend + computeEnrichmentMeanMedian(temp_map.getMap(ARGS.tool_dir), class_pat, ids, ARGS.color_map) + elif(ARGS.choice_map == utils.Model.ENGRO2): + temp_map = utils.Model.ENGRO2_no_legend + computeEnrichmentMeanMedian(temp_map.getMap(ARGS.tool_dir), class_pat, ids, ARGS.color_map) + else: + computeEnrichmentMeanMedian(core_map, class_pat, ids, ARGS.color_map) + + + enrichment_results = computeEnrichment(class_pat, ids) + for i, j, comparisonDict, max_z_score in enrichment_results: + map_copy = copy.deepcopy(core_map) + temp_thingsInCommon(comparisonDict, map_copy, max_z_score, i, j) + createOutputMaps(i, j, map_copy) + + if not ERRORS: return + utils.logWarning( + f"The following reaction IDs were mentioned in the dataset but weren't found in the map: {ERRORS}", + ARGS.out_log) + + print('Execution succeded') + +############################################################################### +if __name__ == "__main__": + main() +