Skip to content

Commit e4507ef

Browse files
edkerkedkerk
committed
WIP
1 parent 09a9e24 commit e4507ef

File tree

6 files changed

+302
-4
lines changed

6 files changed

+302
-4
lines changed

‎code/modelCuration/v9_1_0.m‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,12 @@
1919
% %dataDir=fullfile(pwd(),'..','data','modelCuration','v9.1.0'); % No dataDir required for these curations
2020
cd modelCuration
2121

22+
GAM = 55;
23+
% P = 0.461; %Data from Nissen et al. 1997
24+
% NGAM = 1;
25+
%
26+
% model_an = changeGAM(model_an,GAM,NGAM);
27+
2228
%% ========================================================================
2329
% We blocked MDH2 in anaerobic conditions (see details in the anerobicModel
2430
% script) Experiments suggest that AKG needs to produced inside the
@@ -226,7 +232,6 @@
226232
model.S(:,strcmp(model.rxns,'r_0733')) = -model.S(:,strcmp(model.rxns,'r_0733'));
227233
model = setParam(model,'lb',{'r_0732','r_0733'},0);
228234
model = setParam(model,'rev',{'r_0732','r_0733'},1);
229-
%model = setParam(model,'lb','r_0731',-1000);
230235

231236
% There is no evidence for this PFK1 side reaction in yeast
232237
model = removeReactions(model,'r_0887',true,true,true);

‎code/modelTests/anaerobiosis.m‎

Lines changed: 20 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,36 @@
11
clear; close all
22
% Load the model and apply the corrections for *all* models
3-
%model = loadYeastModel;
3+
model902 = loadYeastModel;
44
cd('../modelCuration/')
55
v9_1_0;
66

77
%% Run growth tests
8-
R2 = growth(model);
8+
R2new = growth(model);
9+
R2old = growthOld(model);
10+
% R decreased from 0.9177 to 0.9077
11+
12+
% Repeat with the model *before* generic 9.1.0 curations were done
13+
R2new902 = growth(model902);
14+
R2old902 = growthOld(model902);
15+
% R decreased from 0.9087 to 0.8985
16+
% So the generic 9.1.0 curations improved, but the new anaerobic function
17+
% actually gives worse predictions
918

1019
%% Convert to anaerobic
1120
cd('../otherChanges/')
1221
modelAn = anaerobicModel(model);
22+
modelAnOld = anaerobicModelOld(model);
23+
modelAn902 = anaerobicModel(model902);
24+
modelAnOld902 = anaerobicModelOld(model902);
1325

1426
cd('../modelTests/');
1527
%% flux predictions
16-
R2=anaerobic_flux_predictions(modelAn);
28+
R2fluxnew = anaerobic_flux_predictions(modelAn);
29+
R2fluxold = anaerobic_flux_predictions(modelAnOld);
30+
R2fluxnew902 = anaerobic_flux_predictions(modelAn902);
31+
R2fluxold902 = anaerobic_flux_predictions(modelAnOld902);
32+
% R
33+
1734

1835
%% Plot
1936
plotAnaerobic(modelAn)

‎code/modelTests/anaerobiosis.mlx‎

225 KB
Binary file not shown.

