rmtfleming
commented
3 years ago Hi Havard, please can you make a pull request from a branch on your fork to the opencobra/develop branch? You can use git natively, or instead here are the instructions for contribution: https://opencobra.github.io/cobratoolbox/latest/contributing.html#how-to-contribute using MATLAB.devTools https://opencobra.github.io/MATLAB.devTools/stable/ Regards, Ronan
On Thu, 8 Apr 2021 at 17:59, molversmyr @.***> wrote:
Hi!
I am using the implementation of mCADRE in the latest version of COBRA Toolbox to create context-specific models of an unpublished model. However, in my work doing so I came over something I suspect might be a bug in the function findFluxConsistentSubset. The corresponding file of this function gets updated as a part of commit 0355eeb https://github.com/opencobra/cobratoolbox/commit/0355eeba0b20d5f5c816d7230b8707973f92b38f. The part that I'm confused about is at line 203-205, where an error is thrown if the inputted model does not return a feasible solution.
From the tutorial page for extraction of context-specific models on the COBRA Toolbox webpage https://opencobra.github.io/cobratoolbox/latest/tutorials/tutorialExtractionTranscriptomic.html, the mCADRE-algorithm is explained as:
A set of core reactions is first found and all other reactions are then pruned based on their expression, connectivity to core, and confidence score. Reactions that are not necessary to support the core or defined functionalities are thus removed.
If I understand this correctly, reactions from a ranked set of non-core reactions are sequentially eliminated from the model if certain conditions are satisfied. In the event that the conditions are not satisfied, e.g. that the temporarily reduced model (not containing the current reaction) does not return an optimal solution, the algorithm should move on to the next reaction in the non-core reaction set. However, the newly added lines of code mentioned above does not allow this to happen. This is also the case for the E. coli-core model used in the tutorial page https://opencobra.github.io/cobratoolbox/latest/tutorials/tutorialExtractionTranscriptomic.html. When following the tutorial, the reduction of this model with mCADRE returns an error ("model is infeasible/unbounded") and thus stops when trying to remove reaction number 10, i.e. EX_glu_L(e),
To overcome this, me and my colleagues made a few changes to the pruningModel.m-file which calls the function in question. I have forked the master-branch of this repo, created a new branch and pushed the changes to the fork on my profile in the commit d62f7b8 https://github.com/opencobra/cobratoolbox/commit/d62f7b837db160aa36f07aa8485b6c20adec9dab. In short, a flag was created to check whether the aforementioned error was thrown or not, removing the reaction from the non-core reaction set before moving on to the next reaction if an error is thrown or continuing to the rest of the code in this function if no error is thrown.
I sincerely hope that this is of some help. Best regards, Håvard Molversmyr
**I hereby confirm that I have:
- Tried to solve the issue on my own
- Retried to run my code with the latest version of The COBRA Toolbox
- Checked that a similar issue has not already been opened
Below is the code to reproduce the problem with the E. coli-core model, as well as the output created before and after implementation of our solution: The code:
%% Initialise The COBRA Toolbox
initCobraToolbox(false) %don't update the toolbox
%% Change COBRA solver to use glpk
changeCobraSolver ('glpk', 'all');
%% Load data and model
modelFileName = 'ecoli_core_model.mat';
modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models.
modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded
model = readCbModel(modelFileName);
load('dataEcoli');
%% Set mCADRE options struct
options = 'options_mCADRE';
load(['options_methods' filesep options]);
options.tol = 1e-06;
%% Create mCADRE-model
mCADRE_model = createTissueSpecificModel(model, options);
The original problem:
Processing inputs and ranking reactions...
Generic model passed precursor metabolites test
Pruning reactions...
Reaction no. 1
Attempting to remove reaction O2t...
Warning: No metabolites defined to check the model function
Warning: No metabolites defined to check the model function
Removed non-core inactive reactions
Num. removed: 11 (0 core, 11 non-core); Num. remaining: 70
Reaction no. 2
Attempting to remove reaction ACALDt...
[…]
Reaction no. 10
Attempting to remove reaction EX_glu_L(e)...
Warning: No metabolites defined to check the model function
sol =
struct with fields:
solver: 'glpk' algorithm: 'default' stat: 0 origStat: 110 origStatText: [] time: 0.0086 basis: [] f: NaN v: [] y: [] w: [] s: [] x: []Error using findFluxConsistentSubset (line 206)
model is infeasible/unbounded
Error in pruningModel (line 69)
[fluxConsistentMetBool,fluxConsistentRxnBool] = findFluxConsistentSubset(model_rem,paramConsistency);Error in mCADRE (line 96)
[tissueModel, cRes] = pruningModel(model, rankNonCore, coreRxn, zeroExpRxns, precursorMets, eta, tol);Error in createTissueSpecificModel (line 247)
tissueModel = mCADRE(model, options.ubiquityScore, options.confidenceScores, options.protectedRxns, options.checkFunctionality, options.eta, options.tol);Error in ecoli_core_test (line 35)
mCADRE_model = createTissueSpecificModel(model, options);
Output after implementing our solution:
Processing inputs and ranking reactions...
Generic model passed precursor metabolites test
Pruning reactions...
Reaction no. 1
Attempting to remove reaction O2t...
Warning: No metabolites defined to check the model function
Warning: No metabolites defined to check the model function
Removed non-core inactive reactions
Num. removed: 11 (0 core, 11 non-core); Num. remaining: 70
Reaction no. 2
Attempting to remove reaction ACALDt...
[…]
Reaction no. 10
Attempting to remove reaction EX_glu_L(e)...
Warning: No metabolites defined to check the model function
sol =
struct with fields:
solver: 'glpk' algorithm: 'default' stat: 0 origStat: 110 origStatText: [] time: 0.0131 basis: [] f: NaN v: [] y: [] w: [] s: [] x: []Num. removed: 29 (0 core, 29 non-core); Num. remaining: 49
Reaction no. 11
Attempting to remove reaction EX_akg(e)...
[…]
Reaction no. 53
Attempting to remove reaction TPI...
Warning: No metabolites defined to check the model function
Warning: No metabolites defined to check the model function
No reactions removed
Num. removed: 49 (0 core, 49 non-core); Num. remaining: 0
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/opencobra/cobratoolbox/issues/1752, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQMEOVMUO7MM7W5OKX7HMLTHXOFXANCNFSM42TJBGMA .
--
Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.
Assistant Professor, Division of Systems Biomedicine and Pharmacology, Leiden Academic Centre for Drug Research, Faculty of Science, Leiden University. https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming & H2020 Project Coordinator, Systems Medicine of Mitochondrial Parkinson’s Disease, http://sysmedpd.eu & Senior Lecturer, School of Medicine, National University of Ireland, Galway.
Peer-reviewed publications: https://goo.gl/FZPG23 Mobile: +353 852 109 806 Skype: ronan.fleming
(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)
molversmyr
Hi!
I am using the implementation of mCADRE in the latest version of COBRA Toolbox to create context-specific models of an unpublished model. However, in my work doing so I came over something I suspect might be a bug in the function
findFluxConsistentSubset. The corresponding file of this function gets updated as a part of commit 0355eeba0b20d5f5c816d7230b8707973f92b38f. The part that I'm confused about is at line 203-205, where an error is thrown if the inputted model does not return a feasible solution.From the tutorial page for extraction of context-specific models on the COBRA Toolbox webpage, the mCADRE-algorithm is explained as:
If I understand this correctly, reactions from a ranked set of non-core reactions are sequentially eliminated from the model if certain conditions are satisfied. In the event that the conditions are not satisfied, e.g. that the temporarily reduced model (not containing the current reaction) does not return an optimal solution, the algorithm should move on to the next reaction in the non-core reaction set. However, the newly added lines of code mentioned above does not allow this to happen. This is also the case for the E. coli-core model used in the tutorial page. When following the tutorial, the reduction of this model with mCADRE returns an error ("model is infeasible/unbounded") and thus stops when trying to remove reaction number 10, i.e. _EX_gluL(e),
To overcome this, me and my colleagues made a few changes to the
pruningModel.m-file which calls the function in question. I have forked the master-branch of this repo, created a new branch and pushed the changes to the fork on my profile in the commit d62f7b837db160aa36f07aa8485b6c20adec9dab. In short, a flag was created to check whether the aforementioned error was thrown or not, removing the reaction from the non-core reaction set before moving on to the next reaction if an error is thrown or continuing to the rest of the code in this function if no error is thrown.I sincerely hope that this is of some help. Best regards, Håvard Molversmyr
**I hereby confirm that I have:
Below is the code to reproduce the problem with the E. coli-core model, as well as the output created before and after implementation of our solution:
The code:
The original problem:
Output after implementing our solution: