Skip to content

Main Functions

TackleMMe functions are sorted into three categories:

  1. Model Building and Project Initialization

  2. Single Model Analysis

  3. Model Comparison

Main functions are documented below.

Model Building and Project Initialization

createProject(params)

Creates a project structure ready for the analysis pipeline.

This function initializes a project object containing one or more models, each defined by a set of parameters. The project can then be used as input for single model analysis and model comparison functions.

Input arguments:

Name Type Description Default
params cell

1-by-N cell array. One struct per model. See the Project Initialization page for the full list of available fields. Required fields: modelName, contextSpecificModel.

 required

Output arguments:

Name Type Description
project struct

Initialized project with a models field.

Examples:

params = {struct('modelName', 'model1', 'contextSpecificModel', model)};
project = createProject(params);

params = {struct(...), struct(...)};
project = createProject(params);
Note

Fields marked as (rFastcormics) in the source code are specific to models built with rFASTCORMICS and are optional for other COBRA models. Only modelName and contextSpecificModel are required. Field validation is performed by validateParamsForPipeline.

Source code in scr/classes/project_initialization/createProject.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
function project = createProject(params)
% Creates a project structure ready for the analysis pipeline.
%
% This function initializes a project object containing one or more models,
% each defined by a set of parameters. The project can then be used as input
% for single model analysis and model comparison functions.
%
% Arguments:
%   params (cell): 1-by-N cell array. One struct per model. See the Project Initialization
%       page for the full list of available fields. Required fields:
%       modelName, contextSpecificModel.
%
% Returns:
%   project (struct): Initialized project with a models field.
%
% Examples:
%   ```matlab
%   params = {struct('modelName', 'model1', 'contextSpecificModel', model)};
%   project = createProject(params);
%
%   params = {struct(...), struct(...)};
%   project = createProject(params);
%   ```
%
% Note:
%   Fields marked as (rFastcormics) in the source code are specific to
%   models built with rFASTCORMICS and are optional for other COBRA models.
%   Only `modelName` and `contextSpecificModel` are required. Field
%   validation is performed by `validateParamsForPipeline`.

arguments
    params (1,:) cell
end

% Initialize project
project = struct();
project.models = struct();

% Loop through models
for i = 1:numel(params)
    paramsForModel = params{i};

    % Validate each struct's fields
    paramsForModel = validateParamsForPipeline(paramsForModel);

    % Initiate a struct per model
    project.models.(paramsForModel.modelName) = formatParamsForModel(paramsForModel);
end

end

addModelsToProject(project, params)

Adds one or more models to an already existing project.

This function extends an existing project by adding new models. Each model is validated and formatted before being added. If a model with the same name already exists, the user is prompted to confirm overwriting.

Input arguments:

Name Type Description Default
project struct

Existing project structure created by createProject.

 required
params cell

1-by-N cell array, one struct per model to add. See the Project Initialization page for the full list of available fields. Required fields: modelName, contextSpecificModel.

 required

Output arguments:

Name Type Description
project struct

The input project with new models added.

Examples:

project = addModelsToProject(project, ...
    {struct('modelName', 'model2', 'contextSpecificModel', model2)});

% Add multiple models
project = addModelsToProject(project, {struct(...), struct(...)});
Note

The project format is validated with checkProjectFormat before adding any model. Fields specific to rFASTCORMICS are optional.

Source code in scr/classes/project_initialization/addModelsToProject.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
function project = addModelsToProject(project, params)
% Adds one or more models to an already existing project.
%
% This function extends an existing project by adding new models. Each
% model is validated and formatted before being added. If a model with
% the same name already exists, the user is prompted to confirm
% overwriting.
%
% Arguments:
%   project (struct): Existing project structure created by createProject.
%   params (cell): 1-by-N cell array, one struct per model to add. See
%       the Project Initialization page for the full list of available
%       fields. Required fields: modelName, contextSpecificModel.
%
% Returns:
%   project (struct): The input project with new models added.
%
% Examples:
%   ```matlab
%   project = addModelsToProject(project, ...
%       {struct('modelName', 'model2', 'contextSpecificModel', model2)});
%
%   % Add multiple models
%   project = addModelsToProject(project, {struct(...), struct(...)});
%   ```
%
% Note:
%   The project format is validated with checkProjectFormat before
%   adding any model. Fields specific to rFASTCORMICS are optional.

arguments
    project
    params (1,:) cell
end

% Check whether the project is in the correct format
fprintf("Checking project format before adding model. \n")
checkProjectFormat(project);

% Add extra models one by one
for i = 1:numel(params)
    paramsForModel = params{i};

    % Validate each struct's fields
    paramsForModel = validateParamsForPipeline(paramsForModel);

    % Initiate a struct per model
    if isfield(project.models, paramsForModel.modelName)
        answer = input("Model '" + paramsForModel.modelName + "' already exists. Overwrite? [y/n]: ", 's');
        if ~strcmpi(answer, 'y')
            fprintf("Model '" + paramsForModel.modelName + "' skipped.")
            continue
        end
    end

    project.models.(paramsForModel.modelName) = formatParamsForModel(paramsForModel);

end

end

Single Model Analysis

singleModelAnalysis(project, parameterTable, modelList={}, analyses={}, saveCheckpoint=true, resumeFromCheckpoint=false)

Runs analyses on one or multiple models and stores results.

This function performs one or several analyses on each model in the provided list and stores the results under a timestamped analysis ID in the project structure. A checkpoint system is available to resume after a crash.

Input arguments:

Name Type Description Default
project struct

Project structure created by createProject.

 required
parameterTable table

Parameter table defining analysis settings. See the Single Model Analysis page for details.

 required
modelList cell

Names of the models to analyze. If empty, all models in the project are used.

 required
analyses cell

List of analyses to perform. If empty, all available analyses are run. Valid keys: FBA, FVA, sampling, loopless, kld, singleGeneDeletion, doubleGeneDeletion.

 required
saveCheckpoint logical

Whether to save a checkpoint after each model. Default: true.

 required
resumeFromCheckpoint logical

Whether to resume from the last saved checkpoint. Default: false.

 required

Output arguments:

Name Type Description
project struct

The input project with an analysis field added to each analyzed model.

Examples:

% Run all analyses on all models
project = singleModelAnalysis(project, parameterTable);

% Run FBA and sampling on specific models
project = singleModelAnalysis(project, parameterTable, ...
    {"model1", "model2"}, {"FBA", "sampling"});

% Resume after a crash
project = singleModelAnalysis(project, parameterTable, ...
    resumeFromCheckpoint = true);
Note

When loopless is requested without sampling, a samplingToUse parameter must be provided in the parameter table, referencing a previous sampling analysis ID. The IDs must be listed in the same order as the models in modelList.

Warning

Unimplemented analyses requested in the analyses list are skipped with a console warning.

Source code in scr/classes/single_model_analysis/singleModelAnalysis.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
function project = singleModelAnalysis(project, parameterTable, modelList, analyses, saveCheckpoint, resumeFromCheckpoint)
% Runs analyses on one or multiple models and stores results.
%
% This function performs one or several analyses on each model in the
% provided list and stores the results under a timestamped analysis ID
% in the project structure. A checkpoint system is available to resume
% after a crash.
%
% Arguments:
%   project (struct): Project structure created by createProject.
%   parameterTable (table): Parameter table defining analysis settings.
%       See the Single Model Analysis page for details.
%   modelList (cell): Names of the models to analyze. If empty, all
%       models in the project are used.
%   analyses (cell): List of analyses to perform. If empty, all
%       available analyses are run. Valid keys: FBA, FVA, sampling,
%       loopless, kld, singleGeneDeletion, doubleGeneDeletion.
%   saveCheckpoint (logical): Whether to save a checkpoint after each
%       model. Default: true.
%   resumeFromCheckpoint (logical): Whether to resume from the last
%       saved checkpoint. Default: false.
%
% Returns:
%   project (struct): The input project with an analysis field added
%       to each analyzed model.
%
% Examples:
%   ```matlab
%   % Run all analyses on all models
%   project = singleModelAnalysis(project, parameterTable);
%
%   % Run FBA and sampling on specific models
%   project = singleModelAnalysis(project, parameterTable, ...
%       {"model1", "model2"}, {"FBA", "sampling"});
%
%   % Resume after a crash
%   project = singleModelAnalysis(project, parameterTable, ...
%       resumeFromCheckpoint = true);
%   ```
%
% Note:
%   When loopless is requested without sampling, a samplingToUse
%   parameter must be provided in the parameter table, referencing a
%   previous sampling analysis ID. The IDs must be listed in the same
%   order as the models in modelList.
%
% Warning:
%   Unimplemented analyses requested in the analyses list are skipped
%   with a console warning.
%
arguments
    project struct
    parameterTable table
    modelList {mustBeText} = {}
    analyses cell {mustBeVectorOrEmpty} = {}
    saveCheckpoint logical = true
    resumeFromCheckpoint logical = false
end

%% Check the structure of project
modelList = checkStructureForSingleModelAnalysis(project, modelList); % get valid models

%% Check format of parameter table
checkParamTableFormat(parameterTable)

%% Check that given analysis are part of the list
validAnalyses = {'FBA', 'FVA', 'sampling', 'loopless', 'kld', 'singleGeneDeletion', 'doubleGeneDeletion'};

if isempty(analyses) || (ischar(analyses) && isempty(analyses))
    toPerform = validAnalyses;
else
    isImplemented = ismember(analyses, validAnalyses);
    if ~all(isImplemented)
        notImplemented = analyses(~isImplemented);
        fprintf('The following analyses are not implemented and will not be performed:');
        for i = 1:length(notImplemented)
            fprintf(' - %s', notImplemented{i});
        end
    end
    toPerform = analyses(isImplemented);
end

% check that the loopless ID given are valid analysisIDs 
% check that the given analysis IDs for the loopless are actually
% available in the corresponding models 
if any("loopless" == analyses) && ~any("sampling" == analyses)
    % check that we have as many samplingIDs given as models
    analysisIDSampling = parameterTable.Value{parameterTable.Parameter == "samplingToUse" & parameterTable.Analysis == "loopless"};
    analysisID = split(string(analysisIDSampling),"/");
    assert(length(analysisID) == length(modelList));
    % check that the sampling analysisIDs are in the same order as the
    % models
    for m = 1:numel(modelList)
        mod = string(modelList(m));
        anaID = analysisID(m);
        if isfield(project.models.(mod),"analysis")
            if ~ismember(anaID, fieldnames(project.models.(mod).analysis))
                error("The analysisID (%s) for model %s can not be found! Check again the analysis IDs (samplingToUse parameter) you specified in the defaultParametersAnalysis.csv. The analysisIDs should be in the same order as their corresponding models are given in the modelList.",...
                      anaID,mod)
            end
        end
    end
end

% Checkpoint: optionally resume from last saved model
checkpointFile = 'singleModelAnalysis_checkpoint.mat';
startIdx = 1;

if resumeFromCheckpoint && exist(checkpointFile, 'file')
    loaded = load(checkpointFile);
    project = loaded.project;
    startIdx = loaded.i + 1;
    if startIdx <= numel(modelList)
        fprintf('Checkpoint loaded — resuming at model %d/%d (%s)\n', ...
                startIdx, numel(modelList), modelList{startIdx});
    end
end

%% Run analysis
for i = startIdx:numel(modelList)
    name = modelList{i};
    fprintf('Analysis running for model %s (%d/%d)\n', name, i, numel(modelList));

    if ~isfield(project.models.(name), 'analysis')
        project.models.(name).analysis = struct();
    end

    % Analysis id
    id = ['analysis_' char(datetime("now", "Format", "yyyyMMdd_HHmm"))];
    project.models.(name).analysis.(id) = struct();

    % select analysisID used for sampling in case loopless is run 
    if any("loopless" == analyses) & ~any("sampling" == analyses)
        parameterTable.Value{parameterTable.Parameter == "samplingToUse" & parameterTable.Analysis == "loopless"} = analysisID(i);
    else % if the sampling is run in the same analysis run, that sampling will be automatically used for the loopless sampling
        parameterTable.Value{parameterTable.Parameter == "samplingToUse" & parameterTable.Analysis == "loopless"} = id;
    end
    % Store the settings
    project.models.(name).analysis.(id).parameters = parameterTable;

    % Checkpoint: try/catch to save state before crash propagates
    try
        project = performAnalysis(project, parameterTable, name, toPerform, id);
    catch ME
        fprintf('Error on model %s : %s', name, ME.message);
        fprintf('Stack:');
        for s = 1:numel(ME.stack)
            fprintf('  %s (line %d)', ME.stack(s).name, ME.stack(s).line);
        end
        fprintf('Last valid checkpoint: model %d (%s).\n', i-1, modelList{i-1});
        fprintf('Partial project returned. Relaunch with resumeFromCheckpoint=true to resume.\n');
        return;
    end

    % Checkpoint: save after each successful model
    if saveCheckpoint
        save(checkpointFile, 'project', 'i', '-v7.3');
        fprintf('Checkpoint saved (model %d/%d)\n', i, numel(modelList));
    end
    % End checkpoint save

end

%% Checkpoint: clean up if all models completed
if saveCheckpoint && exist(checkpointFile, 'file')
    delete(checkpointFile);
    fprintf('All models analyzed — checkpoint deleted.');
end

end

addAnalysisToExistingOne(project, parameterTable, modelName, analyses, analysisId)

Adds analyses to an already existing analysis run.

This function adds one or more analyses to an existing analysis field identified by its analysis ID, without creating a new timestamped entry. It is useful for supplementing a previous run with additional analyses or re-running existing ones with updated parameters.

Input arguments:

Name Type Description Default
project struct

Project structure with an existing analysis field.

 required
parameterTable table

Parameter table containing settings for the analyses to add.

 required
modelName char

Name of the model.

 required
analyses char

Analysis key(s) to add (e.g. 'FBA', 'sampling'). Can be a single string or a cell array of strings.

 required
analysisId char

Existing analysis ID to add to (e.g. 'analysis_20240815_1430').

 required

Output arguments:

Name Type Description
project struct

The input project with updated analysis results.

Examples:

% Add gene deletion analyses to an existing run
project = addAnalysisToExistingOne(project, parameterTable, ...
    'model1', {'singleGeneDeletion', 'doubleGeneDeletion'}, ...
    'analysis_20240815_1430');

% Re-run FVA with updated parameters
project = addAnalysisToExistingOne(project, parameterTable, ...
    'model1', 'FVA', 'analysis_20240815_1430');
Warning

When re-running an existing analysis, results and the corresponding rows in the stored parameters table are overwritten after user confirmation. Parameters for other analyses are preserved.

Source code in scr/classes/single_model_analysis/addAnalysisToExistingOne.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
function project = addAnalysisToExistingOne(project, parameterTable, modelName, analyses, analysisId)
% Adds analyses to an already existing analysis run.
%
% This function adds one or more analyses to an existing analysis field
% identified by its analysis ID, without creating a new timestamped
% entry. It is useful for supplementing a previous run with additional
% analyses or re-running existing ones with updated parameters.
%
% Arguments:
%   project (struct): Project structure with an existing analysis field.
%   parameterTable (table): Parameter table containing settings for the
%       analyses to add.
%   modelName (char): Name of the model.
%   analyses (char): Analysis key(s) to add (e.g. 'FBA', 'sampling').
%       Can be a single string or a cell array of strings.
%   analysisId (char): Existing analysis ID to add to (e.g.
%       'analysis_20240815_1430').
%
% Returns:
%   project (struct): The input project with updated analysis results.
%
% Examples:
%   ```matlab
%   % Add gene deletion analyses to an existing run
%   project = addAnalysisToExistingOne(project, parameterTable, ...
%       'model1', {'singleGeneDeletion', 'doubleGeneDeletion'}, ...
%       'analysis_20240815_1430');
%
%   % Re-run FVA with updated parameters
%   project = addAnalysisToExistingOne(project, parameterTable, ...
%       'model1', 'FVA', 'analysis_20240815_1430');
%   ```
%
% Warning:
%   When re-running an existing analysis, results and the corresponding
%   rows in the stored parameters table are overwritten after user
%   confirmation. Parameters for other analyses are preserved.
%
arguments
    project struct
    parameterTable table
    modelName {mustBeText}
    analyses {mustBeText}
    analysisId {mustBeText}
end

if ischar(analyses)
    analyses = {analyses};
elseif isstring(analyses)
    analyses = cellstr(analyses);
end

validModels = checkStructureForSingleModelAnalysis(project, {modelName});
if isequal(char(validModels(1)), modelName)
    if isfield(project.models.(modelName), 'analysis')
        if isfield(project.models.(modelName).analysis, analysisId)
            if isstruct(project.models.(modelName).analysis.(analysisId))
                for i = 1:numel(analyses)
                    test = char(analyses(i));
                    % Filter parameterTable to only keep rows for the current test
                    testRows = strcmp(parameterTable.Analysis, test);
                    testParams = parameterTable(testRows, :);
                    if isempty(testParams)
                        error('No parameters corresponding to %s were found in the parameterTable', test);
                    else
                        if isfield(project.models.(modelName).analysis.(analysisId), test)
                            % check par rapport à table
                            if isfield(project.models.(modelName).analysis.(analysisId), 'parameters')                            
                                % Also filter existing parameters for the same test
                                existingParams = project.models.(modelName).analysis.(analysisId).parameters;
                                existingTestRows = strcmp(existingParams.Analysis, test);
                                existingTestParams = existingParams(existingTestRows, :);

                                % If no existing params for this test, just add them
                                if isempty(existingTestParams)
                                    answer = input(sprintf('Analysis "%s" already exists without given parameters. Overwrite existing analysis and store new parameters? (y/n): ', analysisId), 's');
                                    if ~strcmpi(answer, 'y')
                                        error('Analysis adding cancelled by user.');
                                    else
                                        project.models.(modelName).analysis.(analysisId).parameters = [existingParams; testParams];
                                        project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);
                                    end

                                elseif isequal(existingTestParams, testParams)
                                    dispValue = evalc('disp(testParams)');
                                    answer = input(sprintf('Analysis "%s" already exists with identical parameters:\n%s\nOverwrite? (y/n): ', test, dispValue), 's');
                                    if ~strcmpi(answer, 'y')
                                        error('Analysis adding cancelled by user.');
                                    else
                                        project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);
                                    end

                                else
                                    % Parameters differ — show differences
                                    commonParams = intersect(existingTestParams.Parameter, testParams.Parameter);
                                    differingParams = {};

                                    for i = 1:numel(commonParams)
                                        param = commonParams{i};
                                        oldVal = existingTestParams.Value{strcmp(existingTestParams.Parameter, param)};
                                        newVal = testParams.Value{strcmp(testParams.Parameter, param)};
                                        if ~strcmp(oldVal, newVal)
                                            differingParams{end+1} = param;
                                        end
                                    end

                                    fprintf('The following parameters differ between the existing and new analysis "%s":\n', analysisId);
                                    fprintf('  %-25s  %-25s  %-25s ', 'Parameter', 'Existing Value', 'New Value');
                                    fprintf('\n  %s\n', repmat('-', 1, 77));
                                    fprintf('');
                                    for i = 1:numel(differingParams)
                                        param = differingParams{i};
                                        oldVal = existingTestParams.Value{strcmp(existingTestParams.Parameter, param)};
                                        newVal = testParams.Value{strcmp(testParams.Parameter, param)};
                                        fprintf('  %-25s  %-25s  %-25s ', param, oldVal, newVal);
                                    end
                                    fprintf('');

                                    answer = input('\nOverwrite existing analysis with new parameters? (y/n): ', 's');
                                    if ~strcmpi(answer, 'y')
                                        error('Analysis adding cancelled by user.');
                                    else
                                        % Overwrite only the test-related rows in the parameters table
                                        otherRows = ~strcmp(existingParams.Analysis, test);
                                        otherParams = existingParams(otherRows, :);
                                        project.models.(modelName).analysis.(analysisId).parameters = [otherParams; testParams];

                                        % Re-run and add analysis
                                        project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);
                                    end
                                end
                            else
                                warning('No parameters table found. Continuing and storing given parameters.')
                                project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);
                                project.models.(modelName).analysis.(analysisId).parameters = testParams;
                            end
                        else
                            fprintf('Analysis "%s" does not exist yet. It will be run and stored with its corresponding parameters.\n', test);

                            if isfield(project.models.(modelName).analysis.(analysisId), 'parameters')                            
                                % Also filter existing parameters for the same test
                                existingParams = project.models.(modelName).analysis.(analysisId).parameters;
                                existingTestRows = strcmp(existingParams.Analysis, test);
                                existingTestParams = existingParams(existingTestRows, :);

                                % If no existing params for this test, just add them
                                if isempty(existingTestParams)
                                    project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);
                                    project.models.(modelName).analysis.(analysisId).parameters = [existingParams; testParams];

                                elseif isequal(existingTestParams, testParams)
                                    project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);

                                else
                                    % Parameters differ
                                    % Overwrite only the test-related rows in the parameters table
                                    otherRows = ~strcmp(existingParams.Analysis, test);
                                    otherParams = existingParams(otherRows, :);
                                    project.models.(modelName).analysis.(analysisId).parameters = [otherParams; testParams];

                                    % Re-run and add analysis
                                    project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);

                                end
                            else
                                warning('No parameters table found. Continuing and storing given parameters.')
                                project = performAnalysis(project, parameterTable, modelName, {test}, analysisId);
                                project.models.(modelName).analysis.(analysisId).parameters = testParams;
                            end
                        end
                    end
                end
            else
                error('project.models.%s.analysis.%s must be a structure.', modelName, analysisId);
            end
        else
            error('project.models.%s.analysis.%s does not exist. Please choose an existing analysis or run singleModelAnalysis.', modelName, analysisId)
        end
    else
        error('There is no analysis field for project.models.%s. Please run singleModelAnalysis first.', modelName)
    end
else
    warning('project.models.%s was not found.', modelName);
end

% checking that project.models.(modelName).analysis.analysisId existe et
% que c'est bien une struct
% si oui, boucler sur les news analyses demandées, et regarder s'il y a
% déjà un champ au nom des analyses demandées
% si oui, regarder si parametres changent entre l'ancienne table et la
% nouvelle
% si oui, demander confirmation avant d'overwrite (display anciens params versus nouveaux), remplacer anciens params
% par nouveaux dans ceux stockés si confirmation, sinon abort.
% Si pas d'ancienne analysis correspondante, on update les infos params
% dans la table (on enlève anciennement stockées car pas utilisées, on met les nouvelles), et on run l'analyse, on stocke les results. 

% Think about checking if parameterTable isequal to the one present. If
% yes, ask for confirmation (displaying changes) before overwriting.

end

writeAnalysisReport(project, modelName, analysisId, varargin)

Generates a PDF report summarizing analysis results.

This function creates a PDF report for a single analysis run, including model characteristics, exchange fluxes, and pathway- or metabolite-level flux details. It requires that FBA and FVA have been performed on the model.

Input arguments:

Name Type Description Default
project struct

Project structure containing analysis results.

 required
modelName char

Name of the model to report on.

 required
analysisId char

Analysis ID (e.g. 'analysis_20240815_1430').

 required
path char

Output directory for the PDF file. Default: current folder. Passed as name-value pair.

 required
pathwaysOfInterest cell

Pathway names to include in the report. Default: empty. Passed as name-value pair.

 required
metsOfInterest cell

Metabolite IDs to include in the report. Default: empty. Passed as name-value pair.

 required

Output arguments:

Name Type Description

None. The PDF file is saved to the specified directory.

Examples:

% Generate a report with default settings
writeAnalysisReport(project, 'model1', 'analysis_20240815_1430', ...
    'path', './results/reports/');

% Include specific pathways and metabolites
writeAnalysisReport(project, 'model1', 'analysis_20240815_1430', ...
    'path', './results/reports/', ...
    'pathwaysOfInterest', {'Glycolysis', 'TCA cycle'}, ...
    'metsOfInterest', {'glc_D', 'o2', 'ac'});
Warning

FBA and FVA must have been performed on the model. If these analyses are not present, the function will error.

Note

The output file is named AnalysisReport.pdf. If the output directory does not exist, it is created automatically. Pathway names must match the subSystems field of the model. Metabolite IDs are matched using a prefix pattern (e.g. 'glc_D' matches 'glc_D[e]').

Source code in scr/classes/single_model_analysis/writeAnalysisReport.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
function [] = writeAnalysisReport(project, modelName, analysisId, varargin)
% Generates a PDF report summarizing analysis results.
%
% This function creates a PDF report for a single analysis run, including
% model characteristics, exchange fluxes, and pathway- or metabolite-level
% flux details. It requires that FBA and FVA have been performed on the
% model.
%
% Arguments:
%   project (struct): Project structure containing analysis results.
%   modelName (char): Name of the model to report on.
%   analysisId (char): Analysis ID (e.g. 'analysis_20240815_1430').
%   path (char): Output directory for the PDF file. Default: current
%       folder. Passed as name-value pair.
%   pathwaysOfInterest (cell): Pathway names to include in the report.
%       Default: empty. Passed as name-value pair.
%   metsOfInterest (cell): Metabolite IDs to include in the report.
%       Default: empty. Passed as name-value pair.
%
% Returns:
%   None. The PDF file is saved to the specified directory.
%
% Examples:
%   ```matlab
%   % Generate a report with default settings
%   writeAnalysisReport(project, 'model1', 'analysis_20240815_1430', ...
%       'path', './results/reports/');
%
%   % Include specific pathways and metabolites
%   writeAnalysisReport(project, 'model1', 'analysis_20240815_1430', ...
%       'path', './results/reports/', ...
%       'pathwaysOfInterest', {'Glycolysis', 'TCA cycle'}, ...
%       'metsOfInterest', {'glc_D', 'o2', 'ac'});
%   ```
%
% Warning:
%   FBA and FVA must have been performed on the model. If these analyses
%   are not present, the function will error.
%
% Note:
%   The output file is named <modelName>_AnalysisReport_<analysisId>.pdf.
%   If the output directory does not exist, it is created automatically.
%   Pathway names must match the subSystems field of the model. Metabolite
%   IDs are matched using a prefix pattern (e.g. 'glc_D' matches 'glc_D[e]').
%
p = inputParser;
addParameter(p,'path','');
addParameter(p,'pathwaysOfInterest',{});
addParameter(p,'metsOfInterest',{});
parse(p,varargin{:});

path = p.Results.path;
pathwaysOfInterest = p.Results.pathwaysOfInterest;
metsOfInterest = p.Results.metsOfInterest;

% Formatting inputs
path = char(path);

if isempty(pathwaysOfInterest)
    pathwaysOfInterest = {};
else
    pathwaysOfInterest = cellstr(string(pathwaysOfInterest));
end

if isempty(metsOfInterest)
    metsOfInterest = {};
else
    metsOfInterest = cellstr(string(metsOfInterest));
end

%% Checking existence of path folder
if ~exist(path, 'dir')   
    mkdir(path) % create folder if it doesn't exist
end

%% Checking existence of the required fields? --> FBA and FVA = minimal


%% Shortcut variables
% requirements
model = project.models.(modelName).model;
analysis = project.models.(modelName).analysis.(analysisId);
FBAflux = analysis.FBA.x;
GR = analysis.FBA.f;
minFlux = analysis.FVA.minMaxFluxes.minFlux;
maxFlux = analysis.FVA.minMaxFluxes.maxFlux;

%% Required folders for saving 
import mlreportgen.report.*
import mlreportgen.dom.*

%% Creating and opening the report
reportName = modelName + "_AnalysisReport_" + analysisId;
rpt = Report(reportName, "pdf");
open(rpt)
rpt.Layout.PageNumberFormat = 'n';

%% Title
dateStr = extractAfter(analysisId, "analysis_");
tp = TitlePage;
titleText = Text( ...
    modelName + " Model - Analysis Report " + ...
    string(datetime(dateStr, ...
        'InputFormat','yyyyMMdd_HHmm', ...
        'Format','dd MMMM yyyy, HH:mm')));
titleText.Bold = true;
titleText.FontSize = "24pt";
tp.Title = titleText;

add(rpt, tp);

%% Table of contents
toc = TableOfContents();
toc.Layout.PageNumberFormat = 'n';
add(rpt,toc);

add(rpt, PageBreak);

%% Define Table Style
tableStyle = {ColSep("solid"), ...
              RowSep("solid"), ...
              Border("solid")};

tableHeaderStyle = {BackgroundColor("lightblue"), ...
                      Bold(true), ...
                      mlreportgen.dom.Hyphenation(true)};

TableEntriesStyle = { ...
    FontSize("10pt"), ...
    mlreportgen.dom.InnerMargin("2pt","2pt","2pt","2pt"), ...
    mlreportgen.dom.WhiteSpace("preserve"), ...
    mlreportgen.dom.Hyphenation(false), ... % prevent from cutting words
    HAlign("center"), ...
    VAlign("middle") ...
    };

TableEntriesStyle2 = { ...
    FontSize("8pt"), ...
    mlreportgen.dom.InnerMargin("2pt","2pt","2pt","2pt"), ...
    mlreportgen.dom.WhiteSpace("preserve"), ...
    mlreportgen.dom.Hyphenation(true), ... % prevent from cutting words
    HAlign("center"), ...
    VAlign("middle") ...
    };

TableEntriesStyle3 = { ...
    FontSize("6pt"), ...
    mlreportgen.dom.InnerMargin("2pt","2pt","2pt","2pt"), ...
    mlreportgen.dom.WhiteSpace("preserve"), ...
    mlreportgen.dom.Hyphenation(false), ... % prevent from cutting words
    HAlign("left"), ...
    VAlign("middle") ...
    };

TableEntriesStyle4 = { ...
    FontSize("8pt"), ...
    mlreportgen.dom.InnerMargin("2pt","2pt","2pt","2pt"), ...
    mlreportgen.dom.WhiteSpace("preserve"), ...
    mlreportgen.dom.Hyphenation(true), ... % prevent from cutting words
    mlreportgen.dom.NumberFormat("%0.4f"), ...
    HAlign("left"), ...
    VAlign("middle") ...
    };

TableEntriesStyle5 = { ...
    FontSize("8pt"), ...
    mlreportgen.dom.InnerMargin("2pt","2pt","2pt","2pt"), ...
    mlreportgen.dom.WhiteSpace("preserve"), ...
    mlreportgen.dom.Hyphenation(" "), ... % prevent from cutting words
    mlreportgen.dom.NumberFormat("%0.4f"), ...
    HAlign("left"), ...
    VAlign("middle") ...
    };

lineBreak = Paragraph("");
lineBreak.Style = {OuterMargin("0pt", "0pt", "1pt", "1pt")}; % OuterMargin(left, right, top, bottom)

%HeadingStyle = {OuterMargin('0pt','0pt','0pt','0pt')};

%% Analysis Parameters
add(rpt, Heading(1, "Analysis Parameters"));
add(rpt, lineBreak);

parameters = analysis.parameters;

% Header
headerLabels = parameters.Properties.VariableNames;
% Body
tableBody = table2cell(parameters);

parametersTable = FormalTable(headerLabels, tableBody);

parametersTable.Style = [parametersTable.Style, tableStyle];
parametersTable.Header.Style = [parametersTable.Header.Style, tableHeaderStyle];

parametersTable.TableEntriesStyle = TableEntriesStyle;

parametersTable.Width = "100%";

add(rpt, parametersTable);

%% Model Characteristics
add(rpt, Heading(1, "Model Characteristics"))
add(rpt, lineBreak);

characteristics = struct(); 

% Consistency
nbConsistentRxns = fastcc(model, 1e-4, 0);
if length(nbConsistentRxns) == length(model.rxns)
    characteristics.modelConsistency = "YES";
else
    characteristics.modelConsistency = "NO";
end
characteristics.nbRxns = length(model.rxns);
characteristics.nbConsistentRxns = length(nbConsistentRxns);

% Nb of uptake rxns
[~, upt] = findExcRxns(model);
characteristics.nbUptakeRxns = sum(upt);

% GPR rules
characteristics.nbGPRrules = length(~cellfun(@isempty, model.grRules));

% Nb of core retained in the model
if isfield(project.models.(modelName), 'coreReactions')
    characteristics.nbCoreRxnsAfterDiscretization = length(project.models.(modelName).coreReactions);
    retainedRxns = model.rxns;
    coreRetained = intersect(retainedRxns, project.models.(modelName).coreReactions);
    characteristics.nbCoreRetained = length(coreRetained);
end

characteristicsT = struct2table(characteristics);

% Header
headerLabels = characteristicsT.Properties.VariableNames;
headerLabelsSpaced = regexprep(headerLabels, '([a-z])([A-Z])', '$1 $2');

% Body
tableBody = table2cell(characteristicsT);

characteristicsTable = FormalTable(headerLabelsSpaced, tableBody);

characteristicsTable.Style = [characteristicsTable.Style, tableStyle];
characteristicsTable.Header.Style = [characteristicsTable.Header.Style, tableHeaderStyle];

characteristicsTable.TableEntriesStyle = TableEntriesStyle2;

characteristicsTable.Width = "100%";

add(rpt, characteristicsTable);
add(rpt, PageBreak);

%% Medium
if isfield(project.models.(modelName).settings, 'medium')
    medium = project.models.(modelName).settings.medium.mediumComposition;
    idx = strcmp(parameters.Parameter, "modelReference");
    if any(idx)
        modelRef = parameters.Value{idx};
        if isequal(modelRef, "Recon3D") && all(ismember(["Mets_Recon3D", "ExRxns_Recon3D", "Concentration_uM"], medium.Properties.VariableNames))
            mediumShort = medium(:, ["Mets_Recon3D","ExRxns_Recon3D","Concentration_uM"]);
            if ismember("Flux_mmol_gDW_h", medium.Properties.VariableNames)
                mediumShort.Flux_mmol_gDW_h = medium.Flux_mmol_gDW_h;
            end
        elseif isequal(modelRef, "HumanGEM") && all(ismember(["Mets_HumanGEM", "ExRxns_HumanGEM", "Concentration_uM"], medium.Properties.VariableNames))
            mediumShort = medium(:, ["Mets_HumanGEM","ExRxns_HumanGEM","Concentration_uM"]);
            if ismember("Flux_mmol_gDW_h", medium.Properties.VariableNames)
                mediumShort.Flux_mmol_gDW_h = medium.Flux_mmol_gDW_h;
            end
        else
            mediumShort = '';
        end
    else
        mediumShort = '';
    end
    if ~isempty(mediumShort)
        add(rpt, Heading(1, "Medium Composition"))
        add(rpt, lineBreak);

        headerLabels = mediumShort.Properties.VariableNames;
        tableBody = table2cell(mediumShort);

        mediumTable = FormalTable(headerLabels, tableBody);

        mediumTable.Style = [mediumTable.Style, tableStyle];
        mediumTable.Header.Style = [mediumTable.Header.Style, tableHeaderStyle];

        mediumTable.TableEntriesStyle = TableEntriesStyle3;

        mediumTable.Width = "100%";

        add(rpt, mediumTable);
    end
end

%% Table of Exchangers
add(rpt, Heading(1,"Table of Exchangers"))
add(rpt, lineBreak);

p = Paragraph(sprintf("FBA for objective function: %.4f", GR));
p.Style = {FontSize("9pt")};   
add(rpt, p);
add(rpt, lineBreak);

