opencobra / cobratoolbox Goto Github PK

Minor issue with the latest update to convertSBMLToCobra (adding libSBML 5.11.0 support).

The update commented out a few lines of code that opened progress bars (lines 78, 163, 195). This causes readCbModel to fail when the function then tries to close these handles, as the variables do not exist.

I'm not sure if the comments were introduced to prevent the progress bar appearance during testing, or if this reflects the intention to remove the bars from future releases. Either uncommenting the lines in question, or removing the subsequent calls of close(h) would fix the issue.

Dear Dr. Fleming,

After installing Gurobi, and after initializing Cobra, I run testAll. This gives an error message somewhere half way. I found out that the problem is with testModels. I have noticed three things with this:

•      The folder testModels did not exist with previous versions of Cobra.
  

•      testModels is not a function, like the functions in the other testing folders
  

•      the second line of testModels.m states “if 1 %by default, do not run this script with testAll”.
  

The subModel structure generated by this function does not include a rxnNames field. One fix could be to add something like this:

[~, ~, IB] = intersect(subModel.rxns, model.rxns, 'stable');
subModel.rxnNames = model.rxnNames(IB);

The sampleSCatterMatrix function does not plot the correlation coefficients. If lines 143-144 have the comment symbol removed, it does plot the coefficients.

(Note: You may replace [ ] with [X] to check the box)

I hereby confirm that I have:

• Run git pull and retried to run my code
• Checked that a similar issue hasn't already been filed

[At line 71 there is a spelling error mfileds --> should be mfields]

Reproduction steps or steps that explain the enhancement

1. [First Step]
2. [Second Step]
3. [Other Steps...]

Expected behavior

[Describe the expected behavior here]

Current/observed behavior

[Describe the currently observed behavior here]

Issue nature (N/A for enhancement)

• MATLAB crashed. If yes, please include a crash report and details how you launch MATLAB.
• Performance issue. If yes, please include more details (see above)
• Issue triggered through a specific action. If yes, please describe what you were doing in more detail.

System information

• pathdef.m and startup.m [Provide details here, e.g. which directories are included and where they are stored]
• MATLAB version: [Enter MATLAB version here]
• OS and version: [Enter OS name and OS version here]
• Installed solvers [Provide a list of installed solvers here]

Additional information (N/A for enhancement)

• Problem might be related to external program (e.g., solver) and not The COBRA Toolbox
• Problem started happening recently, didn't happen in an older version of The COBRA Toolbox
• Problem can be reliably reproduced, doesn't happen randomly
• Problem happens for all COBRA models, not only some

The latest version of optimizeCbModel does not return shadow prices for all the metabolites in the model. It returns one less shadow price that there are metabolites.

When I run git submodule update, I get the following error message:

No submodule mapping found in .gitmodules for path 'external/lusolMex64bit'

Dear Coders of Cobra Thermo tools,

I would like to check the thermodynamic consistency of one of my knockout strains and the fluxes given by an FBA/FVA analysis. Therefore I'd need to check out the deltaG matrices for this specific strain to get the concentration differences for each reaction to keep my fluxes as predicted by the FBA. Could you help me to create my knockout strain model with the vonBertalanffy toolbox and restricting it to the q-vector proposed by my FBA?
Maybe you already created a documentation for the tools?

Thank you,

Sebastian

I've been looking at detectDeadEnds, and when setting removeExternalMets to 1, the behavior seems to be buggy.
All examples are from human RECON 2 website.
Right now, setting it to true would remove any metabolite that has at least one exchange reaction, such as strch1, which has other reactions where it breaks down to D-Glucose.
I have a fix for it on my repository, but I'm not sure about other metabolites - for example, should bvite be removed from the list?
bvite has an exchange reaction, a transport reaction (into the cell) and a demand reaction in the cell. As such, it counts as two different metabolites (the same metabolite in two different compartments counts as two different metabolites). I'm not sure if it should be removed form the list when removeExternalMets=1 or not.

Please let me know what behavior would make sense for such metabolites as bvite, and I'll submit a pull request with the appropriate fix.