‎code/modelTests/growthOld.m‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,118 @@
1+
function R2 = growthOld(model_origin,writeOutput)
2+
% This is for growth test: Fig S4c for yeast8 paper
3+
% here we use several chemostat data: 'N-limited aerobic' 'C-limited
4+
% aerobic' 'C-limited anaerobic' 'N-limited anaerobic'
5+
% when simulating N-limited condition, protein content was rescaled, and
6+
% when simulate anaerobic condtion, heme NADH NADP NADPH NAD were rescaled
7+
% to be 0.
8+
9+
funcDir = dbstack('-completenames');
10+
funcDir = regexprep(funcDir(1).file,[funcDir(1).name '\.m'],'');
11+
12+
if nargin<1
13+
model_origin = loadYeastModel;
14+
end
15+
if nargin<2
16+
writeOutput = false;
17+
end
18+
19+
%Load chemostat data:
20+
fid = fopen(fullfile(funcDir,'..','..','data','physiology','chemostatData_Tobias2013.tsv'),'r');
21+
exp_data = textscan(fid,'%f32 %f32 %f32 %f32','Delimiter','\t','HeaderLines',1);
22+
exp_data = [exp_data{1} exp_data{2} exp_data{3} exp_data{4}];
23+
fclose(fid);
24+
exp_data1 = exp_data(1:9,:);
25+
exp_data2 = exp_data(10:20,:);
26+
exp_data3 = exp_data(21:26,:);
27+
exp_data4 = exp_data(27:32,:);
28+
29+
%'N-limited aerobic'
30+
mod_data(1:9,:) = simulateChemostat(model_origin,exp_data(1:9,:),1,'N');
31+
%'C-limited aerobic'
32+
mod_data(10:20,:) = simulateChemostat(model_origin,exp_data(10:20,:),1,'C');
33+
%'C-limited anaerobic'
34+
mod_data(21:26,:) = simulateChemostat(model_origin,exp_data(21:26,:),2,'C');
35+
%'N-limited anaerobic'
36+
mod_data(27:32,:) = simulateChemostat(model_origin,exp_data(27:32,:),2,'N');
37+
38+
% plot the figure
39+
figure
40+
hold on
41+
cols = [215,25,28;253,174,97;171,217,233;44,123,182]/256;
42+
b(1) = plot(exp_data(1:9,4),mod_data(1:9,4),'o','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(2,:));
43+
b(2) = plot(exp_data(10:20,4),mod_data(10:20,4),'s','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(1,:));
44+
b(3) = plot(exp_data(21:26,4),mod_data(21:26,4),'d','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(3,:));
45+
b(4) = plot(exp_data(27:32,4),mod_data(27:32,4),'>','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(4,:));
46+
exp_max = max(exp_data(:,4));
47+
mod_max = max(mod_data(:,4));
48+
lim = max(exp_max,mod_max)+0.05;
49+
xlim([0 lim])
50+
ylim([0 lim])
51+
x=0:0.001:lim;
52+
y = x;
53+
plot(x,y,'--','MarkerSize',6,'Color',[64,64,64]/256)
54+
xlabel('Experimental growth rate [1/h]','FontSize',14,'FontName','Helvetica')
55+
ylabel('In silico growth rate [1/h]','FontSize',14,'FontName','Helvetica')
56+
legend(b,'N-limited aerobic','C-limited aerobic','C-limited anaerobic','N-limited anaerobic','Location','northwest')
57+
58+
% meanerror = sqrt(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)/sqrt(32);
59+
% text(0.25,0.1,['SEM:',num2str(meanerror)])
60+
hold off
61+
R2=corrcoef(exp_data(:,4),mod_data(:,4));
62+
R2=R2(2)^2;
63+
64+
if writeOutput
65+
saveas(gcf,fullfile(funcDir,'..','..','data','testResults','growth.png'));
66+
fid = fopen(fullfile(funcDir,'..','..','data','testResults','growth.md'),'w');
67+
fprintf(fid,'%s\n','## R2 of growth rate prediction');
68+
fprintf(fid,'%.4g\n\n',R2);
69+
fprintf(fid,'%s\n','![Growth curve](growth.png)');
70+
fclose(fid);
71+
end
72+
end
73+
74+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75+
function [mod_data,solresult] = simulateChemostat(model_origin,exp_data,mode1,mode2)
76+
%Relevant positions:
77+
pos(1) = find(strcmp(model_origin.rxns,'r_1714')); %glc
78+
pos(2) = find(strcmp(model_origin.rxns,'r_1992')); %O2
79+
pos(3) = find(strcmp(model_origin.rxns,'r_1654')); %NH3
80+
pos(4) = find(strcmp(model_origin.rxns,'r_2111')); %growth
81+
82+
%Simulate chemostats:
83+
mod_data = zeros(size(exp_data));
84+
solresult = zeros(length(model_origin.rxns),length(exp_data(:,1)));
85+
if mode1 == 2
86+
model_origin = anaerobicModelOld(model_origin);
87+
end
88+
if strcmp(mode2,'N')
89+
% P content in NH3-lim 0.1/h chemostat, per 10.1016/j.femsyr.2005.04.003
90+
model_origin = scaleBioMassOld(model_origin,'protein',0.28,'',false);
91+
% Assume that RNA decreased by the same amount (40%)
92+
model_origin = scaleBioMassOld(model_origin,'RNA',0.0329,'carbohydrate',false);
93+
model_origin = setParam(model_origin,'ub','r_0472',1000); %Glutamate synthase repressed in excess nitrogen
94+
end
95+
for i = 1:length(exp_data(:,1))
96+
model_test= model_origin;
97+
%Fix glucose uptake rate and maximize growth:
98+
for j = 1:length(exp_data(1,:))-1
99+
100+
if abs(exp_data(i,j))==1000
101+
model_test = setParam(model_test,'lb',model_test.rxns(pos(j)),-exp_data(i,j));
102+
else
103+
model_test = setParam(model_test,'eq',model_test.rxns(pos(j)),-exp_data(i,j));
104+
end
105+
end
106+
107+
model_test = setParam(model_test,'obj',model_test.rxns(pos(4)),1);
108+
sol = solveLP(model_test,1);
109+
%Store relevant variables:
110+
try
111+
mod_data(i,:) = abs(sol.x(pos)');
112+
solresult(:,i) = sol.x;
113+
catch
114+
mod_data(i,:) = 0;
115+
solresult(:,i) = 0;
116+
end
117+
end
118+
end

‎code/otherChanges/scaleBioMassOld.m‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
function model = scaleBioMass(model,component,new_value,balance_out,dispOutput)
2+
% scaleBioMass
3+
% Scales the biomass composition
4+
%
5+
% model (struct) metabolic model in COBRA format
6+
% component (str) name of the component to rescale (e.g. "protein")
7+
% new_value (float) new total fraction for said component
8+
% balance_out (str, opt) if chosen, the name of another component with which
9+
% the model will be balanced out so that the total mass remains = 1 g/gDW
10+
% provide empty string '' if this should not be done
11+
% dispOutput (bool, opt) if output from sumBioMass should be displayed (default = true)
12+
%
13+
% model (struct) modified model
14+
%
15+
% Usage: model = scaleBioMass(model,component,new_value,balance_out,dispOutput)
16+
%
17+
18+
if nargin < 5
19+
dispOutput = true;
20+
end
21+
if nargin < 4
22+
balance_out = '';
23+
end
24+
25+
%Measure current composition and rescale:
26+
[~,P,C,R,D,L,I,F] = sumBioMass(model,dispOutput);
27+
content_all = {'carbohydrate','protein','lipid','RNA','DNA','ion','cofactor'};
28+
content_Cap = {'C','P','L','R','D','I','F'};
29+
pos = strcmp(content_all,component);
30+
old_value = eval(content_Cap{pos});
31+
f = new_value / old_value;
32+
model = rescalePseudoReaction(model,component,f);
33+
34+
%Balance out (if desired):
35+
if ~isempty(balance_out)
36+
pos = strcmp(content_all,balance_out);
37+
balance_value = eval(content_Cap{pos});
38+
f = (balance_value - (new_value - old_value)) / balance_value;
39+
model = rescalePseudoReaction(model,balance_out,f);
40+
end
41+
42+
end

‎code/otherChanges/sumBioMassOld.m‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
function [X,P,C,R,D,L,I,F] = sumBioMass(model,dispOutput)
2+
% sumBioMass
3+
% Calculates breakdown of biomass
4+
%
5+
% model (struct) Metabolic model in COBRA format
6+
% dispOutput (bool, opt) If output should be displayed (default = true)
7+
%
8+
% X (float) Total biomass fraction [gDW/gDW]
9+
% P (float) Protein fraction [g/gDW]
10+
% C (float) Carbohydrate fraction [g/gDW]
11+
% R (float) RNA fraction [g/gDW]
12+
% D (float) DNA fraction [g/gDW]
13+
% L (float) Lipid fraction [g/gDW]
14+
% F (float) cofactor [g/gDW]
15+
% I (float) ion [g/gDW]
16+
%
17+
% Usage: [X,P,C,R,D,L,I,F] = sumBioMass(model,dispOutput)
18+
%
19+
% Function adapted from SLIMEr: https://github.com/SysBioChalmers/SLIMEr
20+
%
21+
22+
if nargin < 2
23+
dispOutput = true;
24+
end
25+
26+
%Load original biomass component MWs:
27+
%TODO: compute MW automatically from chemical formulas (check that all components have them first)
28+
fid = fopen('../../data/physiology/biomassComposition_Forster2003.tsv');
29+
Forster2003 = textscan(fid,'%s %s %f32 %f32 %s','Delimiter','\t','HeaderLines',1);
30+
data.mets = Forster2003{1};
31+
data.MWs = double(Forster2003{4});
32+
fclose(fid);
33+
34+
%load additional cofactor/ion MWs:
35+
fid = fopen('../../data/physiology/biomassComposition_Cofactor_Ion.tsv');
36+
CofactorsIons = textscan(fid,'%s %s %f32 %f32 %s %s','Delimiter','\t','HeaderLines',1);
37+
data_new.mets = CofactorsIons{1};
38+
data_new.MWs = double(CofactorsIons{4});
39+
fclose(fid);
40+
for i = 1:length(data_new.mets)
41+
if ~ismember(data_new.mets(i),data.mets)
42+
data.mets = [data.mets; data_new.mets(i)];
43+
data.MWs = [data.MWs; data_new.MWs(i)];
44+
end
45+
end
46+
47+
%Get main fractions:
48+
[P,X] = getFraction(model,data,'P',0,dispOutput);
49+
[C,X] = getFraction(model,data,'C',X,dispOutput);
50+
[R,X] = getFraction(model,data,'R',X,dispOutput);
51+
[D,X] = getFraction(model,data,'D',X,dispOutput);
52+
[L,X] = getFraction(model,data,'L',X,dispOutput);
53+
[I,X] = getFraction(model,data,'I',X,dispOutput);
54+
[F,X] = getFraction(model,data,'F',X,dispOutput);
55+
56+
if dispOutput
57+
disp(['X -> ' num2str(X) ' gDW/gDW'])
58+
% Simulate growth:
59+
sol = solveLP(model,1);
60+
disp(['Growth = ' num2str(sol.f) ' 1/h'])
61+
disp(' ')
62+
end
63+
64+
end
65+
66+
%%
67+
68+
function [F,X] = getFraction(model,data,compType,X,dispOutput)
69+
70+
%Define pseudoreaction name:
71+
rxnName = [compType ' pseudoreaction'];
72+
rxnName = strrep(rxnName,'P','protein');
73+
rxnName = strrep(rxnName,'C','carbohydrate');
74+
rxnName = strrep(rxnName,'N','biomass');
75+
rxnName = strrep(rxnName,'L','lipid backbone');
76+
rxnName = strrep(rxnName,'R','RNA');
77+
rxnName = strrep(rxnName,'D','DNA');
78+
rxnName = strrep(rxnName,'I','ion');
79+
rxnName = strrep(rxnName,'F','cofactor');
80+
81+
%Add up fraction:
82+
rxnPos = strcmp(model.rxnNames,rxnName);
83+
if ~all(rxnPos==0)
84+
isSub = model.S(:,rxnPos) < 0; %substrates in pseudo-rxn
85+
if strcmp(compType,'L')
86+
F = -sum(model.S(isSub,rxnPos)); %g/gDW
87+
else
88+
F = 0;
89+
%Add up all components:
90+
for i = 1:length(model.mets)
91+
pos = strcmp(data.mets,model.mets{i});
92+
if isSub(i) && sum(pos) == 1
93+
if strcmp(compType,'I') || strcmp(compType,'F')
94+
MW = data.MWs(pos);
95+
else
96+
MW = data.MWs(pos)-18;
97+
end
98+
abundance = -model.S(i,rxnPos)*MW/1000;
99+
F = F + abundance;
100+
end
101+
end
102+
end
103+
X = X + F;
104+
105+
if dispOutput
106+
disp([compType ' -> ' num2str(F) ' g/gDW'])
107+
end
108+
else
109+
if dispOutput
110+
disp([compType ' does not exist '])
111+
end
112+
F = 0;
113+
X = X + F;
114+
end
115+
116+
end

0 commit comments

Comments
 (0)