% Exchangers
idxEX = startsWith(model.rxns, 'EX_');
EXrxns = model.rxns(idxEX);

allFluxEX = sortrows(table(EXrxns, FBAflux(idxEX), minFlux(idxEX), maxFlux(idxEX), ...
    FBAflux(idxEX)/GR, minFlux(idxEX)/GR, maxFlux(idxEX)/GR, ...
    'VariableNames', {'ReactionID', 'FBA', 'FVAmin', 'FVAmax', 'NormalizedFBA', 'NormalizedFVAmin', 'NormalizedFVAmax'}), 2);

% Header
headerLabels = allFluxEX.Properties.VariableNames;
% Body
tableBody = table2cell(allFluxEX);

FluxEXTable = FormalTable(headerLabels, tableBody);

FluxEXTable.Style = [FluxEXTable.Style, tableStyle];
FluxEXTable.Header.Style = [FluxEXTable.Header.Style, tableHeaderStyle];

FluxEXTable.TableEntriesStyle = TableEntriesStyle4;

FluxEXTable.Width = "100%";

add(rpt, FluxEXTable);
add(rpt, PageBreak);


%% Fluxes Per Pathway

if ~isempty(pathwaysOfInterest)

    add(rpt, Heading(1,"Fluxes Per Pathway"));
    add(rpt, lineBreak);

    nbPathwaysIn = 0;
    for i = 1:numel(pathwaysOfInterest)

        pathway = pathwaysOfInterest{i};

        if ismember(pathway, model.subSystems)
            nbPathwaysIn = nbPathwaysIn + 1;

            add(rpt, Heading(2, pathway));
            add(rpt, lineBreak);

            pathwayRxns = findRxnsFromSubSystem(model, pathway);
            idxPathway = find(ismember(model.rxns, pathwayRxns));
            pathwayFluxes = table(model.rxns(idxPathway), printRxnFormula(model, pathwayRxns, 0), ...
                FBAflux(idxPathway), minFlux(idxPathway), maxFlux(idxPathway), ...
                FBAflux(idxPathway)/GR, minFlux(idxPathway)/GR, maxFlux(idxPathway)/GR, ...
                'VariableNames', {'ReactionID', 'Formula', 'FBA', 'FVAmin', 'FVAmax', ...
                'NormalizedFBA', 'NormalizedFVAmin', 'NormalizedFVAmax'});

            % Pathway Flux Sum 
            % with FBA
            fluxSumFBA = computeFluxSum(model, analysis, 'FBA', 'pathway', pathway);
            % with sampling
            if isfield(analysis, 'sampling')
                fluxSumSampling = computeFluxSum(model, analysis, 'sampling', 'pathway', pathway);
            end

        else
            fprintf('%s is not part of the model.\n', pathway);
        end

        if exist('pathwayFluxes','var')
            % Header
            headerLabels = pathwayFluxes.Properties.VariableNames;
            % Body
            tableBody = table2cell(pathwayFluxes);

            pathwayFluxesTable = FormalTable(headerLabels, tableBody);

            pathwayFluxesTable.Style = [pathwayFluxesTable.Style, tableStyle];
            pathwayFluxesTable.Header.Style = [pathwayFluxesTable.Header.Style, tableHeaderStyle];

            pathwayFluxesTable.TableEntriesStyle = TableEntriesStyle5;

            pathwayFluxesTable.Width = "100%";

            add(rpt, Heading(3, "Fluxes"));
            add(rpt, lineBreak);
            add(rpt, pathwayFluxesTable);
            add(rpt, lineBreak);

            % Flux Sum with FBA
            add(rpt, Heading(3, "FBA Flux Sum"));
            add(rpt, lineBreak);
            % Header
            headerLabels = fluxSumFBA.Properties.VariableNames;
            % Body
            tableBody = table2cell(fluxSumFBA);

            pathwayFBAFluxSumTable = FormalTable(headerLabels, tableBody);

            pathwayFBAFluxSumTable.Style = [pathwayFBAFluxSumTable.Style, tableStyle];
            pathwayFBAFluxSumTable.Header.Style = [pathwayFBAFluxSumTable.Header.Style, tableHeaderStyle];

            pathwayFBAFluxSumTable.TableEntriesStyle = TableEntriesStyle5;

            pathwayFBAFluxSumTable.Width = "100%";

            add(rpt, pathwayFBAFluxSumTable);

            if exist('fluxSumSampling', 'var')
                add(rpt, Heading(3, "Sampling Flux Sum"));
                add(rpt, lineBreak);
                % Header
                headerLabels = fluxSumSampling.Properties.VariableNames;
                % Body
                tableBody = table2cell(fluxSumSampling);

                pathwaySamplingFluxSumTable = FormalTable(headerLabels, tableBody);

                pathwaySamplingFluxSumTable.Style = [pathwaySamplingFluxSumTable.Style, tableStyle];
                pathwaySamplingFluxSumTable.Header.Style = [pathwaySamplingFluxSumTable.Header.Style, tableHeaderStyle];

                pathwaySamplingFluxSumTable.TableEntriesStyle = TableEntriesStyle5;

                pathwaySamplingFluxSumTable.Width = "100%";

                add(rpt, pathwaySamplingFluxSumTable);
            end

        end
        add(rpt, PageBreak);
    end
    if nbPathwaysIn == 0
        fprintf("None of the given pathways were found in the model.\n Please check the subSystem field of your model.\n");
    end
end

%% Fluxes per metabolite
if ~isempty(metsOfInterest)

    add(rpt, Heading(1,"Fluxes Per Metabolite"));
    add(rpt, lineBreak);

    nbMetsIn = 0;
    for i = 1:numel(metsOfInterest)

        met = metsOfInterest{i};
        pattern = ['^', met, '\['];
        idx = ~cellfun('isempty', regexp(model.mets, pattern, 'once'));
        matchedMets = model.mets(idx);

        if ~isempty(matchedMets)

            nbMetsIn = nbMetsIn + 1;

            metName = string(model.metNames(find(strcmp(model.mets, matchedMets(1)))));
            add(rpt, Heading(2, sprintf('%s (%s)', metName, met)));
            add(rpt, lineBreak);

            [metsRxns, ~] = findRxnsFromMets(model, matchedMets);
            idxMetsRxns = find(ismember(model.rxns, metsRxns));
            %idxMets = find(ismember(model.mets, matchedMets)); % for Flux Sum

            % FBA/FVA for reactions associated with the metabolite
            metFluxes = table(model.rxns(idxMetsRxns), printRxnFormula(model, metsRxns, 0), FBAflux(idxMetsRxns),...
                minFlux(idxMetsRxns), maxFlux(idxMetsRxns), ...
                FBAflux(idxMetsRxns)/GR, minFlux(idxMetsRxns)/GR, maxFlux(idxMetsRxns)/GR, ...
                'VariableNames', {'ReactionID', 'Formula', 'FBA', 'FVAmin', 'FVAmax', 'NormalizedFBA', 'NormalizedFVAmin', 'NormalizedFVAmax'});

            % Metabolite Flux Sum 
            % with FBA
            fluxSumFBA = computeFluxSum(model, analysis, 'FBA', 'metabolite', met);
            % with sampling
            if isfield(analysis, 'sampling')
                fluxSumSampling = computeFluxSum(model, analysis, 'sampling', 'metabolite', met);
            end
        else
            fprintf('No metabolite corresponding to %s was found in the model.\n', met);
        end

        if exist('metFluxes','var')
            % Fluxes
            % Header
            headerLabels = metFluxes.Properties.VariableNames;
            % Body
            tableBody = table2cell(metFluxes);

            metFluxesTable = FormalTable(headerLabels, tableBody);

            metFluxesTable.Style = [metFluxesTable.Style, tableStyle];
            metFluxesTable.Header.Style = [metFluxesTable.Header.Style, tableHeaderStyle];

            metFluxesTable.TableEntriesStyle = TableEntriesStyle5;

            metFluxesTable.Width = "100%";

            add(rpt, Heading(3, "Fluxes"));
            add(rpt, lineBreak);
            add(rpt, metFluxesTable);
            add(rpt, lineBreak);

            % Flux Sum with FBA
            add(rpt, Heading(3, "FBA Flux Sum"));
            add(rpt, lineBreak);
            % Header
            headerLabels = fluxSumFBA.Properties.VariableNames;
            % Body
            tableBody = table2cell(fluxSumFBA);

            metFBAFluxSumTable = FormalTable(headerLabels, tableBody);

            metFBAFluxSumTable.Style = [metFBAFluxSumTable.Style, tableStyle];
            metFBAFluxSumTable.Header.Style = [metFBAFluxSumTable.Header.Style, tableHeaderStyle];

            metFBAFluxSumTable.TableEntriesStyle = TableEntriesStyle5;

            metFBAFluxSumTable.Width = "100%";

            add(rpt, metFBAFluxSumTable);

            if exist('fluxSumSampling', 'var')
                add(rpt, Heading(3, "Sampling Flux Sum"));
                add(rpt, lineBreak);
                % Header
                headerLabels = fluxSumSampling.Properties.VariableNames;
                % Body
                tableBody = table2cell(fluxSumSampling);

                metSamplingFluxSumTable = FormalTable(headerLabels, tableBody);

                metSamplingFluxSumTable.Style = [metSamplingFluxSumTable.Style, tableStyle];
                metSamplingFluxSumTable.Header.Style = [metSamplingFluxSumTable.Header.Style, tableHeaderStyle];

                metSamplingFluxSumTable.TableEntriesStyle = TableEntriesStyle5;

                metSamplingFluxSumTable.Width = "100%";

                add(rpt, metSamplingFluxSumTable);
            end

        end
        add(rpt, PageBreak);
    end
    if nbMetsIn == 0
        fprintf("None of the given metabolites were found in the model.\n Please check the mets field of your model.\n");
    end
end