Thanks

The extractSubNetwork function call is described as

subModel = extractSubNetwork(model,rxnNames,metNames)

the rxnNames variable should be a cell array of reaction IDs (i.e. model.rxns), not names (model.rxnNames).

the netNames variable is the same.

This ambiguity pops up in other places, too. I'll submit issues as I come across them.

There seems to be aproblem with optKnock. If I run the example code from the cobra protocol

model = readCbModel('ecoli_textbook');
model = changeRxnBounds(model, {'EX_o2(e)', 'EX_glc(e)'}, [0 -20], 'l');
selectedRxns = {model.rxns{ [1, 3:5, 7:8, 10, 12, 15:16, 18, 40:41, 44, 46, 48:49, 51, 53:55, 57, 59:62, 64:68, 71:77, 79:83, 85:86, 89:95]}};

options.targetRxn = 'EX_lac-D(e)';
options.vMax = 1000;
options.numDel = 5;
options.numDelSense = 'L';
constrOpt.rxnList = {'Biomass_Ecoli_core_N(w/GAM)-Nmet2', 'ATPM'};
constrOpt.values = [0.05, 8.39];
constrOpt.sense = 'GE';
optKnockSol = OptKnock(model, selectedRxns, options, constrOpt);

I get the following error messages

Expected one output from a curly brace or dot indexing expression, but
there were 0 results.
Error in OptKnock (line 142)
targetRxnIDirrev = rev2irrev{targetRxnID}(1);
Error in optKnock_test (line 14)
optKnockSol = OptKnock(model, selectedRxns, options, constrOpt);

(Note: You may replace [ ] with [X] to check the box)

I hereby confirm that I have:

• [X ] Run git pull and retried to run my code
• [ X] Checked that a similar issue hasn't already been filed

I suggest to change the default branch of the cobratoolbox to develop.
That would avoid pushing PRs to master by default, while it doesn't change the behavior of git pull.

Lines 124-126:

% Number of inner problem primal varibles with no integers associated with
% them & not part of the special constraint set

n_nint = n - n_int - sum(sel_ic);

This assumes there is no intersection between "variables with integers" (KO candidates) associated and "specially constrained".

This is handled differently when Inint is constructed: Inint = selMatrix(~sel_int & ~sel_ic);, where "intersections" are accounted for.

If there is any reaction that is in both sets in the user input, this will trigger an horzcat error in the corresponding matrix construction line:

A_bl = [A_bl; Inint sparse(n_nint,m+n+2n_int+n_ic+2n_m)];

Potential Solution: If there shouldn't be any intersections (Specially constrained reactions aren't KO candidates frequently), this should be tested explicitly over user input.

Highly similar to the problem described here https://groups.google.com/forum/#!searchin/cobra-toolbox/gurobi/cobra-toolbox/ZxAVYQorrqQ/GjKFYSgLJ4MJ

When gapFind uses solveCobraMILP it fails because solveCobraMILP code for "gurobi5" solver doesnt create the x and f variable if no solution is found. This causes the code to fail at the end when it checks if x is empty because x has never been created.
At the beginning of the "gurobi5" specific section of solveCobraMILP it creates a resultgurobi structure with empty x and objval, but this is overwritten a bit later when the gurobi function is run. Also even if it wasnt overwritten, the x and objval is only transferred to x and f if gurobi returns success.
If you alter the code to return an empty x and f, then it fails on returning to gapfill...

Hello,
I worked for many month in a model with an earlier version of cobra and saved all the results on excel spreadsheets, and it never had any issue. However, I changed my laptop and after installing the new Cobra version, the function to open excel file has the following error:

model = xls2model('model1.xlsx')
??? Reference to non-existent field 'metCharges'.

Error in ==> xls2model at 281
model.metCharges = columnVector(model.metCharges); % all others have plural for
vector

I have already fixed a couple of duplications in the metabolite names and Ids, but the issue remains the same. Do you guys have any idea of what is a good solution for this???
Thanks a lot!

please comment
line 86
line 94
line 96

