Merge pull request #2379 from CCThinnes/develop

Functions for MicroMap paper

src/visualization/metabolicCartography/addFluxFromFileWidthAndColor.m

@@ -0,0 +1,149 @@
+function [map, flux2, fluxMap] = addFluxFromFileWidthAndColor(map, csvFilePath)
+    % Visualizes fluxes on a CellDesigner map. Rxn line width is
+    % proportional to flux magnitude. Positive fluxes are displayed in
+    % shades of red, and negative fluxes in shades of indigo. Higher flux
+    % magnitudes have higher hue saturation.
+    %
+    % INPUTS:
+    %   map:           A parsed model structure generated by
+    %                  'transformXML2Map' function.
+    %   csvFilePath:   Path to the CSV file containing reaction IDs in the
+    %                  first column and fluxes in the second column. The first row
+    %                  contains the respective headers.
+    %
+    % OUTPUTS:
+    %   map:           Updated map with reaction fluxes and colors.
+    %   flux2:         Fluxes and line widths through all reactions.
+    %   fluxMap:       List of reactions carrying flux in the map and their
+    %                  corresponding line widths.
+    %
+    % .. Author: - Cyrille C. Thinnes. University of Galway, Ireland, 25/09/2024.
+
+    % Load the CSV data, ignoring headers.
+    fluxData = readtable(csvFilePath, 'ReadVariableNames', false);  
+
+    % Extract reaction IDs and flux values from the first and second columns.
+    reactionIDs = fluxData{:, 1};  % First column contains reaction IDs.
+    fluxValues = fluxData{:, 2};   % Second column contains flux values.
+
+    % Convert reaction IDs to cell array of strings if necessary.
+    if isnumeric(reactionIDs)
+        reactionIDs = cellstr(num2str(reactionIDs));
+    elseif isstring(reactionIDs)
+        reactionIDs = cellstr(reactionIDs);
+    end
+
+    % Obtain the non-zero flux reactions.
+    idx = fluxValues ~= 0;
+    rxn = reactionIDs(idx);  % Get reactions with non-zero flux.
+
+    % Obtain list of reactions carrying flux and their flux values.
+    flux2 = [rxn, num2cell(fluxValues(idx))];
+
+    % Normalize fluxes to give width.
+    absFlux = abs(fluxValues);
+    maxAbsFlux = max(absFlux);
+
+    % Avoid division by zero if maxAbsFlux is zero.
+    if maxAbsFlux == 0
+        rxnWidth = zeros(size(absFlux));
+    else
+        rxnWidth = absFlux / maxAbsFlux;  % Normalized flux values between 0 and 1
+    end
+
+    % Cap rxnWidth at 1
+    rxnWidth(rxnWidth > 1) = 1;
+
+    % Assign line widths based on non-linear gradations that reflect the decimal thresholds.
+    rxnWidthValues = ones(size(rxnWidth)); % Initialize with ones for zero fluxes.
+
+    % Define non-linear ranges for line widths with values reflecting the
+    % decimal thresholds.
+    rxnWidthValues(rxnWidth > 0 & rxnWidth <= 0.2) = 2;
+    rxnWidthValues(rxnWidth > 0.2 & rxnWidth <= 0.5) = 5;
+    rxnWidthValues(rxnWidth > 0.5 & rxnWidth <= 0.8) = 8;
+    rxnWidthValues(rxnWidth > 0.8 & rxnWidth <= 1) = 10;
+
+    flux2 = [flux2 num2cell(rxnWidthValues(idx))];
+
+    % Add specific width to each reaction in the map based on the fluxes.
+    map.rxnWidth = cell(size(map.rxnName));
+    for i = 1:length(map.rxnName)
+        a = find(strcmp(map.rxnName{i}, reactionIDs));
+        if isempty(a)
+            map.rxnWidth{i} = 1;
+        else
+            map.rxnWidth{i} = rxnWidthValues(a);
+        end
+    end
+
+    fluxMap = [map.rxnName, map.rxnWidth];
+
+    % Define color shades (4 shades each for red and
+    % indigo) as specified. Add FF before the hex code to account for the
+    % alpha channel. We used Open Color to inform the color choice
+    % (https://yeun.github.io/open-color/).
+    redShades = {'FFffc9c9', 'FFff8787', 'FFfa5252', 'FFc92a2a'}; % Lightest to darkest
+    indigoShades = {'FFbac8ff', 'FF748ffc', 'FF4c6ef5', 'FF364fc7'}; % Lightest to darkest
+    LIGHTGRAY = 'FFD3D3D3';
+    WHITE = 'FFFFFFFF';
+
+    % Map line widths to shade indices.
+    lineWidths = [1, 2, 5, 8, 10]; % Unique line widths.
+    shadeIndices = [1, 1, 2, 3, 4]; % Map line widths to shades.
+
+    % Create a mapping from line widths to shade indices.
+    lineWidthToShadeIndex = containers.Map(lineWidths, shadeIndices);
+
+    % Initialize reaction colors to LIGHTGRAY for reactions with zero flux.
+    idZeroFlux = fluxValues == 0;
+    rxZeroFlux = reactionIDs(idZeroFlux);
+    indexZeroFlux = ismember(map.rxnName, rxZeroFlux);
+    map.rxnColor(indexZeroFlux, 1) = {LIGHTGRAY}; % LIGHTGRAY
+
+    % Change all nodes' color to WHITE.
+    map.molColor = repmat({WHITE}, size(map.molAlias));
+
+    % For reactions with non-zero flux, assign colors based on line width
+    % and sign.
+    idNonZeroFlux = ~idZeroFlux;
+    rxNonZeroFlux = reactionIDs(idNonZeroFlux);
+    fluxNonZero = fluxValues(idNonZeroFlux);
+    rxnWidthsNonZero = rxnWidthValues(idNonZeroFlux);
+
+    % Now, for each reaction with non-zero flux.
+    for i = 1:length(rxNonZeroFlux)
+        rxnName = rxNonZeroFlux{i};
+        idxMap = find(strcmp(map.rxnName, rxnName));
+        if isempty(idxMap)
+            continue;
+        end
+        lineWidth = rxnWidthsNonZero(i);
+        shadeIdx = lineWidthToShadeIndex(lineWidth); % Get the shade index corresponding to the line width.
+
+        if fluxNonZero(i) > 0
+            % Positive flux, use red shades.
+            color = redShades{shadeIdx};
+        else
+            % Negative flux, use indigo shades.
+            color = indigoShades{shadeIdx};
+        end
+        map.rxnColor{idxMap, 1} = color;
+
+        % Also change color of reactants and products accordingly.
+        % Reactants.
+        reactants = [map.rxnBaseReactantAlias{idxMap}; map.rxnReactantAlias{idxMap}];
+        for j = 1:length(reactants)
+            a = reactants{j};
+            idMol = strcmp(map.molAlias, a);
+            map.molColor{idMol, 1} = color;
+        end
+        % Products.
+        products = [map.rxnBaseProductAlias{idxMap}; map.rxnProductAlias{idxMap}];
+        for j = 1:length(products)
+            a = products{j};
+            idMol = strcmp(map.molAlias, a);
+            map.molColor{idMol, 1} = color;
+        end
+    end
+end

src/visualization/metabolicCartography/addFluxWidthAndColor.m

@@ -0,0 +1,135 @@
+function [map, flux2, fluxMap] = addFluxWidthAndColor(map, reactionIDs, fluxValues)
+    % Function to add flux widths and corresponding color shades to a map based on flux values.
+    % Rxn line width is
+    % proportional to flux magnitude. Positive fluxes are displayed in
+    % shades of red, and negative fluxes in shades of indigo. Higher flux
+    % magnitudes have higher hue saturation.
+    %
+    % INPUTS:
+    %   map:           A parsed model structure generated by
+    %                  'transformXML2Map' function.
+    %   reactionIDs:   Cell array of reaction IDs.
+    %   fluxValues:    Array of flux values corresponding to the reaction IDs.
+    %
+    % OUTPUTS:
+    %   map:           Updated map with reaction fluxes and colors.
+    %   flux2:         Fluxes and line widths through all reactions.
+    %   fluxMap:       List of reactions carrying flux in the map and their
+    %                  corresponding line widths.
+    %
+    % .. Author: - Cyrille C. Thinnes. University of Galway, Ireland, 26/09/2024.
+
+    % Ensure reaction IDs are cell array of strings
+    if isnumeric(reactionIDs)
+        reactionIDs = cellstr(num2str(reactionIDs));
+    elseif isstring(reactionIDs)
+        reactionIDs = cellstr(reactionIDs);
+    end
+
+    % Obtain the non-zero flux reactions.
+    idx = fluxValues ~= 0;
+    rxn = reactionIDs(idx);  % Reactions with non-zero flux.
+
+    % Obtain list of reactions carrying flux and their flux values.
+    flux2 = [rxn, num2cell(fluxValues(idx))];
+
+    % Normalize fluxes to assign widths.
+    absFlux = abs(fluxValues);
+    maxAbsFlux = max(absFlux);
+
+    % Avoid division by zero if maxAbsFlux is zero
+    if maxAbsFlux == 0
+        rxnWidth = zeros(size(absFlux));
+    else
+        rxnWidth = absFlux / maxAbsFlux;  % Normalized between 0 and 1.
+    end
+
+    % Cap rxnWidth at 1
+    rxnWidth(rxnWidth > 1) = 1;
+
+    % Assign line widths based on thresholds.
+    rxnWidthValues = ones(size(rxnWidth)); % Initialize with ones for zero fluxes.
+
+    % Define thresholds for line widths.
+    rxnWidthValues(rxnWidth > 0 & rxnWidth <= 0.2) = 2;
+    rxnWidthValues(rxnWidth > 0.2 & rxnWidth <= 0.5) = 5;
+    rxnWidthValues(rxnWidth > 0.5 & rxnWidth <= 0.8) = 8;
+    rxnWidthValues(rxnWidth > 0.8 & rxnWidth <= 1) = 10;
+
+    flux2 = [flux2, num2cell(rxnWidthValues(idx))];
+
+    % Add specific widths to reactions in the map.
+    map.rxnWidth = cell(size(map.rxnName));
+    for i = 1:length(map.rxnName)
+        a = find(strcmp(map.rxnName{i}, reactionIDs));
+        if isempty(a)
+            map.rxnWidth{i} = 1;  % Default width for reactions not in the flux data.
+        else
+            map.rxnWidth{i} = rxnWidthValues(a);
+        end
+    end
+
+    fluxMap = [map.rxnName, map.rxnWidth];
+
+    % Define color shades with alpha channel (FF for full opacity).
+    redShades = {'FFffc9c9', 'FFff8787', 'FFfa5252', 'FFc92a2a'}; % Lightest to darkest.
+    indigoShades = {'FFbac8ff', 'FF748ffc', 'FF4c6ef5', 'FF364fc7'}; % Lightest to darkest.
+    LIGHTGRAY = 'FFD3D3D3';
+    WHITE = 'FFFFFFFF';
+
+    % Map line widths to shade indices.
+    lineWidths = [1, 2, 5, 8, 10];  % Unique line widths.
+    shadeIndices = [1, 1, 2, 3, 4]; % Corresponding shade indices.
+
+    % Create mapping from line widths to shade indices.
+    lineWidthToShadeIndex = containers.Map(lineWidths, shadeIndices);
+
+    % Initialize reaction colors to LIGHTGRAY for reactions with zero flux.
+    idZeroFlux = fluxValues == 0;
+    rxZeroFlux = reactionIDs(idZeroFlux);
+    indexZeroFlux = ismember(map.rxnName, rxZeroFlux);
+    map.rxnColor(indexZeroFlux, 1) = {LIGHTGRAY};
+
+    % Set all metabolite colors to WHITE.
+    map.molColor = repmat({WHITE}, size(map.molAlias));
+
+    % Assign colors to reactions with non-zero flux.
+    idNonZeroFlux = ~idZeroFlux;
+    rxNonZeroFlux = reactionIDs(idNonZeroFlux);
+    fluxNonZero = fluxValues(idNonZeroFlux);
+    rxnWidthsNonZero = rxnWidthValues(idNonZeroFlux);
+
+    % Loop over reactions with non-zero flux.
+    for i = 1:length(rxNonZeroFlux)
+        rxnName = rxNonZeroFlux{i};
+        idxMap = find(strcmp(map.rxnName, rxnName));
+        if isempty(idxMap)
+            continue;  % Skip if reaction not in map.
+        end
+        lineWidth = rxnWidthsNonZero(i);
+        shadeIdx = lineWidthToShadeIndex(lineWidth);
+
+        if fluxNonZero(i) > 0
+            % Positive flux: use red shades.
+            color = redShades{shadeIdx};
+        else
+            % Negative flux: use indigo shades.
+            color = indigoShades{shadeIdx};
+        end
+        map.rxnColor{idxMap, 1} = color;
+
+        % Update colors of reactants and products.
+        reactants = [map.rxnBaseReactantAlias{idxMap}; map.rxnReactantAlias{idxMap}];
+        for j = 1:length(reactants)
+            a = reactants{j};
+            idMol = strcmp(map.molAlias, a);
+            map.molColor{idMol, 1} = color;
+        end
+        products = [map.rxnBaseProductAlias{idxMap}; map.rxnProductAlias{idxMap}];
+        for j = 1:length(products)
+            a = products{j};
+            idMol = strcmp(map.molAlias, a);
+            map.molColor{idMol, 1} = color;
+        end
+    end
+end

src/visualization/metabolicCartography/uniqueMetabolites.m

@@ -0,0 +1,25 @@
+function uniqueMets = uniqueMetabolites(model)
+% uniqueMetabolites - Identifies unique metabolites by ignoring compartment tags
+%
+% USAGE:
+%
+%    uniqueMets = uniqueMetabolites(model)
+%
+% INPUT:
+%    model:         COBRA model structure containing metabolites in model.mets
+%
+% OUTPUT:
+%    uniqueMets:    Cell array of unique metabolite names, excluding compartment tags
+%
+% .. Authors:
+%    Cyrille Thinnes, University of Galway, 25/10/2024
+
+% Extract metabolites
+mets = model.mets;
+
+% Remove compartment tags, e.g., 'metabolite[c]' -> 'metabolite'
+metsNoCompartment = regexprep(mets, '\[.*?\]', '');
+
+% Identify unique metabolites
+uniqueMets = unique(metsNoCompartment);
+end

src/visualization/metabolicCartography/uniqueSpeciesInMap.m

@@ -0,0 +1,43 @@
+function uniqueSpecies = uniqueSpeciesInMap(mapMicroMap)
+% uniqueSpeciesInMap - Identifies unique metabolites and other species in a CellDesigner map structure
+%
+% USAGE:
+%
+%    uniqueSpecies = uniqueSpeciesInMap(mapMicroMap)
+%
+% INPUT:
+%    mapMicroMap:    Structure containing species from a CellDesigner map
+%                    mapMicroMap.specName contains species names
+%
+% OUTPUT:
+%    uniqueSpecies:  Structure with fields:
+%                    - mets: Unique metabolites (names without compartment tags)
+%                    - nonMets: Unique non-metabolite species
+%
+% .. Authors:
+%       Cyrille Thinnes, University of Galway, 25/10/2024
+
+% Extract species names
+specNames = mapMicroMap.specName;
+
+% Identify metabolites with compartment tags
+metabolites = specNames(contains(specNames, '[') & contains(specNames, ']'));
+
+% Remove compartment tags for metabolites
+metsNoCompartment = regexprep(metabolites, '\[.*?\]', '');
+
+% Find unique metabolites
+uniqueMets = unique(metsNoCompartment);
+
+% Identify non-metabolite species (without compartment tags)
+nonMetabolites = specNames(~contains(specNames, '[') & ~contains(specNames, ']'));
+
+% Find unique non-metabolite species
+uniqueNonMets = unique(nonMetabolites);
+
+% Compile results into the output structure uniqueSpecies
+uniqueSpecies.mets = uniqueMets;
+uniqueSpecies.nonMets = uniqueNonMets;
+
+end
+

src/visualization/metabolicCartography/visualizeFluxFromFile.m

@@ -0,0 +1,32 @@
+function visualizeFluxFromFile(mapXMLFile, fluxCSVFile, outputXMLFile)
+    % This function takes a CellDesigner map XML file and a flux vector table file (e.g., CSV or XLSX),
+    % visualizes the flux on the map, and saves the output map as a CellDesigner XML file.
+    %
+    % INPUTS:
+    %   mapXMLFile:   The input CellDesigner .xml file for the map.
+    %   fluxCSVFile:  The input file containing the flux vector, e.g., CSV
+    %   or XLSX.
+    %   outputXMLFile: (Optional) The output .xml file name for the updated map.
+    %
+    % OUTPUTS:
+    %   The map with visualized flux saved as the specified CellDesigner XML file.
+
+    % .. Author: - Cyrille C. Thinnes. University of Galway, Ireland, 24/09/2024.
+
+    % Set default output XML file name if not provided.
+    if nargin < 3 || isempty(outputXMLFile)
+        outputXMLFile = 'VisualizedFluxOnMap.xml';
+    end
+
+    % Parse the map CellDesigner .xml file into Matlab.
+    [xmlInputMap, mapInputMap] = transformXML2Map(mapXMLFile);
+
+    % Visualise the flux vector to highlight flux magnitude and sign.
+    mapWidthAndColourFlux = addFluxFromFileWidthAndColor(mapInputMap, fluxCSVFile);
+
+    % Reconstitute the CellDesigner map with the flux vector.
+    transformMap2XML(xmlInputMap, mapWidthAndColourFlux, outputXMLFile);
+
+    % Display success message
+    disp(['Map with flux visualisation saved as: ', outputXMLFile]);
+end