close(rpt);
%% Needed to be integrated
% add(rpt, Heading(1,"Shadow Prices"))
% add(rpt, Heading(1,"Flux Sum from FBA"))
% add(rpt, Heading(1,"Flux Sum from Sampling"))
% 
% 
% 
% % Cr�e la figure
% f = figure;
% plot(1:10,(1:10).^2);
% title("Quadratic plot");
% 
% % Ajoute la figure directement au rapport via handle
% add(rpt, Figure(f));
% 
% T = table((1:5)', rand(5,1), 'VariableNames',{'ID','Value'});
% add(rpt, BaseTable(T))
% 
% close(rpt)
% rptview(rpt)

end

Model Comparison

chooseActiveAnalysis(project, modelList, analysisIDs={}, overwriteActive={'all'})

Designates which analysis run to use for model comparison.

This function must be run before modelComparison. Since multiple analysis runs can coexist on the same model, this function selects which one to use by copying it into an 'active' slot. This allows downstream comparison functions to access results without specifying the exact analysis ID for each model.

Input arguments:

Name Type Description Default
project struct

Project structure with completed single model analyses.

 required
modelList cell

Names of the models to define an active analysis for.

 required
analysisIDs cell

Analysis ID to set as active for each model, in the same order as modelList. If empty, the most recent analysis (by timestamp) is automatically selected for each model.

 required
overwriteActive cell

Fields to overwrite in the existing active slot. Use {'all'} (default) for full replacement, or specify individual fields (e.g. {'FBA'}) to selectively overwrite while preserving other results.

 required

Output arguments:

Name Type Description
project struct

Project with an active analysis defined for each model.
activeAnalysisTable table

Summary of the active analysis IDs used per model.

Examples:

% Use the most recent analysis for each model
[project, activeTable] = chooseActiveAnalysis(project, ...
    {'model1', 'model2'});

% Specify explicit analysis IDs
[project, activeTable] = chooseActiveAnalysis(project, ...
    {'model1', 'model2'}, ...
    {'analysis_20240815_1430', 'analysis_20240816_0900'});

% Selectively overwrite only FBA in the active slot
[project, activeTable] = chooseActiveAnalysis(project, ...
    {'model1'}, {'analysis_20240815_1430'}, {'FBA'});
Note

When overwriteActive is set to a specific field (e.g. {'FBA'}), the parameters table is automatically merged: rows corresponding to the overwritten analyses are replaced, while rows for other analyses are preserved.

Source code in scr/classes/model_comparison/functions/chooseActiveAnalysis.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
function [project, activeAnalysisTable] = chooseActiveAnalysis(project, modelList, analysisIDs, overwriteActive)
% Designates which analysis run to use for model comparison.
%
% This function must be run before modelComparison. Since multiple
% analysis runs can coexist on the same model, this function selects which
% one to use by copying it into an 'active' slot. This allows downstream
% comparison functions to access results without specifying the exact
% analysis ID for each model.
%
% Arguments:
%   project (struct): Project structure with completed single model
%       analyses.
%   modelList (cell): Names of the models to define an active analysis
%       for.
%   analysisIDs (cell): Analysis ID to set as active for each model, in
%       the same order as modelList. If empty, the most recent analysis
%       (by timestamp) is automatically selected for each model.
%   overwriteActive (cell): Fields to overwrite in the existing active
%       slot. Use {'all'} (default) for full replacement, or specify
%       individual fields (e.g. {'FBA'}) to selectively overwrite while
%       preserving other results.
%
% Returns:
%   project (struct): Project with an active analysis defined for each
%       model.
%   activeAnalysisTable (table): Summary of the active analysis IDs used
%       per model.
%
% Examples:
%   ```matlab
%   % Use the most recent analysis for each model
%   [project, activeTable] = chooseActiveAnalysis(project, ...
%       {'model1', 'model2'});
%
%   % Specify explicit analysis IDs
%   [project, activeTable] = chooseActiveAnalysis(project, ...
%       {'model1', 'model2'}, ...
%       {'analysis_20240815_1430', 'analysis_20240816_0900'});
%
%   % Selectively overwrite only FBA in the active slot
%   [project, activeTable] = chooseActiveAnalysis(project, ...
%       {'model1'}, {'analysis_20240815_1430'}, {'FBA'});
%   ```
%
% Note:
%   When overwriteActive is set to a specific field (e.g. {'FBA'}), the
%   parameters table is automatically merged: rows corresponding to the
%   overwritten analyses are replaced, while rows for other analyses are
%   preserved.

    arguments
        project
        modelList (1,:) cell
        analysisIDs (1,:) cell = {}
        overwriteActive (1,:) cell = {'all'}
    end

    % Convert inputs to string arrays for internal use
    modelList = string(modelList);
    overwriteActive = string(overwriteActive);
    if isempty(analysisIDs)
        analysisIDs = strings(1,0);
    else
        analysisIDs = string(analysisIDs);
    end

    % Validate overwriteActive
    if ismember("all", overwriteActive) && numel(overwriteActive) > 1
        error("overwriteActive cannot contain 'all' together with other field names.");
    end

    %% Check format of models/analysisIDs and defining active analyses
    if isempty(analysisIDs)
        checkProjectFormat(project, modelList)

        % Defining the most recent analysis as active in case no
        % analysisIDs were provided
        for m = 1:numel(modelList)
            modelShortcut = project.models.(modelList(m));

            if isfield(modelShortcut, 'analysis')
                analysisFields = string(fieldnames(modelShortcut.analysis));

                % Check that analysisFields is not empty
                if isempty(analysisFields)
                    error("There is no analysis entry for this model: '" + modelList(m) + "'. Run singleModelAnalysis function in order to create an analysis field for the specified models.");
                end

                % Search for slots created by singleModelAnalysis
                isA = startsWith(analysisFields, "analysis_");

                % Check that there is at least one analysis_ field
                if ~any(isA)
                    error("There is no analysis_ entry for this model: '" + modelList(m) + "'. Run singleModelAnalysis function in order to create an analysis field for the specified models.");
                end

                % Validate date format for analysis_ fields
                analysisFieldsA = analysisFields(isA);
                dateStrs = extractAfter(analysisFieldsA, "analysis_");
                isValidDate = false(size(dateStrs));

                for k = 1:numel(dateStrs)
                    try
                        datetime(dateStrs(k), 'InputFormat', 'yyyyMMdd_HHmm');
                        isValidDate(k) = true;
                    catch
                        % not a valid date — will be reported below
                    end
                end

                % Report invalid IDs
                if any(~isValidDate)
                    invalidIds = analysisFieldsA(~isValidDate);
                    fprintf("Warning: The following analysis IDs in model '%s' are not valid dates (yyyyMMdd_HHmm) and will be ignored: %s", ...
                        modelList(m), strjoin(invalidIds, ", "));
                end

                % Check that there is at least one valid date
                if ~any(isValidDate)
                    error("There is no valid analysis_ entry (yyyyMMdd_HHmm) for this model: '" + modelList(m) + "'. Run singleModelAnalysis function in order to create an analysis field for the specified models.");
                end

                % compute time difference between now and the time the
                % analysis were created (only for valid dates)
                timeDiff = NaN(size(analysisFields));
                validIdx = find(isA);
                validIdx = validIdx(isValidDate);
                timeDiff(validIdx) = minutes(datetime("now") - ...
                    datetime(dateStrs(isValidDate), 'InputFormat', 'yyyyMMdd_HHmm'));
                [~, idx] = min(timeDiff); % choose the most recently performed one

                analysisIDs(m) = analysisFields(idx);

            else
                error("There is no analysis field for this model: '" + modelList(m) + "'. Run singleModelAnalysis function in order to create an analysis field for the specified models.");
            end
        end

    else
        checkProjectFormat(project, modelList, analysisIDs)
    end

    %% Creating/updating active field
    for m = 1:numel(modelList)
        modelShortcut = project.models.(modelList(m));
        srcAnalysis = modelShortcut.analysis.(analysisIDs(m));
        srcFields = fieldnames(srcAnalysis);

        if ~isfield(modelShortcut.analysis, 'active') || isequal(overwriteActive, "all")
            % (Re)initializing active field — full replacement
            modelShortcut.analysis.active = struct();

            % Copy all fields from the selected analysis into active
            for k = 1:numel(srcFields)
                fieldName = srcFields{k};
                modelShortcut.analysis.active.(fieldName) = srcAnalysis.(fieldName);
                % Add analysisId to structs and tables only
                if isstruct(modelShortcut.analysis.active.(fieldName))
                    modelShortcut.analysis.active.(fieldName).analysisId = analysisIDs(m);
                end
            end

        else
            % Selective replacement — only overwrite specified fields
            for i = 1:numel(overwriteActive)
                fieldName = overwriteActive(i);

                % Check that the field exists in the source analysis
                if ~isfield(srcAnalysis, fieldName)
                    error("Field '%s' does not exist in analysis '%s' for model '%s'. Cannot overwrite.", ...
                        fieldName, analysisIDs(m), modelList(m));
                end

                % Report whether it's a new field or an overwrite
                if ~isfield(modelShortcut.analysis.active, fieldName)
                    warning("Field '%s' does not exist in active analysis for model '%s'. It will be added from analysis '%s'.", ...
                        fieldName, modelList(m), analysisIDs(m));
                else
                    fprintf("Overwriting '%s' in active analysis for model '%s' with analysis '%s'.", ...
                        fieldName, modelList(m), analysisIDs(m));
                end

                % Copy the field
                modelShortcut.analysis.active.(fieldName) = srcAnalysis.(fieldName);
                % Add analysisId to structs and tables only
                if isstruct(modelShortcut.analysis.active.(fieldName)) 
                    modelShortcut.analysis.active.(fieldName).analysisId = analysisIDs(m);
                end
            end

            % Merge parameters (if parameters was not explicitly overwritten)
            if ~ismember("parameters", overwriteActive)
                if isfield(srcAnalysis, 'parameters') && ...
                   isfield(modelShortcut.analysis.active, 'parameters')
                    % Determine which analyses are being replaced
                    parametersReplace = setdiff(overwriteActive, "parameters");
                    if ~isempty(parametersReplace)
                        oldParams = modelShortcut.analysis.active.parameters;
                        newParams = srcAnalysis.parameters;

                        % Keep rows from oldParams that are NOT being replaced
                        idxKeep = ~ismember(string(oldParams.Analysis), parametersReplace);
                        oldParams = oldParams(idxKeep, :);

                        % Take rows from newParams that ARE being replaced
                        idxAdd = ismember(string(newParams.Analysis), parametersReplace);
                        newParamsToAdd = newParams(idxAdd, :);

                        % Combine
                        modelShortcut.analysis.active.parameters = [oldParams; newParamsToAdd];
                    end
                end
            end
        end

        % Write back to project (MATLAB structs are pass-by-value)
        project.models.(modelList(m)) = modelShortcut;
    end

    activeAnalysisTable = getActiveAnalysisIDTable(project,modelList);

end

modelComparison(project, modelList, referenceModel, identifier=string(datetime('now', 'Format', '_yyyyMMdd_HHmmss')), analyses="structuralComparison")

Compares multiple models on structural, functional, and sampling levels.

This function runs a set of comparative analyses on the specified models. Three types of comparison are available: structural (presence or absence of reactions, metabolites, and genes), functional (FBA and FVA flux differences), and sampling (solution space comparison). Structural comparison is always run first as a prerequisite for the others.

Input arguments:

Name Type Description Default
project struct

Project structure with single model analyses completed and active analyses set via chooseActiveAnalysis.

 required
modelList string

Names of the models to compare.

 required
referenceModel string

Reference model used to compute relative reaction presence.

 required
identifier string

Postfix appended to the comparison name. Default: current timestamp.

 required
analyses string

Analyses to perform. Valid values: structuralComparison, functionalComparison, samplingComparison, IDAREoutput. Default: structuralComparison.

 required

Output arguments:

Name Type Description
project struct

Project with a comparisons field containing all results and plots.
comparisonName string

Name of the created comparison.

Examples:

% Run only the structural comparison (default)
[project, compName] = modelsComparison(project, ...
    ["model1", "model2"], "model1");

% Run structural and functional comparisons
[project, compName] = modelsComparison(project, ...
    ["model1", "model2"], "model1", "batchA", ...
    ["structuralComparison", "functionalComparison"]);

% Run all three comparisons
[project, compName] = modelsComparison(project, ...
    ["model1", "model2", "model3"], "model1", "fullRun", ...
    ["structuralComparison", "functionalComparison", "samplingComparison"]);
Note

The comparison name is built as model1_vs_model2_vs_...__identifier. Models are ordered by their appearance in project.models. If a comparison with the same name already exists and the structural analysis was already run, only the newly requested analyses are performed.

Warning

If a comparison with the same name exists but uses a different reference model, the user is prompted to confirm overwriting.

Source code in scr/classes/model_comparison/modelComparison.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
function [project, comparisonName] = modelsComparison(project, modelList, referenceModel, identifier, analyses)
% Compares multiple models on structural, functional, and sampling levels.
%
% This function runs a set of comparative analyses on the specified
% models. Three types of comparison are available: structural (presence
% or absence of reactions, metabolites, and genes), functional (FBA and
% FVA flux differences), and sampling (solution space comparison).
% Structural comparison is always run first as a prerequisite for the
% others.
%
% Arguments:
%   project (struct): Project structure with single model analyses
%       completed and active analyses set via chooseActiveAnalysis.
%   modelList (string): Names of the models to compare.
%   referenceModel (string): Reference model used to compute relative
%       reaction presence.
%   identifier (string): Postfix appended to the comparison name.
%       Default: current timestamp.
%   analyses (string): Analyses to perform. Valid values:
%       structuralComparison, functionalComparison, samplingComparison,
%       IDAREoutput. Default: structuralComparison.
%
% Returns:
%   project (struct): Project with a comparisons field containing all
%       results and plots.
%   comparisonName (string): Name of the created comparison.
%
% Examples:
%   ```matlab
%   % Run only the structural comparison (default)
%   [project, compName] = modelsComparison(project, ...
%       ["model1", "model2"], "model1");
%
%   % Run structural and functional comparisons
%   [project, compName] = modelsComparison(project, ...
%       ["model1", "model2"], "model1", "batchA", ...
%       ["structuralComparison", "functionalComparison"]);
%
%   % Run all three comparisons
%   [project, compName] = modelsComparison(project, ...
%       ["model1", "model2", "model3"], "model1", "fullRun", ...
%       ["structuralComparison", "functionalComparison", "samplingComparison"]);
%   ```
%
% Note:
%   The comparison name is built as model1_vs_model2_vs_...__identifier.
%   Models are ordered by their appearance in project.models. If a
%   comparison with the same name already exists and the structural
%   analysis was already run, only the newly requested analyses are
%   performed.
%
% Warning:
%   If a comparison with the same name exists but uses a different
%   reference model, the user is prompted to confirm overwriting.

arguments
    project        (1,1) struct
    modelList      (1,:) string
    referenceModel (1,1) string
    identifier     (1,1) string = string(datetime('now', 'Format', '_yyyyMMdd_HHmmss'))
    analyses       (1,:) string {mustBeMember(analyses, {'structuralComparison', 'functionalComparison', 'samplingComparison', 'IDAREoutput'})} = "structuralComparison"
end

%% Check project and models format for comparison

% Check first the reference model
checkProjectFormat(project, referenceModel)

% Check then models and analysis to compare
% active analysis of each model will be used for functionalComparison
checkProjectFormat(project, modelList, repmat({'active'}, 1, numel(modelList)))

%% Reorder models according to their order of appearance in project.models
order = string(fieldnames(project.models));
[~, idx] = ismember(modelList, order);
[~, sortIdx] = sort(idx);

modelListOrdered = modelList(sortIdx);

%% Give the comparison the name of all compared models associated with a given identifier
comparisonName = join(modelListOrdered, "_vs_") + "__" + identifier;

%% Create comparisons slot if not already existing
if ~isfield(project, 'comparisons')
    project.comparisons = struct();
elseif isfield(project, 'comparisons') && ~isstruct(project.comparisons)
    warning('project.comparisons exists but is not a struct (current class: %s).', ...
            class(project.comparisons));
    reply = '';
    while ~ismember(lower(strtrim(reply)), {'y', 'n', 'yes', 'no'})
        reply = input('Delete project.comparisons and reinitialize as an empty struct? [y/n] ', 's');
    end
    if startsWith(lower(strtrim(reply)), 'y')
        project.comparisons = struct();
        disp('project.comparisons reinitialized as an empty struct.');
    else
        error('Operation cancelled by user. project.comparisons left unchanged (class: %s).', ...
              class(project.comparisons));
    end
end

%% Comparison
% Structural comparison is needed for functional and sampling comparisons,
% therefore we check whether it was already run, to not run
% it again if it was already created

if isfield(project.comparisons, comparisonName)

    % TO ADD: checking of the format
    % If format correct, it means there is a field called
    % "referenceModel" which is part of project.comparisons

    if isequal(referenceModel, project.comparisons.(comparisonName).referenceModel)
        if isfield(project.comparisons.(comparisonName), 'structuralAnalysisStatus') && (project.comparisons.(comparisonName).structuralAnalysisStatus == 1)
            % Structural analysis already run
            disp("Structural comparison already run.");
            % Perform the other analyses if specified in analyses

            % Functional comparison
            if any(matches(analyses, "functionalComparison"))
                disp("Running functional comparison.");
                project.comparisons.(comparisonName).functionalComparison.plots = functionalComparison(project, comparisonName);
            end

            % Sampling comparison
            if any(matches(analyses, "samplingComparison"))
                disp("Running sampling comparison.");
                project = samplingComparison(project, comparisonName);
            end
        else
            % Structural analysis not run yet, mandatory
            project = runAllComparisons(project, modelListOrdered, referenceModel, comparisonName, analyses);
        end
    else
        % WARNING: field exists with a different referenceModel
        warning('Comparison "%s" already exists with referenceModel = "%s" (requested: "%s").', ...
                comparisonName, project.comparisons.(comparisonName).referenceModel, referenceModel);

        reply = '';
        while ~ismember(lower(strtrim(reply)), {'y', 'n', 'yes', 'no'})
            reply = input('Overwrite the existing comparison with the new reference model? [y/n] ', 's');
        end

        if startsWith(lower(strtrim(reply)), 'y')
            % YES: overwrite and rerun everything
            project = runAllComparisons(project, modelListOrdered, referenceModel, comparisonName, analyses);
        else
            % NO: stop, user must choose a new identifier
            error('Operation cancelled by user. Comparison "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
                  comparisonName, project.comparisons.(comparisonName).referenceModel);
        end
    end

else
    % Comparison does not exist yet: initialize and run everything
    project = runAllComparisons(project, modelListOrdered, referenceModel, comparisonName, analyses);
end

    % get the IDs for each of the analysis
    project.comparisons.(comparisonName).comparedAnalysisID = getActiveAnalysisIDTable(project,modelList);


end

showFigure(figHandle)

Creates a copy of a figure safely.

This function duplicates a figure given by its handle. It attempts several copy strategies in sequence: direct copyobj, children-only copyobj, and finally save-and-reopen via a temporary .fig file. This ensures compatibility with figures containing UIAxes, tables, clustergrams, and other complex objects.

Input arguments:

Name Type Description Default
figHandle figure

Handle of the figure to duplicate.

 required

Output arguments:

Name Type Description
newFig figure

Handle of the newly created figure copy.

Examples:

newFig = showFigure(plots.funct.import);
Note

The function tries three strategies in order of preference: 1. Direct copyobj of the entire figure 2. copyobj of children only into a new figure 3. Save to temporary .fig file and reopen

Source code in scr/classes/model_comparison/functions/showFigure.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
function newFig = showFigure(figHandle)
% Creates a copy of a figure safely.
%
% This function duplicates a figure given by its handle. It attempts
% several copy strategies in sequence: direct copyobj, children-only
% copyobj, and finally save-and-reopen via a temporary .fig file. This
% ensures compatibility with figures containing UIAxes, tables,
% clustergrams, and other complex objects.
%
% Arguments:
%   figHandle (figure): Handle of the figure to duplicate.
%
% Returns:
%   newFig (figure): Handle of the newly created figure copy.
%
% Examples:
%   ```matlab
%   newFig = showFigure(plots.funct.import);
%   ```
%
% Note:
%   The function tries three strategies in order of preference:
%   1. Direct copyobj of the entire figure
%   2. copyobj of children only into a new figure
%   3. Save to temporary .fig file and reopen

    try

        % tf = any(arrayfun(@(x) isa(x, 'matlab.ui.control.UIAxes'), figHandle.Children));
        % if tf 
        %     error("This is a ui figure it needs to be saved and reopened!")
        % end
        figure(copyobj(figHandle, groot));
        return;  % success, exit function

    catch
        try
            % if tf 
            %     error("This is a ui figure it needs to be saved and reopened!")
            % end
            children = get(figHandle, 'Children');  % get figure children
            newFig = figure;                        % create new figure
            copyobj(children, newFig);              % attempt copy
            set(newFig, 'Visible', 'on');           % ensure it shows
            return;  % success, exit function
        catch 
            % try
            %     % if tf 
            %     %     error("This is a ui figure it needs to be saved and reopened!")
            %     % end
                % Check input
                if ~ishandle(figHandle) || ~strcmp(get(figHandle,'Type'),'figure')
                    error('Input must be a valid figure handle.');
                end


                % Generate a temporary filename in the system temp folder
                tempFile = [tempname, '.fig'];

                try
                    % Save the figure to the temporary file
                    savefig(figHandle, tempFile);

                    % Open the figure as a new figure
                    newFig = openfig(tempFile, 'reuse'); % opens a new figure window

                    % Ensure it is visible
                    set(newFig, 'Visible', 'on');

                catch ME
                    % If any error occurs, delete the temp file and rethrow
                    if exist(tempFile, 'file')
                        delete(tempFile);
                    end
                    rethrow(ME);
                end

                % Delete the temporary file
                if exist(tempFile, 'file')
                    delete(tempFile);
                end
            % catch
            %     drawnow
            %     frame = getframe(figHandle);
            %     figure('Position', [10 10 1000 1000]) % [left bottom width height]
            %     imshow(frame.cdata)
            % end
        end
    end

end

getRxnIDs(project, referenceModel, pattern)

Finds reactions matching a pattern across model fields.

This function filters reactions for visualization using a regex pattern. The pattern can match subsystems, genes, reaction names, or metabolites. Multiple patterns can be combined with & (AND) or | (OR) operators.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelsComparison.

 required
referenceModel string

Reference model to search in.

 required
pattern string

Regex pattern(s) to match. Use & to require multiple patterns simultaneously (e.g. "^lac_. & ^Glycolysis.") or | to match any of several patterns (e.g. "^lac. | ^pyr.").

 required

Output arguments:

Name Type Description
rxnID cell

Reaction indices in the reference model for each pattern.
producing cell

Logical vector indicating whether matched metabolites are produced or consumed in each reaction, based on stoichiometry.
matchedAll cell

Actual strings matched by each pattern (reaction, gene, metabolite, or subsystem names).

Examples:

% Find reactions in Glycolysis involving lactate
[rxnID, producing, matched] = getRxnIDs(project, "model1", ...
    "^lac_.* & ^Glycolysis.*");

% Find reactions involving lactate or pyruvate
[rxnID, producing, matched] = getRxnIDs(project, "model1", ...
    "^lac.* | ^pyr.*");

% Find reactions in a single subsystem
[rxnID, producing, matched] = getRxnIDs(project, "model1", ...
    "^TCA cycle.*");
Note

If no match is found in model fields, the function searches the gene ID dictionary (settings.dico) to match gene names, symbols, or other identifiers. A pattern cannot contain both & and | operators.

Source code in scr/classes/model_comparison/functions/getRxnIDs.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
function [rxnID, producing,matchedAll] = getRxnIDs(project, referenceModel, pattern)
% Finds reactions matching a pattern across model fields.
%
% This function filters reactions for visualization using a regex pattern.
% The pattern can match subsystems, genes, reaction names, or metabolites.
% Multiple patterns can be combined with & (AND) or | (OR) operators.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelsComparison.
%   referenceModel (string): Reference model to search in.
%   pattern (string): Regex pattern(s) to match. Use & to require
%       multiple patterns simultaneously (e.g. "^lac_.* & ^Glycolysis.*")
%       or | to match any of several patterns (e.g. "^lac.* | ^pyr.*").
%
% Returns:
%   rxnID (cell): Reaction indices in the reference model for each
%       pattern.
%   producing (cell): Logical vector indicating whether matched
%       metabolites are produced or consumed in each reaction, based on
%       stoichiometry.
%   matchedAll (cell): Actual strings matched by each pattern (reaction,
%       gene, metabolite, or subsystem names).
%
% Examples:
%   ```matlab
%   % Find reactions in Glycolysis involving lactate
%   [rxnID, producing, matched] = getRxnIDs(project, "model1", ...
%       "^lac_.* & ^Glycolysis.*");
%
%   % Find reactions involving lactate or pyruvate
%   [rxnID, producing, matched] = getRxnIDs(project, "model1", ...
%       "^lac.* | ^pyr.*");
%
%   % Find reactions in a single subsystem
%   [rxnID, producing, matched] = getRxnIDs(project, "model1", ...
%       "^TCA cycle.*");
%   ```
%
% Note:
%   If no match is found in model fields, the function searches the
%   gene ID dictionary (settings.dico) to match gene names, symbols, or
%   other identifiers. A pattern cannot contain both & and | operators.

    arguments
        project (1,1) struct
        referenceModel (1,1) string
        pattern (1,:) string
    end

    model = project.models.(referenceModel).model;
    dict = project.models.(referenceModel).settings.dico;
    rxnID = {};
    producing = {};
    matchedAll = {};

    for n= 1:numel(pattern)

        singlePat = pattern(n);

        assert(~(contains(singlePat, "|") & contains(singlePat, "&")),...
               'The pattern can not contain both & + |, this function is only written for either or !')

        if contains(singlePat, "|")
            patCond =  strtrim(strsplit(singlePat, "|"));
        elseif contains(singlePat, "&")
            patCond =  strtrim(strsplit(singlePat, "&"));  
        else
            patCond = strtrim(singlePat);
        end

        % filter out all fields that can not be matched to pattern 
        fieldsToCheck = structfun(@(x) (ischar(x) || isstring(x) || iscell(x) ) ...
                                        && ((size(x,2) == 1 && size(x,1) == size(model.mets,1) || size(x,1) == size(model.rxns,1)|| size(x,1) == size(model.genes,1))),...
                                  model);
        namesToCheck = fieldnames(model);
        namesToNotCheck = namesToCheck(~fieldsToCheck);
        slotsCheck = rmfield(model,namesToNotCheck);

        [rxnsIDsAll,producingMetAll,matched] = cellfun(@(pat) findPatternsInStruct(pat,slotsCheck,model,dict),...
                                               patCond, 'UniformOutput', false);

        if contains(singlePat, "&")
            commonRxns = rxnsIDsAll{1};
            commonProd = producingMetAll{1};
            for i = 2:numel(rxnsIDsAll)
                [commonRxns, ia] = intersect(commonRxns, rxnsIDsAll{i});
                commonProd = commonProd(ia);
            end
            resultRxns = commonRxns;
            resultProd = commonProd;

        elseif contains(singlePat, "|")
            allRxns = vertcat(rxnsIDsAll{:});
            allProd = vertcat(producingMetAll{:});
            [resultRxns, ia] = unique(allRxns, 'stable');
            resultProd = allProd;
        else
            resultRxns = rxnsIDsAll{:,:};
            resultProd = producingMetAll{:,:};

        end

        rxnID{n} = resultRxns;
        producing{n} = resultProd;
        matchedAll{n} = string(vertcat(matched{:}));
        assert(~isempty(matchedAll{n}), 'No reactions were found! Check your pattern for spelling mistakes!')
    end
end

visDiffRxnSetActivityFBA(project, compName, rxnSets, rxnSetLabels, referenceModel)

Generates a heatmap showing how active the defined reactions are, by computing the sum of all FBA flux values for each reaction set.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSets cell

Sets of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabels string

Labels for each reaction set, displayed in the figure. Must be the same length as rxnSets.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required

Output arguments:

Name Type Description
fluxSet struct

Flux sum data used to generate the heatmap.
figs figure

Figure object generated and displayed by the function.

Examples:

% Retrieve reactions for two sets, then visualize
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
[fluxSet, figs] = visDiffRxnSetActivityFBA(project, compName, ...
    rxnsMetId, ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visDiffRxnSetActivityFBA.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
function [fluxSet,figs] = visDiffRxnSetActivityFBA(project, compName, rxnSets, rxnSetLabels, referenceModel)
% Generates a heatmap showing how active the defined reactions are, by
% computing the sum of all FBA flux values for each reaction set.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSets (cell): Sets of reaction indices to visualize, typically
%       obtained via getRxnIDs.
%   rxnSetLabels (string): Labels for each reaction set, displayed in
%       the figure. Must be the same length as rxnSets.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%
% Returns:
%   fluxSet (struct): Flux sum data used to generate the heatmap.
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Retrieve reactions for two sets, then visualize
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
%   [fluxSet, figs] = visDiffRxnSetActivityFBA(project, compName, ...
%       rxnsMetId, ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSets (1,:) cell
    rxnSetLabels (1,:) string
    referenceModel (1,1) string
end

%% check project format

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName),'functionalComparison')
    error('The comparison object you specified does not entail a functionalComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["functionalComparison"]).')
end
if ~isfield(project.comparisons.(compName).functionalComparison, 'orderedFba')
    error('The comparisons object does not entail a valid functionalComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["functionalComparison"])')
end

% rxnSets and rxSetLabels need to be the same length
if length(rxnSets) ~= length(rxnSetLabels)
    error('The rxnSets and the rxnSetLabels you specified do not have the same length, you need to give as many Lables as rxnID sets. ')
end

% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%

[fluxSet,figs] = visualizeFluxsum(project, compName,[],rxnSets,rxnSetLabels,"heatmap",...
                                       false,false,"orderedFba","reactions",referenceModel,'on')

end

visDiffRxnSetActivitySampling(project, compName, rxnSets, rxnSetLabels, referenceModel, metIdx)

Generates a heatmap showing how active the defined reactions are across sampling solutions. For each reaction set, the sum of flux values is computed per sample, then averaged over all samples to produce one value per reaction set in the heatmap.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSets cell

Sets of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabels string

Labels for each reaction set, displayed in the figure. Must be the same length as rxnSets.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required
metIdx

Optional metabolite indices for filtering.

 required

Output arguments:

Name Type Description
fluxSet struct

Flux sum data used to generate the heatmap.
figs figure

Figure object generated and displayed by the function.

Examples:

% Retrieve reactions for two sets, then visualize
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
[fluxSet, figs] = visDiffRxnSetActivitySampling(project, ...
    compName, rxnsMetId, ...
    ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visDiffRxnSetActivitySampling.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
function [fluxSet, figs] = visDiffRxnSetActivitySampling(project, compName, rxnSets, rxnSetLabels, referenceModel, metIdx)
% Generates a heatmap showing how active the defined reactions are
% across sampling solutions. For each reaction set, the sum of flux
% values is computed per sample, then averaged over all samples to
% produce one value per reaction set in the heatmap.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSets (cell): Sets of reaction indices to visualize, typically
%       obtained via getRxnIDs.
%   rxnSetLabels (string): Labels for each reaction set, displayed in
%       the figure. Must be the same length as rxnSets.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%   metIdx: Optional metabolite indices for filtering.
%
% Returns:
%   fluxSet (struct): Flux sum data used to generate the heatmap.
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Retrieve reactions for two sets, then visualize
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
%   [fluxSet, figs] = visDiffRxnSetActivitySampling(project, ...
%       compName, rxnsMetId, ...
%       ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSets (1,:) cell
    rxnSetLabels (1,:) string
    referenceModel (1,1) string
end

% check project format: 

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName),'samplingComparison')
    error('The comparison object you specified does not entail a samplingComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"]).')
end
if ~isfield(project.comparisons.(compName).samplingComparison, 'orderedSamples')
    error('The comparisons object does not entail a valid samplingComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"])')
end

% rxnSets and rxSetLabels need to be the same length
if length(rxnSets) ~= length(rxnSetLabels)
    error('The rxnSets and the rxnSetLabels you specified do not have the same length, you need to give as many Lables as rxnID sets. ')
end

% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%

[fluxSet,figs] = visualizeFluxsum(project, compName,metIdx,rxnSets,rxnSetLabels,"heatmap",...
                                       false,false,"orderedllSamples","reactions",referenceModel,'on')

end

visDiffMetSetUsageFBA(project, compName, rxnSets, rxnSetLabels, referenceModel)

Heatmap of metabolite usage across reaction sets based on FBA.

Generates a heatmap showing how much each metabolite participating in the defined reaction sets is used, based on FBA flux values. Usage is defined as the flux sum of all reactions that consume the metabolite within each reaction set.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSets cell

Sets of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabels string

Labels for each reaction set, displayed in the figure. Must be the same length as rxnSets.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required

Output arguments:

Name Type Description
figs figure

Figure object generated and displayed by the function.

Examples:

% Retrieve reactions for two sets, then visualize
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
figs = visDiffMetSetUsageFBA(project, compName, rxnsMetId, ...
    ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visDiffMetSetUsageFBA.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
function [figs] = visDiffMetSetUsageFBA(project, compName, rxnSets, rxnSetLabels, referenceModel)
% Heatmap of metabolite usage across reaction sets based on FBA.
%
% Generates a heatmap showing how much each metabolite participating in
% the defined reaction sets is used, based on FBA flux values. Usage is
% defined as the flux sum of all reactions that consume the metabolite
% within each reaction set.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSets (cell): Sets of reaction indices to visualize, typically
%       obtained via getRxnIDs.
%   rxnSetLabels (string): Labels for each reaction set, displayed in
%       the figure. Must be the same length as rxnSets.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%
% Returns:
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Retrieve reactions for two sets, then visualize
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
%   figs = visDiffMetSetUsageFBA(project, compName, rxnsMetId, ...
%       ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSets (1,:) cell
    rxnSetLabels (1,:) string
    referenceModel (1,1) string
end

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName),'functionalComparison')
    error('The comparison object you specified does not entail a functionalComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["functionalComparison"]).')
end
if ~isfield(project.comparisons.(compName).functionalComparison, 'orderedFba')
    error('The comparisons object does not entail a valid functionalComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["functionalComparison"])')
end

% rxnSets and rxSetLabels need to be the same length
if length(rxnSets) ~= length(rxnSetLabels)
    error('The rxnSets and the rxnSetLabels you specified do not have the same length, you need to give as many Lables as rxnID sets. ')
end


% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%

[fluxSet,figs] = visualizeFluxsum(project, compName,[],rxnSets,rxnSetLabels,"heatmap",...
                                       false,false,"orderedFba","outgoing",referenceModel,'on')

end

visDiffMetSetUsageSampling(project, compName, rxnSets, rxnSetLabels, referenceModel)

Generates a heatmap showing how much each metabolite participating in the defined reaction sets is used, based on sampling solutions. An average is computed over all samples to produce one value per reaction set in the heatmap. Usage is defined as the flux sum of all reactions that consume the metabolite within each reaction set.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSets cell

Sets of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabels string

Labels for each reaction set, displayed in the figure. Must be the same length as rxnSets.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required

Output arguments:

Name Type Description
figs figure

Figure object generated and displayed by the function.

Examples:

% Retrieve reactions for two sets, then visualize
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
figs = visDiffMetSetUsageSampling(project, compName, rxnsMetId, ...
    ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visDiffMetSetUsageSampling.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
function [figs] = visDiffMetSetUsageSampling(project, compName, rxnSets, rxnSetLabels, referenceModel)
% Generates a heatmap showing how much each metabolite participating in
% the defined reaction sets is used, based on sampling solutions. An
% average is computed over all samples to produce one value per reaction
% set in the heatmap. Usage is defined as the flux sum of all reactions
% that consume the metabolite within each reaction set.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSets (cell): Sets of reaction indices to visualize, typically
%       obtained via getRxnIDs.
%   rxnSetLabels (string): Labels for each reaction set, displayed in
%       the figure. Must be the same length as rxnSets.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%
% Returns:
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Retrieve reactions for two sets, then visualize
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
%   figs = visDiffMetSetUsageSampling(project, compName, rxnsMetId, ...
%       ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSets (1,:) cell
    rxnSetLabels (1,:) string
    referenceModel (1,1) string
end

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName),'samplingComparison')
    error('The comparison object you specified does not entail a samplingComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"]).')
end
if ~isfield(project.comparisons.(compName).samplingComparison, 'orderedSamples')
    error('The comparisons object does not entail a valid samplingComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"])')
end

% rxnSets and rxSetLabels need to be the same length
if length(rxnSets) ~= length(rxnSetLabels)
    error('The rxnSets and the rxnSetLabels you specified do not have the same length, you need to give as many Lables as rxnID sets. ')
end


% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%

[~,figs] = visualizeFluxsum(project, compName,[],rxnSets,rxnSetLabels,"heatmap",...
                                       false,false,"orderedSamples","outgoing",referenceModel,'on')

end

visRxnSetVariability(project, compName, rxnSets, rxnSetLabels, referenceModel)

Generates a heatmap showing how active the defined reactions are across sampling solutions. For each reaction set, the sum of flux values is computed per sample, then averaged over all samples to display one value per model and reaction set in the heatmap.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSets cell

Sets of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabels string

Labels for each reaction set, displayed in the figure. Must be the same length as rxnSets.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required

Output arguments:

Name Type Description
fluxSet struct

Flux sum data used to generate the heatmap.
figs figure

Figure object generated and displayed by the function.

Examples:

% Retrieve reactions for two sets, then visualize
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
[fluxSet, figs] = visDiffRxnSetActivitySampling(project, ...
    compName, rxnsMetId, ...
    ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visRxnSetVariability.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
function [fluxSet, figs] = visDiffRxnSetActivitySampling(project, compName, rxnSets, rxnSetLabels, referenceModel)
% Generates a heatmap showing how active the defined reactions are
% across sampling solutions. For each reaction set, the sum of flux
% values is computed per sample, then averaged over all samples to
% display one value per model and reaction set in the heatmap.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSets (cell): Sets of reaction indices to visualize, typically
%       obtained via getRxnIDs.
%   rxnSetLabels (string): Labels for each reaction set, displayed in
%       the figure. Must be the same length as rxnSets.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%
% Returns:
%   fluxSet (struct): Flux sum data used to generate the heatmap.
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Retrieve reactions for two sets, then visualize
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
%   [fluxSet, figs] = visDiffRxnSetActivitySampling(project, ...
%       compName, rxnsMetId, ...
%       ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSets (1,:) cell
    rxnSetLabels (1,:) string
    referenceModel (1,1) string
end

% check project format: 

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName),'samplingComparison')
    error('The comparison object you specified does not entail a samplingComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"]).')
end
if ~isfield(project.comparisons.(compName).samplingComparison, 'orderedSamples')
    error('The comparisons object does not entail a valid samplingComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"])')
end

% rxnSets and rxSetLabels need to be the same length
if length(rxnSets) ~= length(rxnSetLabels)
    error('The rxnSets and the rxnSetLabels you specified do not have the same length, you need to give as many Lables as rxnID sets. ')
end

% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%

[fluxSet,figs] = visualizeFluxsum(project, compName,[],rxnSets,rxnSetLabels,"heatmapSample",...
                                       false,false,"orderedSamples","reactions",referenceModel,'on')

end

visMetSetVariability(project, compName, rxnSets, rxnSetLabels, referenceModel)

Generates a heatmap showing how much each metabolite participating in the defined reaction sets is used, based on all sampling solutions. The resulting heatmap contains one value per metabolite set per sample. Usage is defined as the flux sum of all reactions that consume the metabolite within each reaction set.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSets cell

Sets of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabels string

Labels for each reaction set, displayed in the figure. Must be the same length as rxnSets.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required

Output arguments:

Name Type Description
figs figure

Figure object generated and displayed by the function.

Examples:

% Retrieve reactions for two sets, then visualize
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
figs = visMetSetVariability(project, compName, rxnsMetId, ...
    ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visMetSetVariability.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
function [figs] = visMetSetVariability(project, compName, rxnSets, rxnSetLabels, referenceModel)
% Generates a heatmap showing how much each metabolite participating in
% the defined reaction sets is used, based on all sampling solutions.
% The resulting heatmap contains one value per metabolite set per
% sample. Usage is defined as the flux sum of all reactions that
% consume the metabolite within each reaction set.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSets (cell): Sets of reaction indices to visualize, typically
%       obtained via getRxnIDs.
%   rxnSetLabels (string): Labels for each reaction set, displayed in
%       the figure. Must be the same length as rxnSets.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%
% Returns:
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Retrieve reactions for two sets, then visualize
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, ["Pentose.* & g6p.*"; "Glycolysis.*"]);
%   figs = visMetSetVariability(project, compName, rxnsMetId, ...
%       ["Pentose.* & g6p.*"; "Glycolysis.*"], referenceModel);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSets (1,:) cell
    rxnSetLabels (1,:) string
    referenceModel (1,1) string
end

% check project format

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName),'samplingComparison')
    error('The comparison object you specified does not entail a samplingComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"]).')