In addMetabolite.m the lines 214 to 216 should use metCHEBIID not metChEBIID otherwise the model.metCHEBIID is not updated.

if isfield(model,'metChEBIID')
            model.metChEBIID{end+1,1} = CHEBIID{i};
end

I found the Google Groups thread about the same problem and tried downloading the master branch, but the issue is still present:

>> testAll
Solver set to MPS. This function will output an MPS matrix string for the LP problem
Error in testBuildMPS.m
testBuildMPS did not pass
SWITCH expression must be a scalar or string constant.

Error in changeCobraSolver (line 289)
    switch solverName


Error in testAll (line 79)
             changeCobraSolver(CBT_MIQP_SOLVERx,'MIQP');

Are there any updates on this?

Thanks!

I tried to build a test for the three rFBA functions (in src/rFBA folder) but failed without a proper test model.
I couldn't also find any publicly available rFBA model, even looking at the supplementary information of the corresponding paper.
rFBA was designed by M. Covert and he uploaded a set of functions and their corresponding tests in simtk.org. Those functions perform exactly like the COBRA functions would

• solve static rFBA models,

• solve dynamic rFBA models and

• solve combined rFBA and dFBA models

I suggest to delete the rFBA folder in the cobratoolbox and add in the commit the link to the tested and complete resource by the author of the method:

https://simtk.org/projects/ifba/

(Note: You may replace [ ] with [X] to check the box)

I hereby confirm that I have:

• Run git pull and retried to run my code
• Checked that a similar issue hasn't already been filed

[Please include a short description of problem or enhancement here]

Reproduction steps or steps that explain the enhancement

1. [First Step]
2. [Second Step]
3. [Other Steps...]

Expected behavior

[Describe the expected behavior here]

Current/observed behavior

[Describe the currently observed behavior here]

Issue nature (N/A for enhancement)

• MATLAB crashed. If yes, please include a crash report and details how you launch MATLAB.
• Performance issue. If yes, please include more details (see above)
• Issue triggered through a specific action. If yes, please describe what you were doing in more detail.

System information

• pathdef.m and startup.m [Provide details here, e.g. which directories are included and where they are stored]
• MATLAB version: [Enter MATLAB version here]
• OS and version: [Enter OS name and OS version here]
• Installed solvers [Provide a list of installed solvers here]

Additional information (N/A for enhancement)

• Problem might be related to external program (e.g., solver) and not The COBRA Toolbox
• Problem started happening recently, didn't happen in an older version of The COBRA Toolbox
• Problem can be reliably reproduced, doesn't happen randomly
• Problem happens for all COBRA models, not only some

MOMA does not work with the new gurobi (gurobi6). Below is the error lists

Error using solveCobraQP (line 422)
Unknown solver: gurobi6

Error in MOMA (line 209)
QPsolution = solveCobraQP(QPproblem, 'printLevel', verbFlag-1);

Error in EthanolProduction_GDLS_Mutants_MOMA (line 21)
[solutionDel,solutionWT,totalFluxDiff,solStatus] =
MOMA(modelWT,modelMutant,'max','false')

Hello all,

CobraToolBox's model structure carries redundant information in its structure. Occasionally, this causes problems when it comes to model maintenance and coherence between different representations of the same information.
It would be convenient to remove the following elements:

rev - this information can be generated using the bounds of the reactions, and created on the fly if needed.

rules and rxnGeneMat - both these structures should also be removed as they can be generated from grRules

Best regards,
Alberto

@laurentheirendt

CPLEX version 12.6.2
UBUNTU version 16.04
MATLAB version 2015b

Small example (code):
% 1.1 Load model generic
genericPath = '~/work/sbgCloud/data/models/unpublished/Recon3/';
genericFile = '2016_07_13_Recon3.mat';
load(strcat(genericPath, genericFile));
model = modelRecon3;
clear modelRecon3 modelRecon3model
% test feasability
solFBA_init = optimizeCbModel(model);

Result:
Warning: The function 'cplexlink1262' returned an mxArray with non-temporary scope.

In Cplex/subsref
In solveCobraLP (line 1304)
In optimizeCbModel (line 213)

HI,

cobratoolbox/solvers/cplexoptimset is an empty file which is in conflict with the cplex function cplexoptimset, which generates a cplex specific options structure.
I think this file (as it is empty) was added unintentional and should probably be removed as it can cause errors when using cplex solvers.

The function fails to export gene-reaction association information correctly when genes with unsupported characters are present in the model. To my experience, unsupported characters are ",", "-" and "|".

All reactions, associated with unsupported characters containing genes, will not have any gene-reaction association information in exported FBC v2 file. The gene list in FBC v2 file, however, is unchanged.

This bug can be circumvented by removing/replacing unsupported characters in model.genes and model.grRules before using writeCbModel function.

In parsing the SBML notes field, this function includes the statement

if (isempty(subSystem))
subSystem = 'Exchange';
end

which assigns the SubSystem 'Exchange' to any reaction that has SUBSYSTEM showing up in its notes field but with no subsystem assigned.

The problem with this is that reading an SBML model without assigned subsystems automatically creates this field in the model structure, and then writing that model again to a new SBML file creates those (empty) fields. The new model will then be read with all reactions being assigned to the subsystem 'Exchange.' This could be fixed by simply removing this block of code, unless there is a compelling reason to keep it.

There is a problem in convertSBMLToCobra, if SBML models do not have kineticLaws (which is quit often the case, especially for FBA networks).
You have to test for existence of kinetic laws before trying to access them.

Error in convertSBMLToCobra (line 122)
parameters = modelSBML.reaction(i).kineticLaw.parameter;

Error in demo_example_1 (line 9)
[model] = convertSBMLToCobra(modelSBML, defaultBound);

The patch is:

% Bug Fix for SBML models without kinetic laws
if isempty(modelSBML.reaction(i).kineticLaw)
parameters = [];
elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter')
parameters = modelSBML.reaction(i).kineticLaw.parameter;
else
parameters = [];
end

Greetings Matthias

Hello,

I used xls file format into COBRA Toolbox using the command 'xls2model' while doing simulation for genome scale metabolic models. However, now it is showing the following error message:

Error using xlsread (line 251)
Worksheet 'Reaction List' not found.

Error in xls2model (line 75)
[~, Strings, rxnInfo] = xlsread(fileName,'Reaction List');

Can you help me out?

Your help will be greatly appreciated and thanks in advance,
Narendran

Dear Laurent, Sylvain
this seems to be an error rather than a warning:

solveCobraLPCPLEX found and returned a solution.
Optimal solution found, Objective 0.11222
Elapsed time is 473.701875 seconds.
[�Warning: FBA problem error: Invalid input !!!!]�

Best Maike

optimizeCbModel is not working on the latest download. I get the following error

Output argument "FBAsolution" (and maybe others) not assigned during
call to "optimizeCbModel".
Error in robustnessAnalysis (line 44)
solMin = optimizeCbModel(tmpModel,'min');
Error in AerobicGlucoseBioMassRA (line 20)
[controlFlux, objFlux] = robustnessAnalysis(model,'EX_o2(e)',100)

It works on the previous download

While gurobi6 is supported as an option in changeCobraSolver.m, only gurobi5 is defined in solveCobraLP.m (starting at line 188). Setting the solver as gurobi6 results in problems in (for example) optimizeCbModel.m

I am trying to use identifyCorrelSets to compute cosets and I am getting cosets of size = 1 independent of the threshold used. What do these cosets mean ?

I am having issues running checkCobraModelUnique on my model - 'cell content indices must be greater than 0' at line 51, when renaming metabolites. By the looks of it, this can be traced back to line 49 - the function is looking through the model's reactions for a given metabolite, which doesn't seem right. Then, in line 51, it's trying to rename a reaction, while it should be trying to rename a metabolite.

LPproblem.A is defined after it is used in optimizeCbModel.m file

