# HG changeset patch # User luca_milaz # Date 1726651265 0 # Node ID 95cda56f01a938987a43cfad810d73c14df4a861 # Parent 985cd6f42be4dca0cc6513c58edaa93eaf259bf6 Uploaded diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/GSOC project submission.html --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/marea_2/GSOC project submission.html Wed Sep 18 09:21:05 2024 +0000 @@ -0,0 +1,65 @@ + + + + + + Google Summer of Code 2024 - COBRAxy: COBRA and MaREA4Galaxy + + + +
+

Google Summer of Code 2024

+

COBRAxy: COBRA and MaREA4Galaxy

+

National Resource for Network Biology (NRNB)

+

Mentors:

+ +

Contributor:

+ + +

Project Description

+

+ The project focused on developing an advanced Galaxy tool that enhances the data mapping capabilities of MaREA4Galaxy. The extension of this framework includes the analysis of fluxomics data, starting from a metabolic model and progressing to the representation of up-regulated fluxes on a metabolic map. This tool enables users to perform constraint-based enrichment analysis of metabolic pathways. +

+

The primary goals of the project were:

+ + +

What I Did

+ + +

Current State and Future Extensions

+

+ Currently, the updated MaREA4Galaxy tool allows users to perform constraint-based enrichment analysis of metabolic pathways using RNA-seq profiles by simulating fluxomics. Additionally, users can compare different sub-populations identified by the clustering tool. The architecture minimizes computational costs by handling cell-specific models through a set of bounds, without storing complete COBRA models, which would contain a large amount of redundant information. +

+

+ The implementation of the "Metabolic Flux Enrichment Analysis" tool did not leave enough time to extend the clustering module to new algorithms such as HDBSCAN, Leiden, and Louvain. This is a potential future extension to consider. Moreover, implementing a more advanced clustering grid search could further optimize clustering results. +

+ +

About the Code

+

+ I worked on the Mercurial repository of MaREA4Galaxy, where this document is stored. I committed all my changes, as shown by the repository history, though without using any Git-like merge operations due to the limitations of the Mercurial interface. +

+ +

Conclusions

+

+ Over the past years, I have focused on biology-related subjects, particularly metabolic fluxes and other omics data such as gene expression datasets. Through this project, I was able to apply the knowledge I have gained in constraint-based modeling, flux sampling, and omics enrichment analysis by expanding the MaREA4Galaxy tool. This experience not only enhanced my programming skills but also deepened my understanding of the real needs of biologists when working with such omics data. +