end
if ~isfield(project.comparisons.(compName).samplingComparison, 'orderedSamples')
    error('The comparisons object does not entail a valid samplingComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"])')
end

% rxnSets and rxSetLabels need to be the same length
if length(rxnSets) ~= length(rxnSetLabels)
    error('The rxnSets and the rxnSetLabels you specified do not have the same length, you need to give as many Lables as rxnID sets. ')
end


% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%

[~,figs] = visualizeFluxsum(project, compName,[],rxnSets,rxnSetLabels,"heatmapSample",...
                                       false,false,"orderedSamples","outgoing",referenceModel,'on')

end

visSingleRxnSamplingDistribution(project, compName, rxnSet, rxnSetLabel, referenceModel, addKLDValues=false)

Generates violin plots showing the distribution of flux values across all sampling solutions for each reaction in the defined set. One violin plot is produced per reaction index in rxnSet.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSet cell

Set of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabel string

Label displayed for this reaction set in the figure.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required
addKLDValues logical

If true, display the Kullback-Leibler divergence significance between sampling distributions on the violin plots. Default: false.

 required

Output arguments:

Name Type Description
fluxSet struct

Flux data used to generate the violin plots.
figs figure

Figure object generated and displayed by the function.

Examples:

% Visualize sampling distribution for a set of reactions
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, "Glycolysis.*");
[fluxSet, figs] = visSingleRxnSamplingDistribution(project, ...
    compName, rxnsMetId, "Glycolysis", referenceModel);

% Include KLD significance values
[fluxSet, figs] = visSingleRxnSamplingDistribution(project, ...
    compName, rxnsMetId, "Glycolysis", referenceModel, true);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visSingleRxnSamplingDistribution.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
