feat/anaerobicgrowth&Nitrogen-limited growth

ComplementaryData/physiology/chemostatData_Tobias2013.tsv

@@ -0,0 +1,33 @@
+GLCxtI	O2xtI	NH3xtI	experimental growth
+5.800000	2.700000	0.400000	0.100000
+4.830000	4.420000	0.420000	0.100000
+3.500000	7.800000	0.610000	0.150000
+4.610000	9.200000	0.740000	0.180000
+5.300000	1000	0.850000	0.200000
+5.670000	8.700000	0.830000	0.200000
+8.000000	8.800000	0.960000	0.250000
+9.450000	9.300000	1.090000	0.280000
+12.680000	8.200000	1.330000	0.340000
+1.150000	2.700000	1000	0.100000
+1.170000	2.500000	1000	0.100000
+1.690000	4.000000	1000	0.150000
+2.260000	5.000000	1000	0.200000
+2.880000	6.500000	1000	0.250000
+3.270000	7.460000	1000	0.270000
+3.290000	7.800000	1000	0.280000
+3.880000	8.000000	1000	0.310000
+6.200000	7.000000	1000	0.330000
+7.890000	6.500000	1000	0.350000
+13.390000	3.000000	1000	0.380000
+2.300000	0.000000	1000	0.030000
+2.860000	0.000000	1000	0.047000
+5.880000	0.000000	1000	0.101000
+11.310000	0.000000	1000	0.190000
+16.750000	0.000000	1000	0.281000
+22.180000	0.000000	1000	0.369000
+4.150000	0.000000	0.230000	0.050000
+8.400000	0.000000	0.480000	0.100000
+8.700000	0.000000	0.550000	0.100000
+12.400000	0.000000	0.920000	0.160000
+15.300000	0.000000	1.330000	0.200000
+17.400000	0.000000	1.620000	0.240000

ComplementaryScripts/modelCuration/scaleBioMass.m → ComplementaryScripts/modelCuration/fixBiomassComp.m

@@ -1,5 +1,5 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% scaleBioMass
+% fixBiomassComp
 % Corrects the stoichiometry coefficients of all pseudo-rxns in an
 % iterative fashion:
 % 1. Switch back to original model's abundance values (Forster et al. 2003)
@@ -8,19 +8,19 @@
 % 3. Rescale carbohydrate fraction (total) to have biomass add up to 1
 % 4. GAM is fitted to simulate chemostat data of S. cerevisiae at low
 %    growth rates (Van Hoek et al. 1988)
-% 
+%
 % Function adapted from SLIMEr: https://github.com/SysBioChalmers/SLIMEr
 %
 % Benjamin Sanchez. Last update: 2018-09-07
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function scaleBioMass
+function fixBiomassComp
 
 %Load model:
 initCobraToolbox
 cd ..
 model = loadYeastModel;
-cd modelCuration
+cd otherChanges
 
 %Load data from Forster 2003:
 fid = fopen('../../ComplementaryData/physiology/biomassComposition_Forster2003.tsv');
@@ -32,7 +32,7 @@
 fclose(fid);
 
 %Compute current composition:
-sumBioMass(model,data);
+sumBioMass(model);
 
 %Switch to Forster 2003 data:
 for i = 1:length(data.mets)
@@ -54,7 +54,7 @@
 end
 
 %Compute current composition:
-sumBioMass(model,data);
+sumBioMass(model);
 
 %Correct with data from biomassComposition_Cofactor_Ion:
 data_original = data;
@@ -66,7 +66,9 @@
 data.groups     = Cofactors{5};
 fclose(fid);
 
+cd ../modelCuration
 model = addBiomassUpdate(model,data);
+cd ../otherChanges
 
 for j = 1:length(data_original.mets)
     if ~ismember(data_original.mets(j),data.mets)
@@ -76,7 +78,7 @@
         data.groups = [data.groups; data_original.groups(j)];
     end
 end
-[X,P,C,R,~,~,~,~] = sumBioMass(model,data);
+[X,P,C,R,~,~,~,~] = sumBioMass(model);
 
 %Correct with data from Lahtvee 2017:
 fid = fopen('../../ComplementaryData/physiology/biomassComposition_Lahtvee2017.tsv');
@@ -92,11 +94,11 @@
     if strcmp(metName,'protein')      %protein fraction
         fP    = abundance/P;        %ratio to scale
         model = rescalePseudoReaction(model,'protein',fP);
-        
+
     elseif strcmp(metName,'RNA')      %RNA fraction
         fR    = abundance/R;        %ratio to scale
         model = rescalePseudoReaction(model,'RNA',fR);
-        
+
     else    %Some extra carbohydrates:
         modelPos = strcmp(model.mets,metID);
         compPos  = strcmp(data.mets,metID);
@@ -107,39 +109,21 @@
         end
     end
 end
-[X,P,C,R,~,~,~,~] = sumBioMass(model,data);
+[X,P,C,R,~,~,~,~] = sumBioMass(model);
 
 %Balance out mass with carbohydrate content:
 delta = X - 1;           %difference to balance
 fC    = (C - delta)/C;
 model = rescalePseudoReaction(model,'carbohydrate',fC);
-sumBioMass(model,data);
+sumBioMass(model);
 
 %Fit GAM:
 model = fitGAM(model);
-sumBioMass(model,data);
+sumBioMass(model);
 
 %Finally, save model:
 cd ..
 saveYeastModel(model)
 cd modelCuration
 
 end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-function model = rescalePseudoReaction(model,metName,f)
-
-rxnName = [metName ' pseudoreaction'];
-rxnPos  = strcmp(model.rxnNames,rxnName);
-for i = 1:length(model.mets)
-    S_ir   = model.S(i,rxnPos);
-    isProd = strcmp(model.metNames{i},[metName ' [cytoplasm]']);
-    if S_ir ~= 0 && ~isProd
-        model.S(i,rxnPos) = f*S_ir;
-    end
-end
-
-end
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ComplementaryScripts/modelTests/growth.m