`% Fill in the RHS vector if not provided
if ~isfield(model,'b')
warning('LP problem has no defined b in S*v=b. b should be defined, for now we assume b=0')
LPproblem.b=zeros(size(LPproblem.A,1),1);
else
LPproblem.b = model.b;
end

% Rest of the LP problem
[m,n] = size(model.S);
LPproblem.A = model.S;
LPproblem.c = model.c;
LPproblem.lb = model.lb;
LPproblem.ub = model.ub;`

Hi there,

I was using the cobra Toolbox with cplex_direct as a solver, and I ran across an annoying warning on a Win 8.1 Machine:

Warning: File 'clone1.log' not found. 
  In solveCobraLP at 867
  In optimizeCbModel at 159
  In reduceModel at 102
  In pFBA at 89 

Since I was doing repeated optimization in a loop, this flooded my command window with useless messages. After a quick tour on Google, I discovered why this code, starting line 867 in cobratoolbox-master\solvers\solveCobraLP.m, bugged:

        if exist('clone1.log','file')
            delete('clone1.log')
        end

It suffers from the inconsistency of the exist() function in MATLAB for relative paths (see here).
A good fix I found (on the mathworks forum here) is the following:

exist('clone1.log', 'file')                % ==> 2
exist(fullfile(cd, 'clone1.log'), 'file')  % ==> 0

Which, applied to our problem, gives:

        if exist(fullfile(cd, 'clone1.log'), 'file')
            delete('clone1.log')
        end

Which ensures you are using an absolute path, and removes the warning.

Hope this helps,

Pierre

I'm just wondering, whether this is intended behaviour or not.
When performing checkMassAndChargeBalance, a missing formula for a metabolite will not lead to reactions containing this metabolite being flagged as unbalanced. Is this intended, if, it should be documented, if not (and I personally would support this) I think the function should be altered to indicate all reactions with metabolites with missing formulas as unbalanced.
My personal stance would be: If a user is aware of placeholder metabolites in their model he should be knowing which reactions contain those, but not flagging them can easily lead to a user missing those. Especially as the function, while flagging the metabolites, does not indicate any issued with those reactions (except if turning the printlevel to 1, but those messages tend to be easily missed in a larger script).
I'm happy to implement a modification, but I would first want to know if the assumption, that reactions with metabolites that miss a formula should be flagged as unbalanced or not.

The unit "millimole per gram dry weight per hour" mmol/gDW/h = 10^-3 mol ⋅ (gDW)^-1 ⋅ (3600 s)^-1 uses an incorrect multiplier. Instead of 3600 seconds per hour, it produces the multiplier 0.000277777777777778:

<unitDefinition id="mmol_per_gDW_per_hr">
  <listOfUnits>
    <unit kind="mole" exponent="1" scale="-3" multiplier="1" offset="0"/>
    <unit kind="gram" exponent="-1" scale="0" multiplier="1" offset="0"/>
    <unit kind="second" exponent="-1" scale="0" multiplier="0.000277777777777778" offset="0"/>
  </listOfUnits>
</unitDefinition>

The reason is that instead of multiplying the numbers of minutes and seconds they are are divided, see

There is a bug in line 544 of readSBML:
metTmp=[metTmp,'[',modelSBML.species(1).compartment,']'];

It should be species(i), not species(1). This bug has the effect that all metabolites end up in the same compartment, which is the compartment of the first metabolite in the model.

When running the first tutorial examples of Orth with the Ecoli core model, I get the following error message when calling checkObjective(model):

Error using fprintf
Function is not defined for sparse inputs.

Error in checkObjective (line 29)
fprintf('%6.4g\t\t%s\t\t%i\t%s\t%i\n',Sij,model.mets{objMetInd(m)},objMetInd(m),model.rxns{objRxnInd(n)},objRxnInd(n))

When I wanted to sample my model I needed to first make it irreversible. However after applying convertToIrreversible function (which worked OK) I couldn't use the obtained model_Irrev for the sampleCbModel.m as it didn't have the grRules field in it. I updated the convertToIrreversible.m function to added that field to the newly created model.