+
+ + + diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/README.md --- a/marea_2/README.md Wed Sep 18 09:13:24 2024 +0000 +++ b/marea_2/README.md Wed Sep 18 09:21:05 2024 +0000 @@ -1,6 +1,6 @@ -# Official repository for the MaREA 2.0 toolset -> MaREA 2.0 (Metabolic Reaction Enrichment Analysis) is a user-friendly tool that allows a user to characterize and to graphically compare groups of samples with different transcriptional regulation of metabolism, -as estimated from cross-sectional RNA-seq data and flux simulatio. The tool is available as plug-in for the widely-used Galaxy platform for comparative genomics and bioinformatics analyses. +# Official repository for the MaREA toolset +> MaREA (Metabolic Reaction Enrichment Analysis) is a user-friendly tool that allows a user to characterize and to graphically compare groups of samples with different transcriptional regulation of metabolism, +as estimated from cross-sectional RNA-seq data. The tool is available as plug-in for the widely-used Galaxy platform for comparative genomics and bioinformatics analyses. ## Useful links: - Paper: https://pubmed.ncbi.nlm.nih.gov/32373287/ diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/custom_data_generator.py diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/custom_data_generator.xml diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/flux_simulation.py diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/flux_simulation.xml diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/flux_to_map.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/marea_2/flux_to_map.py Wed Sep 18 09:21:05 2024 +0000 @@ -0,0 +1,1055 @@ +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, ImageDraw, ImageFont +from typing import Tuple, Union, Optional, List, Dict +import matplotlib.pyplot as plt + +ERRORS = [] +########################## argparse ########################################## +ARGS :argparse.Namespace +def process_args() -> 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 = '1vs1', + choices = ['manyvsmany', 'onevsrest', 'onevsmany']) + + parser.add_argument( + '-pv' ,'--pValue', + type = float, + default = 0.1, + help = 'P-Value threshold (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"]) + + args :argparse.Namespace = parser.parse_args() + args.net = True + + 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 fix_style(l :str, col :Optional[str], width :str, dash :str) -> str: + """ + 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. + + Args: + l : current style string of an SVG element + col : new value for the "stroke" style property + width : new value for the "stroke-width" style property + dash : new value for the "stroke-dasharray" style property + + Returns: + str : the fixed style string + """ + tmp = l.split(';') + flag_col = False + flag_width = False + flag_dash = False + for i in range(len(tmp)): + if tmp[i].startswith('stroke:'): + tmp[i] = 'stroke:' + col + flag_col = True + if tmp[i].startswith('stroke-width:'): + tmp[i] = 'stroke-width:' + width + flag_width = True + if tmp[i].startswith('stroke-dasharray:'): + tmp[i] = 'stroke-dasharray:' + dash + flag_dash = True + if not flag_col: + tmp.append('stroke:' + col) + if not flag_width: + tmp.append('stroke-width:' + width) + if not flag_dash: + tmp.append('stroke-dasharray:' + dash) + return ';'.join(tmp) + +# The type of d values is collapsed, losing precision, because the dict containst lists instead of tuples, please fix! +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: + """ + Edits the selected SVG map based on the p-value and fold change data (d) and some significance thresholds also passed as inputs. + + Args: + d : dictionary mapping a p-value and a fold-change value (values) to each reaction ID as encoded in the SVG map (keys) + core_map : SVG map to modify + threshold_P_V : threshold for a p-value to be considered significant + threshold_F_C : threshold for a fold change value to be considered significant + max_z_score : highest z-score (absolute value) + + Returns: + ET.ElementTree : the modified core_map + + Side effects: + core_map : mut + """ + maxT = 12 + minT = 2 + grey = '#BEBEBE' + blue = '#6495ed' + red = '#ecac68' + for el in core_map.iter(): + el_id = str(el.get('id')) + if el_id.startswith('R_'): + tmp = d.get(el_id[2:]) + if tmp != None: + p_val :float = tmp[0] + f_c = tmp[1] + z_score = tmp[2] + if p_val < threshold_P_V: + if not isinstance(f_c, str): + if abs(f_c) < ((threshold_F_C - 1) / (abs(threshold_F_C) + 1)): # + col = grey + width = str(minT) + else: + if f_c < 0: + col = blue + elif f_c > 0: + col = red + width = str(max((abs(z_score) * maxT) / max_z_score, minT)) + else: + if f_c == '-INF': + col = blue + elif f_c == 'INF': + col = red + width = str(maxT) + dash = 'none' + else: + dash = '5,5' + col = grey + width = str(minT) + el.set('style', fix_style(el.get('style', ""), col, width, dash)) + return core_map + +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 + Transparent = "#ffffff00" # white, not significant p-value + 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 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(foldChange): + try: + width = max(abs(z_score * Arrow.MAX_W) / maxNumericZScore, Arrow.MIN_W) + + except ZeroDivisionError: pass + + #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 = "result") + +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"] + 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 a Kolmogorov-Smirnov test on the provided data. + - Z-score of the difference between means of the two datasets. + """ + # Perform Kolmogorov-Smirnov test + ks_statistic, p_value = st.ks_2samp(dataset1Data, dataset2Data) + + # Calculate means and standard deviations + mean1 = np.mean(dataset1Data) + mean2 = np.mean(dataset2Data) + std1 = np.std(dataset1Data, ddof=1) + std2 = np.std(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" + tmp :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) + avg = fold_change(avg1, avg2) + if not isinstance(z_score, str) and max_z_score < abs(z_score): max_z_score = abs(z_score) + tmp[reactId] = [float(p_value), avg, z_score, avg1, avg2] + except (TypeError, ZeroDivisionError): continue + + return tmp, max_z_score + +def computeEnrichment(metabMap :ET.ElementTree, class_pat :Dict[str, List[List[float]]], ids :List[str]) -> None: + """ + Compares clustered data based on a given comparison mode and applies enrichment-based styling on the + provided metabolic map. + + Args: + metabMap : SVG map to modify. + class_pat : the clustered data. + ids : ids for data association. + + + Returns: + None + + Raises: + sys.exit : if there are less than 2 classes for comparison + + Side effects: + metabMap : mut + ids : mut + """ + 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') + + if ARGS.comparison == "manyvsmany": + for i, j in it.combinations(class_pat.keys(), 2): + #TODO: these 2 functions are always called in pair and in this order and need common data, + # some clever refactoring would be appreciated. + comparisonDict, max_z_score = compareDatasetPair(class_pat.get(i), class_pat.get(j), ids) + temp_thingsInCommon(comparisonDict, metabMap, max_z_score, i, j) + + elif ARGS.comparison == "onevsrest": + for single_cluster in class_pat.keys(): + t :List[List[List[float]]] = [] + for k in class_pat.keys(): + if k != single_cluster: + t.append(class_pat.get(k)) + + rest :List[List[float]] = [] + for i in t: + rest = rest + i + + comparisonDict, max_z_score = compareDatasetPair(class_pat.get(single_cluster), rest, ids) + temp_thingsInCommon(comparisonDict, metabMap, max_z_score, single_cluster) + + 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) + temp_thingsInCommon(comparisonDict, metabMap, max_z_score, ARGS.control, otherDataset) + +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)") + 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) + 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.median(np.array(value), axis=1), 6) for key, value in class_pat.items()} + means = {key: np.round(np.mean(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="result"), colormap) + save_colormap_image(min_flux_means, max_flux_means, utils.FilePath("Color map mean", ext=utils.FileFormat.PNG, prefix="result"), colormap) + + cmap = plt.get_cmap(colormap) + + 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): + isNegative = medians[key][i] < 0 + + # Apply median arrows + apply_arrow(metabMap_median, rxn_id, colors_median[rxn_id], isNegative) + + isNegative = means[key][i] < 0 + # Apply mean arrows + apply_arrow(metabMap_mean, rxn_id, colors_mean[rxn_id], isNegative) + + # 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): + """ + 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. + + Returns: + None + """ + arrow = Arrow(width=5, 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="result") + utils.writeSvg(svgFilePath, metabMap) + if ARGS.generate_pdf: + pngPath = utils.FilePath(f"PNG Map {map_type} - {key}", ext=utils.FileFormat.PNG, prefix="result") + pdfPath = utils.FilePath(f"PDF Map {map_type} - {key}", ext=utils.FileFormat.PDF, prefix="result") + convert_to_pdf(svgFilePath, pngPath, pdfPath) + if not ARGS.generate_svg: + os.remove(svgFilePath.show()) + + + + +############################ MAIN ############################################# +def main() -> 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() + + if os.path.isdir('result') == False: os.makedirs('result') + + 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) + + + computeEnrichment(core_map, class_pat, ids) + + # create output files: TODO: this is the same comparison happening in "maps", find a better way to organize this + if ARGS.comparison == "manyvsmany": + for i, j in it.combinations(class_pat.keys(), 2): createOutputMaps(i, j, core_map) + return + + if ARGS.comparison == "onevsrest": + for single_cluster in class_pat.keys(): createOutputMaps(single_cluster, "rest", core_map) + return + + for otherDataset in class_pat.keys(): + if otherDataset != ARGS.control: createOutputMaps(i, j, core_map) + + 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() \ No newline at end of file diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/flux_to_map.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/marea_2/flux_to_map.xml Wed Sep 18 09:21:05 2024 +0000 @@ -0,0 +1,247 @@ + + + marea_macros.xml + + + + numpy + pandas + seaborn + scipy + svglib + pyvips + cairosvg + cobra + lxml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/marea.xml --- a/marea_2/marea.xml Wed Sep 18 09:13:24 2024 +0000 +++ b/marea_2/marea.xml Wed Sep 18 09:21:05 2024 +0000 @@ -1,4 +1,4 @@ - + marea_macros.xml diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/marea_macros.xml --- a/marea_2/marea_macros.xml Wed Sep 18 09:13:24 2024 +0000 +++ b/marea_2/marea_macros.xml Wed Sep 18 09:21:05 2024 +0000 @@ -140,29 +140,6 @@ url = {https://doi.org/10.1016/j.csbj.2020.04.008}, } - - @article{di2022integrate, - title={INTEGRATE: Model-based multi-omics data integration to characterize multi-level metabolic regulation}, - author={Di Filippo, Marzia and Pescini, Dario and Galuzzi, Bruno Giovanni and Bonanomi, Marcella and Gaglio, Daniela and Mangano, Eleonora and Consolandi, Clarissa and Alberghina, Lilia and Vanoni, Marco and Damiani, Chiara}, - journal={PLoS computational biology}, - volume={18}, - number={2}, - pages={e1009337}, - year={2022}, - publisher={Public Library of Science San Francisco, CA USA} - } - - - @article{galuzzi2024adjusting, - title={Adjusting for false discoveries in constraint-based differential metabolic flux analysis}, - author={Galuzzi, Bruno G and Milazzo, Luca and Damiani, Chiara}, - journal={Journal of Biomedical Informatics}, - volume={150}, - pages={104597}, - year={2024}, - publisher={Elsevier} - } - diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/ras_generator.xml --- a/marea_2/ras_generator.xml Wed Sep 18 09:13:24 2024 +0000 +++ b/marea_2/ras_generator.xml Wed Sep 18 09:21:05 2024 +0000 @@ -1,4 +1,4 @@ - + - Reaction Activity Scores computation marea_macros.xml diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/ras_to_bounds.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/marea_2/ras_to_bounds.py Wed Sep 18 09:21:05 2024 +0000 @@ -0,0 +1,231 @@ +import argparse +import utils.general_utils as utils +from typing import Optional, List +import os +import numpy as np +import pandas as pd +import cobra +import sys +import csv +from joblib import Parallel, delayed, cpu_count + +################################# process args ############################### +def process_args(args :List[str]) -> argparse.Namespace: + """ + Processes command-line arguments. + + Args: + args (list): List of command-line arguments. + + Returns: + Namespace: An object containing parsed arguments. + """ + parser = argparse.ArgumentParser(usage = '%(prog)s [options]', + description = 'process some value\'s') + + parser.add_argument( + '-ms', '--model_selector', + type = utils.Model, default = utils.Model.ENGRO2, choices = [utils.Model.ENGRO2, utils.Model.Custom], + help = 'chose which type of model you want use') + + parser.add_argument("-mo", "--model", type = str, + help = "path to input file with custom rules, if provided") + + parser.add_argument("-mn", "--model_name", type = str, help = "custom mode name") + + parser.add_argument( + '-mes', '--medium_selector', + default = "allOpen", + help = 'chose which type of medium you want use') + + parser.add_argument("-meo", "--medium", type = str, + help = "path to input file with custom medium, if provided") + + parser.add_argument('-ol', '--out_log', + help = "Output log") + + parser.add_argument('-td', '--tool_dir', + type = str, + required = True, + help = 'your tool directory') + + parser.add_argument('-ir', '--input_ras', + type=str, + required = False, + help = 'input ras') + + parser.add_argument('-rs', '--ras_selector', + required = True, + type=utils.Bool("using_RAS"), + help = 'ras selector') + + ARGS = parser.parse_args() + return ARGS + +########################### warning ########################################### +def warning(s :str) -> None: + """ + Log a warning message to an output log file and print it to the console. + + Args: + s (str): The warning message to be logged and printed. + + Returns: + None + """ + with open(ARGS.out_log, 'a') as log: + log.write(s + "\n\n") + print(s) + +############################ dataset input #################################### +def read_dataset(data :str, name :str) -> pd.DataFrame: + """ + Read a dataset from a CSV file and return it as a pandas DataFrame. + + Args: + data (str): Path to the CSV file containing the dataset. + name (str): Name of the dataset, used in error messages. + + Returns: + pandas.DataFrame: DataFrame containing the dataset. + + Raises: + pd.errors.EmptyDataError: If the CSV file is empty. + sys.exit: If the CSV file has the wrong format, the execution is aborted. + """ + 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 + + +def apply_ras_bounds(model, ras_row, rxns_ids): + """ + Adjust the bounds of reactions in the model based on RAS values. + + Args: + model (cobra.Model): The metabolic model to be modified. + ras_row (pd.Series): A row from a RAS DataFrame containing scaling factors for reaction bounds. + rxns_ids (list of str): List of reaction IDs to which the scaling factors will be applied. + + Returns: + None + """ + for reaction in rxns_ids: + if reaction in ras_row.index and pd.notna(ras_row[reaction]): + rxn = model.reactions.get_by_id(reaction) + scaling_factor = ras_row[reaction] + rxn.lower_bound *= scaling_factor + rxn.upper_bound *= scaling_factor + +def process_ras_cell(cellName, ras_row, model, rxns_ids, output_folder): + """ + Process a single RAS cell, apply bounds, and save the bounds to a CSV file. + + Args: + cellName (str): The name of the RAS cell (used for naming the output file). + ras_row (pd.Series): A row from a RAS DataFrame containing scaling factors for reaction bounds. + model (cobra.Model): The metabolic model to be modified. + rxns_ids (list of str): List of reaction IDs to which the scaling factors will be applied. + output_folder (str): Folder path where the output CSV file will be saved. + + Returns: + None + """ + model_new = model.copy() + apply_ras_bounds(model_new, ras_row, rxns_ids) + bounds = pd.DataFrame([(rxn.lower_bound, rxn.upper_bound) for rxn in model_new.reactions], index=rxns_ids, columns=["lower_bound", "upper_bound"]) + bounds.to_csv(output_folder + cellName + ".csv", sep='\t', index=True) + +def generate_bounds(model: cobra.Model, medium: dict, ras=None, output_folder='output/') -> pd.DataFrame: + """ + Generate reaction bounds for a metabolic model based on medium conditions and optional RAS adjustments. + + Args: + model (cobra.Model): The metabolic model for which bounds will be generated. + medium (dict): A dictionary where keys are reaction IDs and values are the medium conditions. + ras (pd.DataFrame, optional): A DataFrame with RAS scaling factors for different cell types. Defaults to None. + output_folder (str, optional): Folder path where output CSV files will be saved. Defaults to 'output/'. + + Returns: + pd.DataFrame: DataFrame containing the bounds of reactions in the model. + """ + rxns_ids = [rxn.id for rxn in model.reactions] + + # Set medium conditions + for reaction, value in medium.items(): + if value is not None: + model.reactions.get_by_id(reaction).lower_bound = -float(value) + + # Perform Flux Variability Analysis (FVA) + df_FVA = cobra.flux_analysis.flux_variability_analysis(model, fraction_of_optimum=0, processes=1).round(8) + + # Set FVA bounds + for reaction in rxns_ids: + rxn = model.reactions.get_by_id(reaction) + rxn.lower_bound = float(df_FVA.loc[reaction, "minimum"]) + rxn.upper_bound = float(df_FVA.loc[reaction, "maximum"]) + + if ras is not None: + Parallel(n_jobs=cpu_count())(delayed(process_ras_cell)(cellName, ras_row, model, rxns_ids, output_folder) for cellName, ras_row in ras.iterrows()) + else: + model_new = model.copy() + apply_ras_bounds(model_new, pd.Series([1]*len(rxns_ids), index=rxns_ids), rxns_ids) + bounds = pd.DataFrame([(rxn.lower_bound, rxn.upper_bound) for rxn in model_new.reactions], index=rxns_ids, columns=["lower_bound", "upper_bound"]) + bounds.to_csv(output_folder + "bounds.csv", sep='\t', index=True) + + +############################# main ########################################### +def main() -> None: + """ + Initializes everything and sets the program in motion based on the fronted input arguments. + + Returns: + None + """ + if not os.path.exists('ras_to_bounds'): + os.makedirs('ras_to_bounds') + + + global ARGS + ARGS = process_args(sys.argv) + + ARGS.output_folder = 'ras_to_bounds/' + + if(ARGS.ras_selector == True): + ras = read_dataset(ARGS.input_ras, "ras dataset") + ras.replace("None", None, inplace=True) + ras.set_index("Reactions", drop=True, inplace=True) + ras = ras.T + ras = ras.astype(float) + + model_type :utils.Model = ARGS.model_selector + if model_type is utils.Model.Custom: + model = model_type.getCOBRAmodel(customPath = utils.FilePath.fromStrPath(ARGS.model), customExtension = utils.FilePath.fromStrPath(ARGS.model_name).ext) + else: + model = model_type.getCOBRAmodel(toolDir=ARGS.tool_dir) + + if(ARGS.medium_selector == "Custom"): + medium = read_dataset(ARGS.medium, "medium dataset") + medium.set_index(medium.columns[0], inplace=True) + medium = medium.astype(float) + medium = medium[medium.columns[0]].to_dict() + else: + df_mediums = pd.read_csv(ARGS.tool_dir + "/local/medium/medium.csv", index_col = 0) + ARGS.medium_selector = ARGS.medium_selector.replace("_", " ") + medium = df_mediums[[ARGS.medium_selector]] + medium = medium[ARGS.medium_selector].to_dict() + + if(ARGS.ras_selector == True): + generate_bounds(model, medium, ras = ras, output_folder=ARGS.output_folder) + else: + generate_bounds(model, medium, output_folder=ARGS.output_folder) + + pass + +############################################################################## +if __name__ == "__main__": + main() \ No newline at end of file diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/ras_to_bounds.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/marea_2/ras_to_bounds.xml Wed Sep 18 09:21:05 2024 +0000 @@ -0,0 +1,110 @@ + + + + marea_macros.xml + + + + numpy + pandas + cobra + lxml + joblib + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff -r 985cd6f42be4 -r 95cda56f01a9 marea_2/rps_generator.xml --- a/marea_2/rps_generator.xml Wed Sep 18 09:13:24 2024 +0000 +++ b/marea_2/rps_generator.xml Wed Sep 18 09:21:05 2024 +0000 @@ -1,4 +1,4 @@ - + - Reaction Propensity Scores computation marea_macros.xml