@@ -0,0 +1,97 @@
+%This is for growth test: Fig S4c for yeast8 paper
+% here we use several chemostat data: 'N-limited aerboic' 'C-limited
+% aerobic' 'C-limited anaerobic' 'N-limited anaerobic'
+% when simulating N-limited condition, protein content was rescaled, and
+% when simulate anaerobic condtion, heme NADH NADP NADPH NAD were rescaled
+% to be 0.
+% Feiran Li -2018-08-25
+% Feiran Li -Last update: 2019-09-24
+
+initCobraToolbox
+cd ..
+model = loadYeastModel;
+model_origin = model;
+cd otherChanges/
+
+%Load chemostat data:
+fid = fopen('../../ComplementaryData/physiology/chemostatData_Tobias2013.tsv','r');
+exp_data = textscan(fid,'%f32 %f32 %f32  %f32','Delimiter','\t','HeaderLines',1);
+exp_data = [exp_data{1} exp_data{2} exp_data{3} exp_data{4}];
+fclose(fid);
+%'N-limited aerboic'
+exp_data1 = exp_data(1:9,:);
+%'C-limited aerobic'
+exp_data2 = exp_data(10:20,:);
+%'C-limited anaerobic'
+exp_data3 = exp_data(21:26,:);
+%'N-limited anaerobic'
+exp_data4 = exp_data(27:32,:);
+
+mod_data1 = simulateChemostat(model_origin,exp_data1,1,'N');
+mod_data2 = simulateChemostat(model_origin,exp_data2,1,'C');
+mod_data3 = simulateChemostat(model_origin,exp_data3,2,'C');
+mod_data4 = simulateChemostat(model_origin,exp_data4,2,'N');
+
+cd ../modelTests/
+% plot the figure
+figure
+hold on
+cols = [215,25,28;253,174,97;171,217,233;44,123,182]/256;
+b(1) = plot(exp_data1(:,4),mod_data1(:,4),'o','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(2,:));
+b(2) = plot(exp_data2(:,4),mod_data2(:,4),'s','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(1,:));
+b(3) = plot(exp_data3(:,4),mod_data3(:,4),'d','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(3,:));
+b(4) = plot(exp_data4(:,4),mod_data4(:,4),'>','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(4,:));
+exp_max = max(exp_data2(:,4));
+mod_max = max(mod_data1(:,4));
+lim = max(exp_max,mod_max)+0.05;
+xlim([0 lim])
+ylim([0 lim])
+x=0:0.001:lim;
+y = x;
+plot(x,y,'--','MarkerSize',6,'Color',[64,64,64]/256)
+xlabel('Experimental growth rate [1/h]','FontSize',14,'FontName','Helvetica')
+ylabel('In silico growth rate [1/h]','FontSize',14,'FontName','Helvetica')
+legend(b,'N-limited aerboic','C-limited aerobic','C-limited anaerobic','N-limited anaerobic','Location','northwest')
+meanerror = sum(([exp_data1(:,4);exp_data2(:,4);exp_data3(:,4);exp_data4(:,4)]-[mod_data1(:,4);mod_data2(:,4);mod_data3(:,4);mod_data4(:,4)]).^2)/32;
+text(0.4,0.1,['meanerror:',num2str(meanerror*100),'%'])
+hold off
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [mod_data,solresult] = simulateChemostat(model_origin,exp_data,mode1,mode2)
+model = model_origin;
+%Relevant positions:
+pos(1) = find(strcmp(model.rxns,'r_1714')); %glc
+pos(2) = find(strcmp(model.rxns,'r_1992')); %O2
+pos(3) = find(strcmp(model.rxns,'r_1654')); %NH3
+pos(4) = find(strcmp(model.rxns,'r_2111')); %growth
+
+%Simulate chemostats:
+mod_data = zeros(size(exp_data));
+solresult = zeros(length(model.rxns),length(exp_data(:,1)));
+if strcmp(mode2,'N')
+    model = scaleBioMass(model,'protein',0.289);
+    model = scaleBioMass(model,'lipid',0.048);
+    model = scaleBioMass(model,'RNA',0.077,'carbohydrate');
+end
+if mode1 == 2
+    model = anaerobicModel(model);
+end
+for i = 1:length(exp_data(:,1))
+    model_test= model;
+    %Fix glucose uptake rate and maxmize growth:
+    for j = 1:length(exp_data(1,:))-1
+
+        if abs(exp_data(i,j))==1000
+            model_test = changeRxnBounds(model_test,model_test.rxns(pos(j)),-exp_data(i,j),'l');
+        else
+            model_test = changeRxnBounds(model_test,model_test.rxns(pos(j)),-exp_data(i,j),'b');
+        end
+    end
+
+    model_test = changeObjective(model_test,model_test.rxns(pos(4)),+1);
+    sol        = optimizeCbModel(model_test,'max');
+    %Store relevant variables:
+    mod_data(i,:) = abs(sol.x(pos)');
+    solresult(:,i) = sol.x;
+end
+end

ComplementaryScripts/otherChanges/anaerobicModel.m

@@ -2,10 +2,28 @@
 % anaerobicModel.m
 % Converts model to anaerobic
 %
-% Benjam�n J. S�nchez
+% Benjamin J. Sanchez
+% Feiran Li - 2019-09-24
+% Feiran Li - Last update: 2019-10-02 modify the order of changes
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-%1st change: Changes media to anaerobic (no O2 uptake and allows sterol
+function model = anaerobicModel(model)
+
+%1th 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');
+
+%2nd change: Removes the requirement of heme a in the biomass equation
+%            (not used under aerobic conditions)
+mets = {'s_3714[c]','s_1198[c]','s_1203[c]','s_1207[c]','s_1212[c]'};
+[~,met_index] = ismember(mets,model.mets);
+model.S(met_index,strcmp(model.rxns,'r_4598')) = 0;
+
+%3st change: Changes media to anaerobic (no O2 uptake and allows sterol
 %            and fatty acid exchanges)
 model.lb(strcmp(model.rxns,'r_1992')) = 0;        %O2
 model.lb(strcmp(model.rxns,'r_1757')) = -1000;    %ergosterol
@@ -16,11 +34,7 @@
 model.lb(strcmp(model.rxns,'r_2137')) = -1000;    %ergosta-5,7,22,24(28)-tetraen-3beta-ol
 model.lb(strcmp(model.rxns,'r_2189')) = -1000;    %oleate
 
-%2nd change: Removes the requirement of heme a in the biomass equation
-%            (not used under aerobic conditions)
-model.S(strcmp(model.mets,'s_3714[c]'),strcmp(model.rxns,'r_4041')) = 0;
-
-%3rd change: Blocked pathways for proper glycerol production
+%4rd change: Blocked pathways for proper glycerol production
 %Block oxaloacetate-malate shuttle (not present in anaerobic conditions)
 model.lb(strcmp(model.rxns,'r_0713')) = 0; %Mithocondria
 model.lb(strcmp(model.rxns,'r_0714')) = 0; %Cytoplasm
@@ -29,12 +43,25 @@
 %Block 2-oxoglutarate + L-glutamine -> 2 L-glutamate (alternative pathway)
 model.ub(strcmp(model.rxns,'r_0472')) = 0;
 
-%4th 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
-cd ../modelCuration
-model = changeBiomass(model,P,GAM,NGAM);
-cd ../otherChanges
+end
+
+%%
+
+function model = changeGAM(model,GAM,NGAM)
+
+bioPos = strcmp(model.rxnNames,'biomass pseudoreaction');
+for i = 1:length(model.mets)
+    S_ix  = model.S(i,bioPos);
+    isGAM = sum(strcmp({'ATP [cytoplasm]','ADP [cytoplasm]','H2O [cytoplasm]', ...
+        'H+ [cytoplasm]','phosphate [cytoplasm]'},model.metNames{i})) == 1;
+    if S_ix ~= 0 && isGAM
+        model.S(i,bioPos) = sign(S_ix)*GAM;
+    end
+end
+
+if nargin >1
+    pos = strcmp(model.rxnNames,'non-growth associated maintenance reaction');%NGAM
+    model = setParam(model,'eq',model.rxns(pos),NGAM);% set both lb and ub to be mu
+end
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+end

ComplementaryScripts/otherChanges/rescalePseudoReaction.m

@@ -0,0 +1,29 @@
+function model = rescalePseudoReaction(model,metName,f)
+  % rescalePseudoReaction
+  %   Rescales a specific pseudoreaction by a given factor
+  %
+  %   model      (struct) the yeast GEM
+  %   metName    (str) name of the component to rescale (e.g. "protein")
+  %   f          (float) fraction to use for rescaling
+  %
+  %   model    (struct) the (rescaled) yeast GEM
+  %
+  %   Usage: model = rescalePseudoReaction(model,metName,f)
+  %
+
+if strcmp(metName,'lipid')
+    model = rescalePseudoReaction(model,'lipid backbone',f);
+    model = rescalePseudoReaction(model,'lipid chain',f);
+else
+    rxnName = [metName ' pseudoreaction'];
+    rxnPos  = strcmp(model.rxnNames,rxnName);
+    for i = 1:length(model.mets)
+        S_ir   = model.S(i,rxnPos);
+        isProd = strcmp(model.metNames{i},[metName ' [cytoplasm]']);
+        if S_ir ~= 0 && ~isProd
+            model.S(i,rxnPos) = f*S_ir;
+        end
+    end
+end
+
+end

ComplementaryScripts/otherChanges/scaleBioMass.m

@@ -0,0 +1,33 @@
+function model = scaleBioMass(model,component,new_value,balance_out)
+  % scaleBioMass
+  %   Scales the biomass composition
+  %
+  %   model          (struct) metabolic model in COBRA format
+  %   component      (str) name of the component to rescale (e.g. "protein")
+  %   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
+  %
+  %   model          (struct) modified model
+  %
+  %   Usage: model = scaleBioMass(model,component,new_value,balance_out)
+  %
+
+%Measure current composition and rescale:
+[~,P,C,R,D,L,I,F] = sumBioMass(model);
+content_all = {'carbohydrate','protein','lipid','RNA','DNA','ion','cofactor'};
+content_Cap = {'C','P','L','R','D','I','F'};
+pos         = strcmp(content_all,component);
+old_value   = eval(content_Cap{pos});
+f           = new_value / old_value;
+model       = rescalePseudoReaction(model,component,f);
+
+%Balance out (if desired):
+if nargin > 3
+    pos           = strcmp(content_all,balance_out);
+    balance_value = eval(content_Cap{pos});
+    f             = (balance_value - (new_value - old_value)) / balance_value;
+    model         = rescalePseudoReaction(model,balance_out,f);
+end
+
+end