% Build final structure
modelIrrev.S = modelIrrev.S(:,1:cnt);
modelIrrev.ub = columnVector(modelIrrev.ub(1:cnt));
modelIrrev.lb = columnVector(modelIrrev.lb(1:cnt));
modelIrrev.c = columnVector(modelIrrev.c(1:cnt));
modelIrrev.rev = modelIrrev.rev(1:cnt);
modelIrrev.rev = columnVector(modelIrrev.rev == 1);
modelIrrev.rxns = columnVector(modelIrrev.rxns); 
modelIrrev.mets = model.mets;
matchRev = columnVector(matchRev(1:cnt));
modelIrrev.match = matchRev;
if (isfield(model,'b'))
    modelIrrev.b = model.b;
end
if isfield(model,'description')
    modelIrrev.description = [model.description ' irreversible'];
end
if isfield(model,'subSystems')
    modelIrrev.subSystems = model.subSystems(irrev2rev);
end
if isfield(model,'genes')
    modelIrrev.genes = model.genes;
    genemtxtranspose = model.rxnGeneMat';
    modelIrrev.rxnGeneMat = genemtxtranspose(:,irrev2rev)';
    modelIrrev.rules = model.rules(irrev2rev);
    modelIrrev.grRules = model.grRules(irrev2rev); %added to allow model reduction 18/02/2016 AW
end

I encountered the issue with addReaction.m recently in my work with Recon 2. After modifying the model (used removeRxns.m and addReaction.m scripts) I couldn't use the removeRxns.m function anymore as it complained about index exceeding matrix dimensions. Turns out that removeRxns.m were deleting all the information from all the fields in respect to the removed reaction, however addReaction.m did not update certain reaction-related fields in the model (namely: rxnKeggID, rxnConfidenceEcoIDA, rxnConfidenceScores and rxnsboTerm). The same was true for some metabolite-related fields (metCHEBIID, metKeggID, metInchiString, metHepatoNetID, metEHMNID and metHMDB). I have fixed it by including this fields in the if isfield ..end loops in addReaction.m . It is my first post on GitHub so I hope I've put it in a correct place.

% Update model fields
model.rxns{rxnID,1} = rxnName;
if (revFlag)
    model.rev(rxnID,1) = 1;
else
    model.rev(rxnID,1) = 0;
end
model.lb(rxnID,1) = lowerBound;
model.ub(rxnID,1) = upperBound;
model.c(rxnID,1) = objCoeff;

if (isfield(model,'rxnNames'))
    if exist('rxnNameFull','var')
        model.rxnNames{rxnID,1} = rxnNameFull;
    else
        model.rxnNames{rxnID,1} = model.rxns{rxnID};
    end
end
if (isfield(model,'subSystems'))
    model.subSystems{rxnID,1} = subSystem;
end
if isfield(model,'rxnNotes')
    model.rxnNotes{rxnID,1} = '';
end
if isfield(model,'confidenceScores')
    model.confidenceScores{rxnID,1} = '';
end
if isfield(model,'rxnReferences')
    model.rxnReferences{rxnID,1} = '';
end
if isfield(model,'rxnECNumbers')
    model.rxnECNumbers{rxnID,1} = '';
end
% 17/02/2016 AW Fixed four additional fields that were not being updated
if isfield(model,'rxnKeggID')
    model.rxnKeggID{rxnID,1} = '';
end
if isfield(model,'rxnConfidenceEcoIDA')
    model.rxnConfidenceEcoIDA{rxnID,1} = '';
end
if isfield(model,'rxnConfidenceScores')
    model.rxnConfidenceScores{rxnID,1} = '';
end
if isfield(model,'rxnsboTerm')
    model.rxnsboTerm{rxnID,1} = '';
end

% Construct S-matrix column
newMetsCoefs=zeros(0);
for i = 1:length(metaboliteList)
    if (isInModel(i))
        Scolumn(metID(i),1) = stoichCoeffList(i);
    else
        warning(['Metabolite ' metaboliteList{i} ' not in model - added to the model']);
        Scolumn(end+1,1) = stoichCoeffList(i);
        model.mets{end+1,1} = metaboliteList{i};
        newMetsCoefs(end+1) = stoichCoeffList(i);
        if (isfield(model,'metNames'))      %Prompts to add missing info if desired
            model.metNames{end+1,1} = regexprep(metaboliteList{i},'(\[.+\]) | (\(.+\))','') ;
            warning(['Metabolite name for ' metaboliteList{i} ' set to ' model.metNames{end}]);