function [fluxSet, figs] = visSingleRxnSamplingDistribution(project, compName, rxnSet, rxnSetLabel, referenceModel, addKLDValues)
% Generates violin plots showing the distribution of flux values across
% all sampling solutions for each reaction in the defined set. One
% violin plot is produced per reaction index in rxnSet.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSet (cell): Set of reaction indices to visualize, typically
%       obtained via getRxnIDs.
%   rxnSetLabel (string): Label displayed for this reaction set in the
%       figure.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%   addKLDValues (logical): If true, display the Kullback-Leibler
%       divergence significance between sampling distributions on the
%       violin plots. Default: false.
%
% Returns:
%   fluxSet (struct): Flux data used to generate the violin plots.
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Visualize sampling distribution for a set of reactions
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, "Glycolysis.*");
%   [fluxSet, figs] = visSingleRxnSamplingDistribution(project, ...
%       compName, rxnsMetId, "Glycolysis", referenceModel);
%
%   % Include KLD significance values
%   [fluxSet, figs] = visSingleRxnSamplingDistribution(project, ...
%       compName, rxnsMetId, "Glycolysis", referenceModel, true);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSet (1,:) cell
    rxnSetLabel (1,:) string
    referenceModel (1,1) string
    addKLDValues (1,1) logical = false
end

% check project format: 

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName), 'samplingComparison')
    error('The comparison object you specified does not entail a samplingComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"]).')
end
if ~isfield(project.comparisons.(compName).samplingComparison, 'orderedSamples')
    error('The comparisons object does not entail a valid samplingComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"])')
end

% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%

%
rxnLabel    = matlab.lang.makeValidName(rxnSetLabel);  % used as field name

[fluxSet, figs] = visualizeFlux(project, compName, rxnSet, rxnLabel, "orderedSamples", "all", 'on', addKLDValues)

end

visSingleMetSamplingDistribution(project, compName, rxnSet, rxnSetLabel, referenceModel)

Generates a violin plot showing the flux sum value distribution for all metabolites participating in the defined reaction set. Each dot in the violin plot represents one sampling solution.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison. Both must have been run before calling this function.

 required
compName string

Name of the comparison to visualize. Available comparisons can be listed with project.comparisons.

 required
rxnSet cell

A single set of reaction indices to visualize, typically obtained via getRxnIDs.

 required
rxnSetLabel string

Label displayed for this reaction set in the plot.

 required
referenceModel string

Name of the reference model. Must match the reference model used in the specified comparison and when retrieving reaction IDs with getRxnIDs.

 required

Output arguments:

Name Type Description
fluxSet struct

Flux sum data used to generate the violin plot.
figs figure

Figure object generated and displayed by the function.

Examples:

% Retrieve reactions for one set, then visualize
[rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
    referenceModel, "Glycolysis.*");
[fluxSet, figs] = visSingleMetSamplingDistribution(project, ...
    compName, rxnsMetId, "Glycolysis", referenceModel);
Warning

The reference model must be the same one used when calling getRxnIDs to retrieve the reaction indices. Using a different reference model will cause index misalignment.

Source code in scr/classes/model_comparison/functions/visSingleMetSamplingDistribution.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
function [fluxSet, figs] = visSingleMetSamplingDistribution(project, compName, rxnSet, rxnSetLabel, referenceModel)
% Generates a violin plot showing the flux sum value distribution for
% all metabolites participating in the defined reaction set. Each dot
% in the violin plot represents one sampling solution.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison. Both must have been run before calling this
%       function.
%   compName (string): Name of the comparison to visualize. Available
%       comparisons can be listed with project.comparisons.
%   rxnSet (cell): A single set of reaction indices to visualize,
%       typically obtained via getRxnIDs.
%   rxnSetLabel (string): Label displayed for this reaction set in the
%       plot.
%   referenceModel (string): Name of the reference model. Must match
%       the reference model used in the specified comparison and when
%       retrieving reaction IDs with getRxnIDs.
%
% Returns:
%   fluxSet (struct): Flux sum data used to generate the violin plot.
%   figs (figure): Figure object generated and displayed by the
%       function.
%
% Examples:
%   ```matlab
%   % Retrieve reactions for one set, then visualize
%   [rxnsMetId, producingMet, matched] = getRxnIDs(project, ...
%       referenceModel, "Glycolysis.*");
%   [fluxSet, figs] = visSingleMetSamplingDistribution(project, ...
%       compName, rxnsMetId, "Glycolysis", referenceModel);
%   ```
%
% Warning:
%   The reference model must be the same one used when calling
%   getRxnIDs to retrieve the reaction indices. Using a different
%   reference model will cause index misalignment.

arguments
    project 
    compName (1,1) string
    rxnSet (1,1) cell
    rxnSetLabel (1,1) string
    referenceModel (1,1) string
end

% check project format: 

if ~isfield(project, 'comparisons')
    error('Project object does not contain a comparison object. After running the singleModels analysis you still have to run modelsComparison before beeing able to use this function!')
end
if ~isfield(project.comparisons, compName)
    error('The comparison name you gave as an input is not available in the object. Check your spelling!')
end
if ~isfield(project.comparisons.(compName),'samplingComparison')
    error('The comparison object you specified does not entail a samplingComparison. Run modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"]).')
end
if ~isfield(project.comparisons.(compName).samplingComparison, 'orderedSamples')
    error('The comparisons object does not entail a valid samplingComparison object rerun the functionalComparison: modelsComparison(project, modelsToCompare, referenceModel, compID, ["samplingComparison"])')
end


% check that the reference model used when retrieving the rxnIDs is the
% same as in the comparison, otherwise the ids are wrongly assigned

if project.comparisons.(compName).referenceModel ~= referenceModel
    error('The reference Model you gave as input (%s) is not the same reference Model defined in the specified comparison object (%s). Check the spelling! The reference model when retrieving the rxnIDS with getRxnIDs needs to be the same as the one specified here! Otherwise the idx are misaligned.  "%s" keeps referenceModel = "%s".Choose a different identifier to create a new comparison.', ...
          referenceModel, project.comparisons.(compName).referenceModel);
end

%
rxnLabel    = matlab.lang.makeValidName(rxnSetLabel);  % used as field name

[fluxSet,figs] = visualizeFluxsum(project, compName, [], rxnSet, rxnLabel, "violin", ...
                                       false, false, "orderedSamples", "outgoing", referenceModel, 'on')

end

visSingleRxnFBA(project, comparisonName, idxToVis, FVA=false, thresholdFlux="none", titlePlots="", visiblePlots="on")

Visualizes FBA and FVA values for selected reactions across models in a comparison. By default, only FBA values are shown as a grouped horizontal bar plot. Optionally, FVA boundaries can be displayed as grey boxes around the FBA dots. Reactions are split into a high-flux and low-flux panel using 1D k-means clustering for readability.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison.

 required
comparisonName string

Name of the comparison to visualize.

 required
idxToVis cell

Indices of the reactions to display in the reference model, typically obtained via getRxnIDs.

 required
options.FVA logical

If true, display FVA boundaries as grey boxes around the FBA dots. Default: false.

 required
options.thresholdFlux string

Flux filtering mode. Valid values: "lower" (positive flux only), "upper" (negative flux only), "none" (non-zero flux only), "all" (include zero-flux reactions). Default: "none".

 required
options.titlePlots string

Custom title for the plots. Used when thresholdFlux is "none" or "all". Default: "".

 required
options.visiblePlots string

Figure visibility, "on" or "off". Default: "on".

 required

Output arguments:

Name Type Description
fig figure

Figure object containing the bar plot and optional table with reaction formulas, medium constraints, model mappings, and gene-protein-reaction rules.

Examples:

% Visualize FBA values for selected reactions
fig = visSingleRxnFBA(project, compName, rxnIDs);

% Include FVA boundaries and show only positive flux reactions
fig = visSingleRxnFBA(project, compName, rxnIDs, ...
    struct('FVA', true, 'thresholdFlux', "lower"));

% Show all reactions including zero-flux ones
fig = visSingleRxnFBA(project, compName, rxnIDs, ...
    struct('thresholdFlux', "all", 'titlePlots', "All reactions"));
Note

When all models share the same medium composition, a table is displayed alongside the plot showing reaction formulas, medium constraints, model-specific reaction mappings, and symbol GPR rules.

