Feat/growth check

ComplementaryScripts/increaseVersion.m

@@ -54,7 +54,7 @@ function increaseVersion(bumpType)
 
 %Include tag and save model:
 model.modelID = ['yeastGEM_v' newVersion];
-saveYeastModel(model,false)
+saveYeastModel(model,false,false)   %only save if model can grow
 
 %Check if any file changed (except for history.md and 1 line in yeastGEM.xml):
 diff   = git('diff --numstat');

ComplementaryScripts/otherChanges/anaerobicModel.m

@@ -1,17 +1,20 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% anaerobicModel.m
-% Converts model to anaerobic
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 function model = anaerobicModel(model)
+% anaerobicModel
+%   Converts model to anaerobic.
+%
+%   Inputs: model           (struct) aerobic model
+%   Output: model           (struct) anaerobic model
+%   
+%   Usage: model = anaerobicModel(model)
+%
 
 %1st change: Refit GAM and NGAM to exp. data, change biomass composition
 GAM   = 30.49;  %Data from Nissen et al. 1997
 P     = 0.461;  %Data from Nissen et al. 1997
 NGAM  = 0;      %Refit done in Jouthen et al. 2012
 
 model = changeGAM(model,GAM,NGAM);
-model = scaleBioMass(model,'protein',P,'carbohydrate');
+model = scaleBioMass(model,'protein',P,'carbohydrate',false);
 
 %2nd change: Removes the requirement of heme a, NAD(PH), coenzyme A in the biomass equation
 %            (not used under anaerobic conditions)

ComplementaryScripts/otherChanges/scaleBioMass.m

@@ -1,4 +1,4 @@
-function model = scaleBioMass(model,component,new_value,balance_out)
+function model = scaleBioMass(model,component,new_value,balance_out,dispOutput)
   % scaleBioMass
   %   Scales the biomass composition
   %
@@ -7,14 +7,19 @@
   %   new_value      (float) new total fraction for said component
   %   balance_out    (str, opt) if chosen, the name of another component with which
   %                  the model will be balanced out so that the total mass remains = 1 g/gDW
+  %   dispOutput     (bool, opt) if output from sumBioMass should be displayed (default = true)
   %
   %   model          (struct) modified model
   %
-  %   Usage: model = scaleBioMass(model,component,new_value,balance_out)
+  %   Usage: model = scaleBioMass(model,component,new_value,balance_out,dispOutput)
   %
 
+if nargin < 5
+    dispOutput = true;
+end
+
 %Measure current composition and rescale:
-[~,P,C,R,D,L,I,F] = sumBioMass(model);
+[~,P,C,R,D,L,I,F] = sumBioMass(model,dispOutput);
 content_all = {'carbohydrate','protein','lipid','RNA','DNA','ion','cofactor'};
 content_Cap = {'C','P','L','R','D','I','F'};
 pos         = strcmp(content_all,component);

ComplementaryScripts/otherChanges/sumBioMass.m

@@ -1,23 +1,28 @@
-function [X,P,C,R,D,L,I,F] = sumBioMass(model)
+function [X,P,C,R,D,L,I,F] = sumBioMass(model,dispOutput)
   % sumBioMass
   %   Calculates breakdown of biomass
   %
-  %   model    (struct) Metabolic model in COBRA format
+  %   model         (struct) Metabolic model in COBRA format
+  %   dispOutput    (bool, opt) If output should be displayed (default = true)
   %
-  %   X         (float) Total biomass fraction [gDW/gDW]
-  %   P         (float) Protein fraction [g/gDW]
-  %   C         (float) Carbohydrate fraction [g/gDW]
-  %   R         (float) RNA fraction [g/gDW]
-  %   D         (float) DNA fraction [g/gDW]
-  %   L         (float) Lipid fraction [g/gDW]
-  %   F         (float) cofactor [g/gDW]
-  %   I         (float) ion [g/gDW]
+  %   X             (float) Total biomass fraction [gDW/gDW]
+  %   P             (float) Protein fraction [g/gDW]
+  %   C             (float) Carbohydrate fraction [g/gDW]
+  %   R             (float) RNA fraction [g/gDW]
+  %   D             (float) DNA fraction [g/gDW]
+  %   L             (float) Lipid fraction [g/gDW]
+  %   F             (float) cofactor [g/gDW]
+  %   I             (float) ion [g/gDW]
   %
-  %   Usage: [X,P,C,R,D,L,I,F] = sumBioMass(model)
+  %   Usage: [X,P,C,R,D,L,I,F] = sumBioMass(model,dispOutput)
   %
   %   Function adapted from SLIMEr: https://github.com/SysBioChalmers/SLIMEr
   %
 
+if nargin < 2
+    dispOutput = true;
+end
+
 %Load original biomass component MWs:
 %TODO: compute MW automatically from chemical formulas (check that all components have them first)
 fid = fopen('../../ComplementaryData/physiology/biomassComposition_Forster2003.tsv');
@@ -40,26 +45,27 @@
 end
 
 %Get main fractions:
-[P,X] = getFraction(model,data,'P',0);
-[C,X] = getFraction(model,data,'C',X);
-[R,X] = getFraction(model,data,'R',X);
-[D,X] = getFraction(model,data,'D',X);
-[L,X] = getFraction(model,data,'L',X);
-[I,X] = getFraction(model,data,'I',X);
-[F,X] = getFraction(model,data,'F',X);
-
-disp(['X -> ' num2str(X) ' gDW/gDW'])
+[P,X] = getFraction(model,data,'P',0,dispOutput);
+[C,X] = getFraction(model,data,'C',X,dispOutput);
+[R,X] = getFraction(model,data,'R',X,dispOutput);
+[D,X] = getFraction(model,data,'D',X,dispOutput);
+[L,X] = getFraction(model,data,'L',X,dispOutput);
+[I,X] = getFraction(model,data,'I',X,dispOutput);
+[F,X] = getFraction(model,data,'F',X,dispOutput);
 
-% Simulate growth:
-sol = optimizeCbModel(model);
-disp(['Growth = ' num2str(sol.f) ' 1/h'])
-disp(' ')
+if dispOutput
+    disp(['X -> ' num2str(X) ' gDW/gDW'])
+    % Simulate growth:
+    sol = optimizeCbModel(model);
+    disp(['Growth = ' num2str(sol.f) ' 1/h'])
+    disp(' ')
+end
 
 end
 
 %%
 
-function [F,X] = getFraction(model,data,compType,X)
+function [F,X] = getFraction(model,data,compType,X,dispOutput)
 
 %Define pseudoreaction name:
 rxnName = [compType ' pseudoreaction'];
@@ -72,7 +78,6 @@
 rxnName = strrep(rxnName,'I','ion');
 rxnName = strrep(rxnName,'F','cofactor');
 
-
 %Add up fraction:
 rxnPos = strcmp(model.rxnNames,rxnName);
 if ~all(rxnPos==0)
@@ -96,11 +101,15 @@
         end
     end
     X = X + F;
-
-    disp([compType ' -> ' num2str(F) ' g/gDW'])
+
+    if dispOutput
+        disp([compType ' -> ' num2str(F) ' g/gDW'])
+    end
 else
-    disp([compType ' do not exist '])
-        F = 0;
+    if dispOutput
+        disp([compType ' does not exist '])
+    end
+    F = 0;
     X = X + F;
 end
 

ComplementaryScripts/saveYeastModel.m

@@ -1,19 +1,26 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% saveYeastModel(model,upDATE)
-% Saves model as a .xml, .txt and .yml file. Also updates complementary
-% files (boundaryMets.txt, README.md and dependencies.txt).
+function saveYeastModel(model,upDATE,allowNoGrowth)
+% saveYeastModel
+%   Saves model as a .xml, .txt and .yml file. Also updates complementary
+%   files (boundaryMets.txt, README.md and dependencies.txt).
+%
+%   Inputs: model           (struct) Model to save (NOTE: must be COBRA format)
+%           upDATE          (bool, opt) If updating the date in the README file
+%                           is needed (default true)
+%           allowNoGrowth   (bool, opt) if saving should be allowed whenever
+%                           the model cannot grow, returning a warning (default
+%                           = true), otherwise will error.
+%   
+%   Usage: saveYeastModel(model,upDATE,allowNoGrowth)
 %
-% model     model structure to save (note: must be in COBRA format)
-% upDATE    logical =true if updating the date in the README file is needed
-%           (opt, default true)
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-function saveYeastModel(model,upDATE)
 
 if nargin < 2
     upDATE = true;
 end
 
+if nargin < 3
+    allowNoGrowth = true;
+end
+
 %Get and change to the script folder, as all folders are relative to this
 %folder
 scriptFolder = fileparts(which(mfilename));
@@ -45,6 +52,10 @@ function saveYeastModel(model,upDATE)
     error('Model should be a valid SBML structure. Please fix all errors before saving.')
 end
 
+%Check if model can grow:
+checkGrowth(model,'aerobic',allowNoGrowth)
+checkGrowth(model,'anaerobic',allowNoGrowth)
+
 %Update .xml, .txt and .yml models:
 copyfile('tempModel.xml','../ModelFiles/xml/yeastGEM.xml')
 delete('tempModel.xml');
@@ -109,4 +120,39 @@ function saveYeastModel(model,upDATE)
 %Switch back to original folder
 cd(currentDir)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+end
+
+%%
+
+function checkGrowth(model,condition,allowNoGrowth)
+%Function that checks if the model can grow or not using COBRA under a
+%given condition (aerobic or anaerobic). Will either return warnings or
+%errors depending on allowNoGrowth.
+
+if strcmp(condition,'anaerobic')
+    cd otherChanges
+    model = anaerobicModel(model);
+    cd ..
+end
+try
+    xPos = strcmp(model.rxnNames,'growth');
+    sol  = optimizeCbModel(model);
+    if sol.v(xPos) < 1e-6
+        dispText = ['The model is not able to support growth under ' ...
+                    condition ' conditions. Please ensure the model can grow'];
+    end
+catch
+    dispText = ['The model yields an infeasible simulation using COBRA ' ...
+                'under ' condition ' conditions. Please ensure the model ' ...
+                'can be simulated with COBRA'];
+end
+
+if exist('dispText','var')
+    if allowNoGrowth
+        warning([dispText ' before opening a PR.'])
+    else
+        error([dispText ' before comitting.'])
+    end
+end
+
+end