%             model.metNames(end) = cellstr(input('Enter complete metabolite name, if available:', 's'));
        end
        if (isfield(model,'metFormulas'))
            model.metFormulas{end+1,1} = '';
            warning(['Metabolite formula for ' metaboliteList{i} ' set to ''''']);
%             model.metFormulas(end) = cellstr(input('Enter metabolite chemical formula, if available:', 's'));
        end
        if isfield(model,'metCHEBIID')
            model.metCHEBIID{end+1,1} = ''; %changed from metChEBIID to CHEBIID as this is how this field is named in Recon2 AW
        end
        if isfield(model,'metKeggID')
            model.metKeggID{end+1,1} = ''; %changed AW
        end
        if isfield(model,'metPubChemID')
            model.metPubChemID{end+1,1} = '';
        end
        if isfield(model,'metInchiString')
            model.metInchiString{end+1,1} = ''; %changed AW
        end
        if isfield(model,'metCharge')
            model.metCharge(end+1,1) = 0;
        end
        if isfield(model,'metHepatoNetID')
            model.metHepatoNetID{end+1,1} = ''; %added AW
        end
        if isfield(model,'metEHMNID')
            model.metEHMNID{end+1,1} = ''; %added AW
        end
        if isfield(model,'metHMDB')
            model.metHMDB{end+1,1} = ''; %added AW
        end
    end
end

Hi, I'm trying to implement non-linear objective functions for FBA using cobra toolbox. Looks like optimizeCbModelNLP() is the one I need but it didn't work. I use fmincon() in MATLAB as the nonlinear solver and based on script optimizeCbModelNLP, i created for myself a FBA model with max-Biomass-per-flux-unit (doi:10.1038/msb4100162) as the objective function.
The problem is that by using 100 or even 500 initial random points, the optimization solutions still not converged (values of solution.f were different from each run though they're very close together). There are also warnings shown up 'Matrix is close to singular or badly scaled...'. Any suggestion for that?

The original GIMME article (http://dx.plos.org/10.1371/journal.pcbi.1000082) uses a weighted penalty when re-enabling flux in reactions that do not have genes expressed.
The penalty is weighted by expression.

GIMME algorithm, as implemented in createTissueSpecificModel.m uses binary expression, and then the penalty becomes unweighted.

Why is the difference? Is there another version of GIMME that implements the weighted penalty?

Line 242 of optimizeCbModel is currently

LPproblem2.ub = [model.ub;10000_ones(2_nRxns,1)];

This hard coding of a max upper bound of 10000 for unconstrained reactions can lead to infeasible solutions or other unexpected behavior for simulations that have large (but not unconstrained) upper limits.

When using gurobi 5.60, such behavior can be fixed by relaxing this upper bound by changing line 242 to:

LPproblem2.ub = [model.ub;inf_ones(2_nRxns,1)];

I do not know how other solvers manage "inf".

Hi,

I was attempting to build my own COBRA model from scratch (modeling glycolysis for class), entering in equations using the addReaction() function and ran into problems with reversible/irreversible reaction bounds. Basically, I used the flag to add reversible reactions, but at least in Matlab that doesn't work -- it seems to be a bug. The flag does change the direction of the reaction arrow when printed to the command window and update the model.rev vector, but it does not update the model.lb vector (which is what actually specifies the reversibility of a reaction). I looked at the code for addReaction() and believe I found the problem and fixed it.

Here is the relevant code within the addReaction function:

0104 % Missing arguments
0105 if (nargin < 6 | isempty(lowerBound))
0106     if (oldRxnFlag)
0107         lowerBound = model.lb(rxnID);
0108     else
0109         if (revFlag)
0110             lowerBound = min(model.lb);
0111             if isempty(lowerBound)
0112                 lowerBound=-1000;
0113             end
0114         else
0115             lowerBound = 0;
0116         end
0117     end
0118 end

In line 110, the lower bound is set to the minimum value in the vector of lower bounds already set in the model. In lines 111-113, the lower bound is only set to -1000 if there is already a reaction with lower bound -1000. In other words, the first reaction entered must be a reversible reaction in order for all subsequent reversible reactions to have the correct lower bound. If on the other hand, the first reaction entered is irreversible (as in my case when adding the glycolysis reactions) with lowerBound = 0, the lower bound would be set to the lower bound of that first reaction (lowerBound = 0). Thus, adding reversible reactions (lowerBound = -1000) after adding the first reaction as irreversible (lowerBound = 0) will result in an error for the lower bound of all subsequent reversible reactions.

I believe this bug would be fixed by adding another conditional statement that only uses the minimum of the model.lb vector as the new lower bound if there is a reversible reaction already present (lines 110-112):

0104 % Missing arguments
0105 if (nargin < 6 | isempty(lowerBound))
0106     if (oldRxnFlag)
0107         lowerBound = model.lb(rxnID);
0108     else
0109         if (revFlag)
0110             if sum(model.rev) ~= 0
0111                 lowerBound = min(model.lb);
0112             end
0113             if isempty(lowerBound)
0114                 lowerBound=-1000;
0115             end
0116         else
0117             lowerBound = 0;
0118         end
0119     end
0120 end

I made this change in my version of addReaction(), and it has so far worked without any bugs. This change still ensures that the minimum reaction bound is set to the present minimum (i.e. if the minimum bound was specified to be something other than -1000 this will preserve that), which I believe was the goal of the original code.

Thanks,

James

example -
DM_C00459[c] C00459[c] ->
DM_C00458[c] C00458[c] ->
DM_C00286[c] C00286[c] ->

When attempting to test solvers, MILP test fails. Upon investigation, a syntax error is at fault due to an incorrect merge of git commit. See here:

Removing the offending lines solves the problem and the test proceeds as intended.

Breaks pFBA (line 96).

Hi all,

I am using recon2 model to test COBRA. I try to optimize the model for biomass reaction obtaining first solution. Then I set some constrains (such as no glucose import) and optimize the second model for biomass reaction. When I compare two solutions the Objective Function Value f remain untouched as well as import level of glucose. Why?

these are the commands I used:
run initCobraToolbox.m
load('/home/matteo/Desktop/Gia/generic GEM/Recon2.v04.mat')
reconS=optimizeCbModel(modelR204, 'max')
reconNOGLU=changeRxnBounds(modelR204, 'EX_Glc(e)', 0, 'l')
reconNOGLUS=optimizeCbModel(reconNOGLU, 'max')
printFluxVector(modelR204, [reconNOGLUS.x reconS.x], true, true)

in the resulting list it results : EX_glc(e) -1 -1

test13CFitting is not working.
I have experiences several bugs with solveCobraNLP when trying to use matlab NLP solver

changeCobraSolver('matlab','NLP');

When looking deeply in the code, it seems that several variables are not defined before they are use.

Hope that you will fix the problem

Best,
Adama.

Error code
`Undefined function or variable 'csense'.

Error in solveCobraNLP (line 148)
A1 = [A(csense == 'L',:);-A(csense ==
'G',:)];

Error in fitC13Data (line 76)
NLPsolution = solveCobraNLP(NLPproblem,
'checkNaN', cnan, 'printLevel', printLevel,
'iterationLimit', majorIterationLimit,
'logFile', 'minimize_SNOPT.txt');

Error in testC13Fitting (line 38)
[vout, rout] = fitC13Data(v0,expdata,model,
majorIterationLimit);`

Hi,

I am using Cobra together with parallel toolbox, and the standard in this setting global variables do not get passed from the MATLAB client to the workers. Thus, in every thread the global variables have to be set, e.g. I have to call changeCobraSolver. I think there could be an alternative mechanism for changeCobraSolver which also suits parallel computing.