Source code in scr/classes/model_comparison/functions/visSingleRxnFBA.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
function fig = visSingleRxnFBA(project, comparisonName, idxToVis, options)
% Visualizes FBA and FVA values for selected reactions across models in
% a comparison. By default, only FBA values are shown as a grouped
% horizontal bar plot. Optionally, FVA boundaries can be displayed as
% grey boxes around the FBA dots. Reactions are split into a high-flux
% and low-flux panel using 1D k-means clustering for readability.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison.
%   comparisonName (string): Name of the comparison to visualize.
%   idxToVis (cell): Indices of the reactions to display in the
%       reference model, typically obtained via getRxnIDs.
%   options.FVA (logical): If true, display FVA boundaries as grey
%       boxes around the FBA dots. Default: false.
%   options.thresholdFlux (string): Flux filtering mode. Valid values:
%       "lower" (positive flux only), "upper" (negative flux only),
%       "none" (non-zero flux only), "all" (include zero-flux
%       reactions). Default: "none".
%   options.titlePlots (string): Custom title for the plots. Used when
%       thresholdFlux is "none" or "all". Default: "".
%   options.visiblePlots (string): Figure visibility, "on" or "off".
%       Default: "on".
%
% Returns:
%   fig (figure): Figure object containing the bar plot and optional
%       table with reaction formulas, medium constraints, model
%       mappings, and gene-protein-reaction rules.
%
% Examples:
%   ```matlab
%   % Visualize FBA values for selected reactions
%   fig = visSingleRxnFBA(project, compName, rxnIDs);
%
%   % Include FVA boundaries and show only positive flux reactions
%   fig = visSingleRxnFBA(project, compName, rxnIDs, ...
%       struct('FVA', true, 'thresholdFlux', "lower"));
%
%   % Show all reactions including zero-flux ones
%   fig = visSingleRxnFBA(project, compName, rxnIDs, ...
%       struct('thresholdFlux', "all", 'titlePlots', "All reactions"));
%   ```
%
% Note:
%   When all models share the same medium composition, a table is
%   displayed alongside the plot showing reaction formulas, medium
%   constraints, model-specific reaction mappings, and symbol GPR rules.

    arguments
        project
        comparisonName (1,1) string
        idxToVis  (1,1) cell
        options.FVA  (1,1) logical = false
        %options.reducedCost (1,1) logical = false
        options.thresholdFlux (1,1) string {mustBeMember(options.thresholdFlux, ["lower", "upper", "none", "all"])} = "none" 
        options.titlePlots = ""
        options.visiblePlots = "on"
    end

    idxToVis = idxToVis{1};
    modelList = project.comparisons.(comparisonName).modelNames;
    referenceModel = project.comparisons.(comparisonName).referenceModel;

    replacementValue = "analysis.active.FBA.v"; % get the fba solution values
    orderedFbaMatrix = getOrderedFeatureMatrix(project, modelList, referenceModel, "rxns", replacementValue);
    replacementValue = "mappedDiscretizedRxns"; % get the fba solution values
    orderedMappingRxnMatrix = getOrderedFeatureMatrix(project, modelList, referenceModel, "rxns", replacementValue);

    if options.thresholdFlux == "upper"
        getExchangeRxnsIdx = intersect(find(sum(orderedFbaMatrix, 2) ~=0 & sum(orderedFbaMatrix < 0, 2) ~= 0), ...
                                          idxToVis);  
        titleWord = "Negative flux reactions";
    elseif options.thresholdFlux == "lower"
        getExchangeRxnsIdx = intersect(find(sum(orderedFbaMatrix,2) ~=0 & sum(orderedFbaMatrix > 0,2) ~= 0), ...
                                                      idxToVis); 
        titleWord = "Positive flux reactions";
    elseif options.thresholdFlux == "none"
        getExchangeRxnsIdx = intersect(find(sum(orderedFbaMatrix, 2) ~=0), ...
                                                      idxToVis);
        titleWord = options.titlePlots;
    elseif options.thresholdFlux == "all"
        getExchangeRxnsIdx = idxToVis;
        titleWord = options.titlePlots;
    else
        error("Wrong value for chosen threshold. Possible values: 'lower', 'upper', 'none', 'all'.")
    end

    orderedFbaMatrixEx = orderedFbaMatrix(getExchangeRxnsIdx, :);
    orderedMappingRxnMatrixEx = orderedMappingRxnMatrix(getExchangeRxnsIdx, :);
    rxnNames = project.models.(referenceModel).model.rxns(getExchangeRxnsIdx);

    if options.FVA
        replacementValue = "analysis.active.FVA.minMaxFluxes.maxFlux"; % get the fba solution values
        orderedFvaMaxMatrix = getOrderedFeatureMatrix(project, modelList, referenceModel, "rxns", replacementValue);
        replacementValue = "analysis.active.FVA.minMaxFluxes.minFlux"; % get the fba solution values
        orderedFvaMinMatrix = getOrderedFeatureMatrix(project, modelList, referenceModel, "rxns", replacementValue);
        orderedFvaMaxMatrixEx = orderedFvaMaxMatrix(getExchangeRxnsIdx, :);
        orderedFvaMinMatrixEx = orderedFvaMinMatrix(getExchangeRxnsIdx, :);

        if options.thresholdFlux == "all"
            idxFVAVar = ~(all(orderedFvaMaxMatrixEx == 0, 2) & all(orderedFvaMaxMatrixEx == 0, 2));
            orderedFvaMaxMatrixEx = orderedFvaMaxMatrixEx(idxFVAVar, :);
            orderedFvaMinMatrixEx = orderedFvaMinMatrixEx(idxFVAVar, :);
            rxnNames = rxnNames(idxFVAVar);
            orderedFbaMatrixEx = orderedFbaMatrixEx(idxFVAVar, :);
            orderedMappingRxnMatrixEx = orderedMappingRxnMatrixEx(idxFVAVar, :);
        end

        % if options.reducedCost
        % 
        %     models = rmfield(project.models,setdiff(fieldnames(project.models), modelList));
        % 
        %     if all(structfun(@(x) isfield(x.analysis.active.FBA, "w"), models)) && all(structfun(@(x) ismember('one', x.analysis.active.parameters.Value(x.analysis.active.parameters.Analysis == "FBA")), models))
        % 
        %         for m = modelList'
        %             nRxns = length(project.models.(m).model.rxns);
        %             w = project.models.(m).analysis.active.FBA.w;  % 7350×1
        %             project.models.(m).analysis.active.FBA.reducedCost = w(2*nRxns+1 : 3*nRxns);    % net flux variables ← use this one
        %         end
        % 
        %         replacementValue = "analysis.active.FBA.reducedCost"; % get the fba solution values
        %         orderedReducedCostMatrix = getOrderedFeatureMatrix(project, modelList, referenceModel, "rxns", replacementValue);
        %         orderedReducedCostMatrixEx = orderedReducedCostMatrix(getExchangeRxnsIdx, :);
        %         if options.thresholdFlux == "all"
        %             orderedReducedCostMatrixEx = orderedReducedCostMatrixEx(idxFVAVar, :);
        %         end
        %     else
        %         error("The reduced costs could not be found in the .w slot of the FBA analysis. In case you did not use the 'one' minNorm parameter for the FBA, shadowprices might be stored elsewere.")
        %     end
        % end

    end

    rxnFormulas = printRxnFormula(project.models.(referenceModel).model, 'rxnAbbrList', rxnNames, 'printFlag', false);

    % medium composition 
    % do the models have the same  medium ? 
    % if so I can add it as column to the plot, otherwise it is not added
    models = rmfield(project.models, setdiff(fieldnames(project.models), modelList));
    mediaForModels = structfun(@(x) x.settings.medium, models);
    ref = fieldnames(models);
    ref = models.(ref{1,1}).settings.medium;
    mediumIsEqualBetweenModels = all(arrayfun(@(x) isequaln(x, ref), mediaForModels));

    if mediumIsEqualBetweenModels
        if isfield(ref,'mediumComposition') %in case there is no medium defined
            rxnNames = string(rxnNames);
            % checking which column of the medium composition table entails
            % the rxn names -> just checking for the column with most
            % overlap
            isText = varfun(@(x) iscell(x) || isstring(x) || ischar(x), ...
                            ref.mediumComposition, ...
                            "OutputFormat", "uniform");

            textTable = ref.mediumComposition(:, isText);
            nMatches = varfun(@(x) sum(ismember(string(x), rxnNames)), textTable);
            [maxMatches, bestColumn] = max(nMatches{:,:});
            bestColumnName = ref.mediumComposition.Properties.VariableNames{bestColumn};

            mediumConstrained = ismember(rxnNames, ref.mediumComposition.(bestColumn)); %| ...
                             % ismember(rxns_names, ref.manual_set_boundaries.unwanted_export) | ...
                             % ismember(rxns_names, ref.manual_set_boundaries.unwanted_import);
        end

        orderedLb = getOrderedFeatureMatrix(project, referenceModel, ...
                                             referenceModel, "rxns", "model.lb");
        orderedUb = getOrderedFeatureMatrix(project, referenceModel, ...
                                             referenceModel, "rxns", "model.ub");
        orderedUb = orderedUb(getExchangeRxnsIdx, :);
        orderedLb = orderedLb(getExchangeRxnsIdx, :);

        if options.FVA && options.thresholdFlux == "all"
            orderedLb = orderedLb(idxFVAVar, :);
            orderedUb = orderedUb(idxFVAVar, :);
        end

        % get rxn gene rules to add to the table
        symbolGprRules = string(cellfun(@(rxnName)getRxnSymbolRule(project.models.(referenceModel), ...
                                                       rxnName), string(rxnNames), 'UniformOutput', false));

        T = table(rxnFormulas, mediumConstrained,...
            join(string(orderedMappingRxnMatrixEx), "|", 2), symbolGprRules, ...
                  'VariableNames', ["Reaction Formula", "Medium Constrained", ...
                                    join(string(project.comparisons.(comparisonName).modelNames), "_"), ...
                                    "Symbol GPR Rules"], 'RowNames', rxnNames);
        T = T(flip(string(T.Properties.RowNames)), :);
        T.lb = flip(orderedLb);
        T.ub = flip(orderedUb);
    end

    % if options.shadowPrice
    %     replacementValue = "analysis.active.FBA.basis.dual"; % get the fba solution values
    %     % shadow prices are measured for every metabolite therefore mapped according to the mets field
    %     ordered_shadowPrices_matrix = getOrderedFeatureMatrix(project, modelList, referenceModel, "mets", replacementValue);
    % end

    %%% Parameters for the figure 

    %%% =========================
    % Threshold & masks
    % =========================
    mu = mean(orderedFbaMatrixEx, 2);

    x  = abs(mu);

    % 1D k-means to minimize within-group distances

    assert(length(x) >1,...
           'Only one reaction value was returned after filtering for the active reactions in the FBA solution. In order to visualize this plot, you need to formulate more rxns in rxnID or to apply other threholds (try "all" threshold).')

    [idx, C] = kmeans(x, 2, 'Replicates', 10);

    % Identify low-flux (dense) cluster
    [~, lowCluster] = min(C);

    isLow = idx == lowCluster;
    isHigh = ~isLow;

    % Derive lowRange for documentation
    % lowRange = [-2 2];
    lowRange = [-max(x(isLow)), max(x(isLow))];

    isLow  = mu > lowRange(1) & mu < lowRange(2);
    isHigh = ~isLow;

    dataLow  = orderedFbaMatrixEx(isLow, :);
    dataHigh = orderedFbaMatrixEx(isHigh, :);

    rxnNamesLow  = rxnNames(isLow);
    rxnNamesHigh = rxnNames(isHigh);

    %%% =========================
    % Adaptive relative heights
    % =========================
    alpha = 0.75;
    nLow  = numel(rxnNamesLow);
    nHigh = numel(rxnNamesHigh);

    w = [nLow nHigh].^alpha;
    usableHeight = 1 - 0.08 - 0.10 - 0.03;

    heightLow  = usableHeight * w(1) / sum(w);
    heightHigh = usableHeight * w(2) / sum(w);

    bottomLow  = 0.10;
    bottomHigh = bottomLow + heightLow + 0.03;
    %%%

    fig = uifigure('Name', titleWord + " with Table", ...
                       'Position', [100 100 1000 450], 'Visible', options.visiblePlots);

    plotWidth = 0.52;

    %%% =========================
    % Helper for axes formatting
    % =========================
    fmtAxes = @(ax) set(ax, ...
        'FontSize', 16, ...
        'PositionConstraint', 'outerposition', ...
        'GridColor', [.8 .8 .8], ...
        'GridAlpha', .5);

    if ~options.FVA %&& ~options.reducedCost % specify which of the fields need to be true and false!!!
        % in the case that only the FBA solution should be visualized we
        % use a grouped horizontal barplot to do so

        %%% =========================
        % TOP AXIS — high values
        % =========================
        if heightHigh ~= 0
            axTop = uiaxes(fig, 'Units', 'normalized', ...
                'Position', [0.02 bottomHigh plotWidth heightHigh]);
            fmtAxes(axTop)

            if size(dataHigh, 1) == 1
                dataHigh = [dataHigh; NaN(1, size(dataHigh, 2))];
                dataHighOneBar = 1;
            else
                dataHighOneBar = 0;
            end

            barh(axTop, dataHigh, 'grouped');
            grid(axTop, 'on')

            %axTop.YTick = [];
            %axTop.YColor = 'none';

            yticks(axTop, 1:numel(rxnNamesHigh))
            yticklabels(axTop, strrep(rxnNamesHigh, "_", "\_"))
            title(axTop, titleWord)
            if heightLow == 0
                xlabel(axTop, 'Flux value [µMol/(gDW*h)]')
            end
            if dataHighOneBar
                ylim(axTop, [0.5, size(dataHigh, 1) - 0.4 ])
            end
        end

        %%% =========================
        % BOTTOM AXIS — low values
        % =========================
        if heightLow ~= 0
            axBottom = uiaxes(fig, 'Units', 'normalized', ...
                'Position', [0.02 bottomLow plotWidth heightLow]);
            fmtAxes(axBottom)

            if size(dataLow, 1) == 1
                dataLow = [dataLow; NaN(1, size(dataLow, 2))];
                dataLowOneBar = 1; 
            else
                dataLowOneBar = 0;
            end
            barh(axBottom, dataLow, 'grouped');
            grid(axBottom, 'on')

            yticks(axBottom, 1:numel(rxnNamesLow))
            yticklabels(axBottom, strrep(rxnNamesLow, "_", "\_"))
            if heightHigh == 0
                title(axBottom, titleWord)
            end
            xlabel(axBottom, 'Flux value [µMol/(gDW*h)]')
            if dataLowOneBar
                ylim(axBottom, [0.5, size(dataLow, 1) - 0.4 ])
            end
        end

        %%% =========================
        % Reverse direction if needed
        % =========================
        if options.thresholdFlux == "upper"
            set([axTop axBottom], 'XDir', 'reverse')
        end

        %%% =========================
        % Legend
        if nHigh < nLow
            legend(axBottom, modelList, 'Location', 'best');
        else
            legend(axTop, modelList, 'Location', 'best');
        end

    else

        Q1_low  = orderedFvaMinMatrixEx(isLow, :);
        MED_low = orderedFbaMatrixEx(isLow, :);
        Q3_low = orderedFvaMaxMatrixEx(isLow, :);

        Q1_high  = orderedFvaMinMatrixEx(isHigh, :);
        MED_high = orderedFbaMatrixEx(isHigh, :);
        Q3_high = orderedFvaMaxMatrixEx(isHigh, :);

        % =========================
        % Check if reduced cost coloring is enabled
        % =========================
        %addColorLegend = options.reducedCost;

        % if addColorLegend
        %     cmap = colormap('cool');  % N colors
        %     N = size(cmap, 1);
        % 
        %     % Split reduced cost into high and low
        %     reducedCostLow  = orderedReducedCostMatrixEx(isLow, :);
        %     reducedCostHigh = orderedReducedCostMatrixEx(isHigh, :);
        % 
        %     % Helper function to scale matrix to colormap indices
        %     scaleToCmap = @(mat, valMin, valMax) round((mat - valMin) / (valMax - valMin) * (N-1)) + 1;
        % 
        %     % Get global min/max from full matrix
        %     valMin = min(orderedReducedCostMatrixEx(:));
        %     valMax = max(orderedReducedCostMatrixEx(:));
        % 
        %     if valMin == valMax
        %         % Entire matrix has one unique value — all same color
        %         scaledIdxLow = ones(size(reducedCostLow));
        %         scaledIdxHigh = ones(size(reducedCostHigh));
        %     else
        %         % Check each subset individually
        %         if length(unique(reducedCostLow(:))) == 1
        %             scaledIdxLow = ones(size(reducedCostLow));
        %         else
        %             scaledIdxLow = scaleToCmap(reducedCostLow, valMin, valMax);
        %         end
        % 
        %         if length(unique(reducedCostHigh(:))) == 1
        %             scaledIdxHigh = ones(size(reducedCostHigh));
        %         else
        %             scaledIdxHigh = scaleToCmap(reducedCostHigh, valMin, valMax);
        %         end
        %     end
        % else
        %     addColorLegend = 0;
        % end

       %%

        % --- TOP AXES: high flux ---
        %%% =========================
        % Plot high reactions
        %%% =========================
        if nHigh > 0
            axTop = uiaxes(fig, 'Units', 'normalized', 'Position', [0.02 bottomHigh plotWidth heightHigh]);
            hold(axTop, 'on'); 
            fmtAxes(axTop);

            % --- Grid & font ---
            grid(axTop, 'on')
            axTop.GridColor = [0.8 0.8 0.8];
            axTop.GridAlpha = 0.5;
            axTop.FontSize = 16;

            nGroups = nHigh;
            nPerGroup = size(Q1_high, 2);
            boxHeight = 0.2; 
            groupSep = 1;
            greyLevels = linspace(0.9, 0.6, nPerGroup);
            greyRGBs = [greyLevels', greyLevels', greyLevels'];

            for i = 1:nGroups
                yBase = i*groupSep;
                for j = 1:nPerGroup
                    offset = (j-(nPerGroup+1)/2)*(boxHeight*1.3);
                    rectangle(axTop, 'Position', [Q1_high(i,j), yBase+offset-boxHeight/2, Q3_high(i,j)-Q1_high(i,j), boxHeight],...
                              'FaceColor', greyRGBs(j, :), 'EdgeColor', 'none');
                    % Median dot
                    % if addColorLegend
                    %     plot(axTop, MED_high(i, j), yBase+offset, 'o', ...
                    %     'MarkerSize', 5, ...
                    %     'MarkerFaceColor', cmap(scaledIdxHigh(i, j), :), ...
                    %     'MarkerEdgeColor', cmap(scaledIdxHigh(i, j), :));
                    % else
                        plot(axTop, MED_high(i, j), yBase+offset, 'o', 'MarkerSize', 5, 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k')
                    %end
                end
            end

            yticks(axTop, 1:nGroups*groupSep)
            yticklabels(axTop, strrep(rxnNamesHigh, "_", "\_"))
            title(axTop, titleWord)
        end

        %%% =========================
        % Plot low reactions
        %%% =========================
        if nLow > 0

            axBottom = uiaxes(fig, 'Units', 'normalized', 'Position', [0.02 bottomLow plotWidth heightLow]);
            hold(axBottom, 'on'); 
            fmtAxes(axBottom);

            % --- Grid & font ---
            grid(axBottom, 'on')
            axBottom.GridColor = [0.8 0.8 0.8];
            axBottom.GridAlpha = 0.5;
            axBottom.FontSize = 16;

            nGroups = nLow;
            nPerGroup = size(Q1_low, 2);
            boxHeight = 0.2; 
            groupSep = 1;
            greyLevels = linspace(0.9, 0.6, nPerGroup);
            greyRGBs = [greyLevels', greyLevels', greyLevels'];

            for i = 1:nGroups
                yBase = i*groupSep;
                for j = 1:nPerGroup
                    offset = (j-(nPerGroup+1)/2)*(boxHeight*1.3);
                    rectangle(axBottom, 'Position', [Q1_low(i,j), yBase+offset-boxHeight/2, Q3_low(i,j)-Q1_low(i,j), boxHeight], ...
                              'FaceColor', greyRGBs(j,:), 'EdgeColor', 'none');
                    % Median dot
                    % if addColorLegend
                    %     plot(axBottom, MED_low(i,j), yBase+offset, 'o', ...
                    %     'MarkerSize', 5, ...
                    %     'MarkerFaceColor', cmap(scaledIdxLow(i, j), :), ...
                    %     'MarkerEdgeColor', cmap(scaledIdxLow(i, j), :));
                    % else
                        plot(axBottom, MED_low(i, j), yBase+offset, 'o', 'MarkerSize', 5, 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k')
                %     end
                end
            end

            yticks(axBottom, 1:nGroups*groupSep)
            yticklabels(axBottom, strrep(rxnNamesLow, "_", "\_"))
            xlabel(axBottom, 'Flux value [µMol/(gDW*h)]')
        end

        %%% =========================
        % Set axes limits based on MED only
        %%% =========================
        if exist('axTop', 'var')
            % X limits
            xMinHigh = min(MED_high(:));
            xMaxHigh = max(MED_high(:));
            % Add small padding
            xPad = (xMaxHigh - xMinHigh) * 0.05;
            % xlim(axTop, [xMinHigh-xPad, xMaxHigh+xPad]);

            % Y limits
            yMinHigh = 0.5;  % first y
            yMaxHigh = nHigh * groupSep + 0.5;
            ylim(axTop, [yMinHigh, yMaxHigh]);
        end

        if exist('axBottom','var')
            % X limits
            xMinLow = min(MED_low(:));
            xMaxLow = max(MED_low(:));
            xPad = (xMaxLow - xMinLow)*0.05;
            if xMinLow-xPad == xMaxLow-xPad
                xlim(axBottom, [-1, 1]);
            else
                xlim(axBottom, [xMinLow-xPad, xMaxLow+xPad]);
            end

            % Y limits
            yMinLow = 0.5;
            yMaxLow = nLow * groupSep + 0.5;
            ylim(axBottom, [yMinLow, yMaxLow]);
        end

        % --- Add colorbar if needed
        % if addColorLegend
        %     % Choose axTop if it exists, else axBottom
        %     if exist('axTop', 'var')
        %         cbAx = axTop;
        %     else
        %         cbAx = axBottom;
        %     end
        % 
        %     colormap(cbAx, cmap);          % Set colormap for chosen axis
        %     % --- Set caxis safely
        %     if valMin == valMax
        %         caxis(cbAx, [valMin-0.5, valMax+0.5]);  % small padding if single value
        %     else
        %         caxis(cbAx, [valMin, valMax]);
        %     end
        % 
        %     cb = colorbar(cbAx);           % Attach colorbar
        %     cb.Label.String = ...
        %        'Reduced Cost';
        %     cb.FontSize = 14;
        % end

        % --- Grey patches legend (for models)
        % Choose axTop if it exists, else axBottom
        if exist('axTop', 'var')
            lgdAx = axTop;
        else
            lgdAx = axBottom;
        end

        hGrey = gobjects(nPerGroup, 1);
        for j = 1:nPerGroup
            hGrey(j) = patch(lgdAx, NaN, NaN, greyRGBs(j,:), 'EdgeColor', 'none');
        end

        % Add legend
        lgd = legend(lgdAx, hGrey, modelList, 'Location', 'northeastoutside');
        % lgd = legend(lgdAx, hGrey, modelList, 'Location','northeast');
        lgd.FontSize = 16;
        lgd.Box = 'off';
        lgd.Color = 'none';
        lgd.Title.String = "Models";

    end

    % --- Reverse direction if needed ---
    if options.thresholdFlux == "upper"
        set([axTop axBottom], 'XDir', 'reverse')
    end

    %%% =========================
    % Table
    % =========================
    % resort the table according to the order in the plots
    if exist('T', 'var')
        T = T([flip(string(rxnNamesHigh)); flip(string(rxnNamesLow))], :);
        tbl = uitable(fig, ...
                        'Data', T, ...
                        'ColumnName', T.Properties.VariableNames, ...
                        'Units', 'normalized', ...
                        'Position', [plotWidth+0.05 0.10 0.40 0.85], ... % width increased from 0.30 → 0.40
                        'FontSize', 16, ...
                        'ColumnWidth', 'auto');
         tbl.FontSize = 16; 
    end

end

visualizeSamplingLandscape(project, comparison_name, rxn_to_visualize="biomass_reaction", dim_reduction_type="UMAP", pcs_vis=[1,2], sampling_feature="flux", num_clusters=0, pcs_used_dim_red=0, perform_kmeans=0, thinning=10, n_neighbors=50, overwrite=0, visible_plot="on")

Visualizes sampling solutions in a dimension-reduced space (PCA or UMAP) with reaction flux or flux sum values overlaid as color. PCA is always computed first; UMAP is then applied on the selected principal components. Optionally, k-means clustering can be performed on the PCA-reduced data.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison.

 required
comparison_name string

Name of the comparison containing the sampling results to visualize.

 required
rxn_to_visualize string

Reaction name whose flux values are displayed as color in the reduced space. Default: "biomass_reaction".

 required
options.dim_reduction_type string

Dimension reduction method, "PCA" or "UMAP". Default: "UMAP".

 required
options.pcs_vis 1-by-2 array

Principal components to display when using PCA. Default: [1, 2].

 required
options.sampling_feature string

Feature space for reduction, "flux" (per-reaction) or "fluxsum" (per-metabolite). Default: "flux".

 required
options.num_clusters numeric

Number of k-means clusters. If 0, defaults to the number of unique model labels. Default: 0.

 required
options.pcs_used_dim_red numeric

Number of PCs fed into UMAP. If 0, automatically determined to reach 70 percent cumulative variance. Default: 0.

 required
options.perform_kmeans numeric

If 1, run k-means clustering on the PCA-reduced data. Default: 0.

 required
options.thinning numeric

Subsampling interval to reduce computational load. Every nth sample is kept. Default: 10.

 required
options.n_neighbors numeric

UMAP n_neighbors parameter. Default: 50.

 required
options.overwrite numeric

If 1, overwrite existing dimension reduction results. Default: 0.

 required
options.visible_plot string

Figure visibility, "on" or "off". Default: "on".

 required

Output arguments:

Name Type Description
fig_out struct

Struct containing the generated figures. Fields include "label" (model labels scatter), "cluster" (k-means cluster scatter, if performed), and a field named after the visualized reaction (flux value colored scatter).

Examples:

% Default UMAP visualization with biomass reaction flux as color
fig_out = visualizeSamplingLandscape(project, compName);

% Use PCA with custom components and fluxsum features
fig_out = visualizeSamplingLandscape(project, compName, ...
    "EX_glc(e)", struct('dim_reduction_type', "PCA", ...
    'pcs_vis', [1, 3], 'sampling_feature', "fluxsum"));

% Run UMAP with k-means clustering
fig_out = visualizeSamplingLandscape(project, compName, ...
    "biomass_reaction", struct('perform_kmeans', 1, ...
    'num_clusters', 3, 'thinning', 5));
Note

PCA is always computed on z-scored, zero-variance-filtered samples. UMAP is applied on the first numPCs principal components of the thinned data. K-means clustering quality is assessed via silhouette score and label homogeneity.

Warning

UMAP requires the UMAP toolbox, which is compatible with MATLAB R2019a through R2024b only. R2025a and later remove Java access to MATLAB figures, breaking the toolbox.

Source code in scr/classes/model_comparison/functions/visualizeSamplingLandscape.m 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
function fig_out = visualizeSamplingLandscape(project, comparison_name, rxn_to_visualize, options)
% Visualizes sampling solutions in a dimension-reduced space (PCA or
% UMAP) with reaction flux or flux sum values overlaid as color. PCA is
% always computed first; UMAP is then applied on the selected principal
% components. Optionally, k-means clustering can be performed on the
% PCA-reduced data.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison.
%   comparison_name (string): Name of the comparison containing the
%       sampling results to visualize.
%   rxn_to_visualize (string): Reaction name whose flux values are
%       displayed as color in the reduced space. Default:
%       "biomass_reaction".
%   options.dim_reduction_type (string): Dimension reduction method,
%       "PCA" or "UMAP". Default: "UMAP".
%   options.pcs_vis (1-by-2 array): Principal components to display
%       when using PCA. Default: [1, 2].
%   options.sampling_feature (string): Feature space for reduction,
%       "flux" (per-reaction) or "fluxsum" (per-metabolite). Default:
%       "flux".
%   options.num_clusters (numeric): Number of k-means clusters. If 0,
%       defaults to the number of unique model labels. Default: 0.
%   options.pcs_used_dim_red (numeric): Number of PCs fed into UMAP.
%       If 0, automatically determined to reach 70 percent cumulative
%       variance. Default: 0.
%   options.perform_kmeans (numeric): If 1, run k-means clustering on
%       the PCA-reduced data. Default: 0.
%   options.thinning (numeric): Subsampling interval to reduce
%       computational load. Every nth sample is kept. Default: 10.
%   options.n_neighbors (numeric): UMAP n_neighbors parameter.
%       Default: 50.
%   options.overwrite (numeric): If 1, overwrite existing dimension
%       reduction results. Default: 0.
%   options.visible_plot (string): Figure visibility, "on" or "off".
%       Default: "on".
%
% Returns:
%   fig_out (struct): Struct containing the generated figures. Fields
%       include "label" (model labels scatter), "cluster" (k-means
%       cluster scatter, if performed), and a field named after the
%       visualized reaction (flux value colored scatter).
%
% Examples:
%   ```matlab
%   % Default UMAP visualization with biomass reaction flux as color
%   fig_out = visualizeSamplingLandscape(project, compName);
%
%   % Use PCA with custom components and fluxsum features
%   fig_out = visualizeSamplingLandscape(project, compName, ...
%       "EX_glc(e)", struct('dim_reduction_type', "PCA", ...
%       'pcs_vis', [1, 3], 'sampling_feature', "fluxsum"));
%
%   % Run UMAP with k-means clustering
%   fig_out = visualizeSamplingLandscape(project, compName, ...
%       "biomass_reaction", struct('perform_kmeans', 1, ...
%       'num_clusters', 3, 'thinning', 5));
%   ```
%
% Note:
%   PCA is always computed on z-scored, zero-variance-filtered samples.
%   UMAP is applied on the first numPCs principal components of the
%   thinned data. K-means clustering quality is assessed via silhouette
%   score and label homogeneity.
%
% Warning:
%   UMAP requires the UMAP toolbox, which is compatible with MATLAB
%   R2019a through R2024b only. R2025a and later remove Java access to
%   MATLAB figures, breaking the toolbox.

    arguments
        project 
        comparison_name
        rxn_to_visualize (1,1) string ="biomass_reaction"
        options.dim_reduction_type (1,1) string {mustBeMember(options.dim_reduction_type,["PCA", "UMAP"])} ="UMAP"
        options.pcs_vis (1,2) =[1,2]
        options.sampling_feature (1,1) string {mustBeMember(options.sampling_feature,["flux", "fluxsum"])} ="flux"
        options.num_clusters =0
        options.pcs_used_dim_red =0
        options.perform_kmeans =0
        options.thinning =10
        options.n_neighbors =50
        options.overwrite =0
        options.visible_plot ="on"
    end
    fig_out = struct();
    if options.num_clusters == 0
        options.num_clusters = length(unique(project.comparisons.(comparison_name).samplingComparison.sampleModelLabels));
    end
    reference_model = project.comparisons.(comparison_name).referenceModel;

    sampleModelLabels = project.comparisons.(comparison_name).samplingComparison.sampleModelLabels;
    if options.sampling_feature == "flux"
        ordered_samples = project.comparisons.(comparison_name).samplingComparison.orderedSamples;
        dim_names = project.models.(reference_model).model.rxns;
    elseif options.sampling_feature == "fluxsum"
        ordered_samples = project.comparisons.(comparison_name).samplingComparison.ordered_samples_fluxsum;
        dim_names = project.models.(reference_model).model.mets;
    else
        error("You choose a  non valid argument for the sampling feature! Choose flux or fluxsum, by default the flux per reaction will be portrayed!")
    end

    % create new slot for output of function 
    project.comparisons.(comparison_name).dimension_reduction = struct();

    % run the Principle component analysis
    pc_x = options.pcs_vis(1);
    pc_y = options.pcs_vis(2);

    % filter out dimensions which are always 0
    ordered_samples_filter = ordered_samples(find(any(ordered_samples ~= 0,2)),:);
    ordered_dim_names = dim_names(find(any(ordered_samples ~= 0,2)),:);
    % z-scale the matrix so that all features have the same influence on
    % the dimension reduction
    samples_matrix = ordered_samples_filter';
    scaled_samples = zscore(samples_matrix, 0, 1); 
    % peform pca
    [pca_samp.coeff,pca_samp.score,pca_samp.latent,pca_samp.tsquared,pca_samp.explained] = pca(scaled_samples);
    % determine number of PCs used in downstream analysis or use the number
    % of PCs specified in the options parameter as function input
    if options.pcs_used_dim_red == 0
        cumulativeVariance = cumsum(pca_samp.explained);
        exp_var = 70;
        numPCs = find(cumulativeVariance >= exp_var, 1, 'first');
        fprintf('Number of PCs to reach 70%% variance: %d\n', numPCs);
    else
        numPCs = options.pcs_used_dim_red;
        exp_var = sum(pca_samp.explained(1:numPCs));
    end
    project.comparisons.(comparison_name).dimension_reduction.pca = pca_samp;

    thin = options.thinning;
    % perform umap 
    if options.dim_reduction_type == "UMAP"
        % prepare data for dimension reduction 
        X = pca_samp.score(:,1:numPCs);

        X = X(1:thin:end,:);
        fprintf('Computing umap dimension reduction ...!');
        [reduction,~] = run_umap(X, ...
                                'n_neighbors',options.n_neighbors, ...
                                'min_dist',0.3, ...
                                'verbose', 'none',...  
                                'method', 'C vectorized'); % 'MATLAB vectorized' also works, but much slower! C++ did not work for me
                                % 'MEX' is fastest but needs manual download of files for every execution
        fprintf('...finished!');
        %figure
        %scatter(reduction(:,1),reduction(:,2),20,categorical(labels),'filled')

        project.comparisons.(comparison_name).dimension_reduction.umap.reduction = reduction;
        project.comparisons.(comparison_name).dimension_reduction.umap.n_neighbors = options.n_neighbors;
        project.comparisons.(comparison_name).dimension_reduction.umap.thinning = thin;
        project.comparisons.(comparison_name).dimension_reduction.umap.used_samples_idx = 1:thin:length(sampleModelLabels);
        project.comparisons.(comparison_name).dimension_reduction.umap.used_pcs = numPCs;
    end
    if options.perform_kmeans ==1
        % prepare data for dimension reduction 
        X = pca_samp.score(:,1:numPCs);
        thin = options.thinning;
        X = X(1:thin:end,:);
        labels = sampleModelLabels(1:thin:end);
        kmeans_results = struct();
        [kmeans_results.idx, kmeans_results.C,...
         kmeans_results.sumd,kmeans_results.D] = kmeans(X, options.num_clusters, ...
                                                        'Distance','sqeuclidean', ...
                                                        'Replicates',10, ...
                                                        'MaxIter',500, ...
                                                        'Display','final');

        % figure;
        % scatter(reduction(:,1), reduction(:,2), 20, categorical(idx), 'filled');
        % xlabel('UMAP 1'); ylabel('UMAP 2');
        % title('K-means clusters on PCA space visualized via UMAP');
        % colorbar;

        % check quality of clustering - silhoutte score and homogeneity of
        % the clustering
        [kmeans_results.silhouette,~] = silhouette(X,kmeans_results.idx,'Euclidean', 'MaxIter', 100);
        mean_sil = mean(kmeans_results.silhouette);
        %
        homogen = [];
        for k = unique(kmeans_results.idx)'
            labels_in_cluster = labels(find(kmeans_results.idx == k));
            [s,~,j]=unique(labels_in_cluster);
            f = s{mode(j)};
            m = sum(labels_in_cluster == f);
            homogen = [ homogen, m/length(labels_in_cluster)];
        end
        kmeans_results.homogen = mean(homogen);

        project.comparisons.(comparison_name).dimension_reduction.kmeans = kmeans_results; 
    end

    % Visualization of computed metrics

    if options.dim_reduction_type == "PCA"

         pca_samp = project.comparisons.(comparison_name).dimension_reduction.pca;

         dim1 = pca_samp.score(:,pc_x);
         dim2 = pca_samp.score(:,pc_y);
         labels = sampleModelLabels;


         fig = figure('Position',[20 20 700 300],'Visible',options.visible_plot);
         scatter(dim1,dim2,10,...
                 categorical(labels),'filled')
         add_labels(dim1, dim2, labels)

         title("PC "+ num2str(pc_x) +  " & " +  num2str(pc_y) + " with the model lables!", 'FontSize',18)
         xlabel("PC" + num2str(pc_x) + " var: " + pca_samp.explained(pc_x), 'FontSize',18)
         ylabel("PC" + num2str(pc_y) + " var: " + pca_samp.explained(pc_y), 'FontSize',18)
         hold off
        fig_out.label = fig;

    else

        umap_res = project.comparisons.(comparison_name).dimension_reduction.umap;

        dim1 = umap_res.reduction(:,1);
        dim2 = umap_res.reduction(:,2);
        labels = sampleModelLabels(umap_res.used_samples_idx);
        dim_used = umap_res.used_pcs;


        fig = figure('Position',[20 20 700 300],'Visible',options.visible_plot);

        scatter(dim1,dim2,10,...
                 categorical(labels),'filled')
        add_labels(dim1, dim2, labels)

        title("UMAP with the model lables", 'FontSize',18)
        xlabel("UMAP1 " + num2str(dim_used) + "PCs -> 2D ,var: " + exp_var, 'FontSize',14)
        ylabel("UMAP2 " + num2str(dim_used) + "PCs -> 2D, var: " + exp_var, 'FontSize',14)

        hold off
        fig_out.label = fig;

    end

    if options.perform_kmeans

        if options.dim_reduction_type == "PCA"
             pca_samp = project.comparisons.(comparison_name).dimension_reduction.pca;
             umap_res = project.comparisons.(comparison_name).dimension_reduction.umap;

             dim1 = pca_samp.score(umap_res.used_samples_idx,pc_x);
             dim2 = pca_samp.score(umap_res.used_samples_idx,pc_y);
             labels = "cluster " + string(project.comparisons.(comparison_name).dimension_reduction.kmeans.idx)';

             fig2 = figure('Position',[20 20 700 300],'Visible',options.visible_plot);

             scatter(dim1,dim2,10,...
                     categorical(labels),'filled')
             add_labels(dim1, dim2, labels)

             title("PC "+ num2str(pc_x) +  " & " +  num2str(pc_y) + " with the cluster assignment from the kmeans!", 'FontSize',18)
             xlabel("PC" + num2str(pc_x) + " var: " + pca_samp.explained(pc_x), 'FontSize',18)
             ylabel("PC" + num2str(pc_y) + " var: " + pca_samp.explained(pc_y), 'FontSize',18)
             hold off

        else
             umap_res = project.comparisons.(comparison_name).dimension_reduction.umap;
             labels = project.comparisons.(comparison_name).dimension_reduction.kmeans.idx;

            dim1 = umap_res.reduction(:,1);
            dim2 = umap_res.reduction(:,2);
            labels = "cluster " + string(project.comparisons.(comparison_name).dimension_reduction.kmeans.idx)';
            dim_used = umap_res.used_pcs;


            fig2 = figure('Position',[20 20 700 300],'Visible',options.visible_plot);

            scatter(dim1,dim2,10,...
                     categorical(labels),'filled')
            add_labels(dim1, dim2, labels)

            title("UMAP with the model lables", 'FontSize',18)

            xlabel("UMAP1 " + num2str(dim_used) + "PCs -> 2D,var: " + exp_var, 'FontSize',14)
            ylabel("UMAP2 " + num2str(dim_used) + "PCs -> 2D,var: " + exp_var, 'FontSize',14)

            hold off

        end
        fig_out.cluster = fig2;
    end

    % visualize flux(sum) value on the dim reduction
    [~,rxn_idx] = ismember(rxn_to_visualize,project.models.(reference_model).model.rxns);
     value_color = ordered_samples(rxn_idx,:);



        if options.dim_reduction_type == "PCA"

            pca_samp = project.comparisons.(comparison_name).dimension_reduction.pca;

            dim1 = pca_samp.score(:,pc_x);
            dim2 = pca_samp.score(:,pc_y);
            labels = sampleModelLabels;

             fig3 = figure('Position',[20 20 700 300],'Visible',options.visible_plot);

             scatter(dim1,dim2,10,...
                     value_color,'filled')
             add_labels(dim1, dim2, labels)
             colormap(jet)
             colorbar

             title("PC "+ num2str(pc_x) +  " & " +  num2str(pc_y) + " with the cluster assignment from the kmeans!", 'FontSize',18)
             xlabel("PC" + num2str(pc_x) + " var: " + pca_samp.explained(pc_x), 'FontSize',18)
             ylabel("PC" + num2str(pc_y) + " var: " + pca_samp.explained(pc_y), 'FontSize',18)
             hold off

        else
            value_color = value_color(1:thin:end);
             umap_res = project.comparisons.(comparison_name).dimension_reduction.umap;

            dim1 = umap_res.reduction(:,1);
            dim2 = umap_res.reduction(:,2);
            dim_used = umap_res.used_pcs;

            fig3 = figure('Position',[20 20 700 300],'Visible',options.visible_plot);

            scatter(dim1, dim2, 10, value_color, 'filled')
            add_labels(dim1, dim2, labels)
            colormap(jet)
            colorbar

            title("UMAP with the model lables + " + rxn_to_visualize + " flux value" , 'FontSize',18)
            xlabel("UMAP1 " + num2str(dim_used) + "PCs -> 2D,var: " + exp_var, 'FontSize',16)
            ylabel("UMAP2 " + num2str(dim_used) + "PCs -> 2D,var: " + exp_var, 'FontSize',16)

            hold off

        end
    safe_name = matlab.lang.makeValidName(rxn_to_visualize);
    fig_out.(safe_name) = fig3;
end

prepareDataForIDAREVisualization(project, comparisonName, folderPath, options=[])

Prepares project and comparison data for visualization in the IDARE Cytoscape app. Exports all models as SBML XML files (with GPR rules and compartments stripped) and generates reaction and metabolite data tables containing FBA, FVA, sampling, flux sum, and structural presence information for each model in the comparison.

Input arguments:

Name Type Description Default
project struct

Project object from singleModelAnalysis and modelComparison.

 required
comparisonName string

Name of the comparison to export.

 required
folderPath string

Base directory where the timestamped output folder will be created.

 required
options

Reserved for future options. Currently unused.

 required

Examples:

% Export data for IDARE visualization
prepareDataForIDAREVisualization(project, compName, "C:/output/idare");
Note

Creates a timestamped subfolder inside folderPath with two subdirectories: "models" (containing .mat and .xml files for each model plus the reference model) and "data" (containing reaction_data.xlsx and metabolite_data.xlsx). The XML export uses COBRApy via a Python environment to convert MATLAB models to SBML format.

Warning

Requires a configured Python environment with COBRApy installed. The Python path is hardcoded in the exportToXML helper function and must be updated to match your local conda environment path.

Source code in scr/classes/model_comparison/functions/prepareDataForIDAREVisualization.m 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
function prepareDataForIDAREVisualization(project, comparisonName, folderPath, options)
% Prepares project and comparison data for visualization in the IDARE
% Cytoscape app. Exports all models as SBML XML files (with GPR rules
% and compartments stripped) and generates reaction and metabolite data
% tables containing FBA, FVA, sampling, flux sum, and structural
% presence information for each model in the comparison.
%
% Arguments:
%   project (struct): Project object from singleModelAnalysis and
%       modelComparison.
%   comparisonName (string): Name of the comparison to export.
%   folderPath (string): Base directory where the timestamped output
%       folder will be created.
%   options: Reserved for future options. Currently unused.
%
% Examples:
%   ```matlab
%   % Export data for IDARE visualization
%   prepareDataForIDAREVisualization(project, compName, "C:/output/idare");
%   ```
%
% Note:
%   Creates a timestamped subfolder inside folderPath with two
%   subdirectories: "models" (containing .mat and .xml files for each
%   model plus the reference model) and "data" (containing
%   reaction_data.xlsx and metabolite_data.xlsx). The XML export uses
%   COBRApy via a Python environment to convert MATLAB models to SBML
%   format.
%
% Warning:
%   Requires a configured Python environment with COBRApy installed.
%   The Python path is hardcoded in the exportToXML helper function and
%   must be updated to match your local conda environment path.

    arguments
        project 
        comparisonName (1,1) string
        folderPath (1,1) string 
        options =[]
    end