-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
v1.2.2
- Loading branch information
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,25 @@ | ||
reaction reactionName A120 A100 A80 A60 G120 G100 G80 G60 X120 X100 X80 X60 | ||
r_1634 acetate exchange -0.512667 -3.172305 -4.343820 -4.020462 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 | ||
r_1654 ammonium exchange -0.030923 -0.247279 -0.389617 -0.452686 -1.507276 -0.993742 -0.936494 -5.786349 -0.002223 -0.036967 -9.244457 -9.478203 | ||
r_1672 carbon dioxide exchange 0.258000 3.215537 4.652600 3.503269 0.807700 2.311900 1.869300 3.000400 0.039452 0.394382 2.488000 3.088000 | ||
r_1718 D-xylose exchange -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.017241 -0.164336 -2.594055 -3.453143 | ||
r_1761 ethanol exchange 0.000000 0.000000 0.000000 0.000000 0.000000 0.136508 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 | ||
r_1793 formate exchange 0.006539 0.000000 0.000000 0.000000 0.000000 0.000000 0.202581 0.000000 0.000000 0.000000 0.012118 0.012554 | ||
r_1808 glycerol exchange -0.000000 -0.000000 -0.000000 -0.000000 -1.003376 -5.226696 -6.285165 -13.274739 -0.000000 -0.000000 -0.000000 -0.000000 | ||
r_1815 glyoxylate exchange 0.186315 0.000000 0.000000 0.000000 0.001652 0.000000 0.207868 0.000000 0.000000 0.000000 0.000000 0.000000 | ||
r_1832 H+ exchange -0.301951 -3.039587 -4.146997 -3.802682 1.105225 0.410663 2.027415 7.103974 -0.000245 -0.002014 8.956486 9.017361 | ||
r_1861 iron(2+) exchange -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 -0.000000 | ||
r_1865 isoamylol exchange 0.000000 0.043541 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 | ||
r_1873 L-alanine exchange 0.000000 0.000000 0.000000 0.000000 0.000000 0.175717 0.000000 4.079272 0.000000 0.000000 0.000000 0.000000 | ||
r_1989 oxaloacetate(2-) exchange 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.499306 0.000000 0.000000 0.000000 0.000000 0.000000 | ||
r_1992 oxygen exchange -0.353445 -2.290040 -3.913066 -2.872298 -1.695843 -3.572956 -5.437861 -9.819611 -0.028052 -0.305925 -6.127170 -6.371272 | ||
r_2005 phosphate exchange -0.002299 -0.016636 -0.022741 -0.024849 -0.008556 -0.047134 -0.047384 -0.060969 -0.000245 -0.002014 -0.014598 -0.023036 | ||
r_2033 pyruvate exchange 0.000000 0.000000 0.000000 0.000000 0.000000 1.435312 1.926621 4.411282 0.000000 0.000000 0.000000 0.000000 | ||
r_2060 sulphate exchange -0.000632 -0.519094 -0.007948 -0.009278 -0.002533 -0.016531 -0.018965 -0.034167 -0.000046 -0.000746 -0.005541 -0.008913 | ||
r_2061 sulphite exchange 0.000000 0.514020 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 | ||
r_2091 urea exchange 0.000000 0.000000 0.000000 0.000000 0.691837 0.000000 0.000000 0.000000 0.000000 0.000000 4.485542 4.520198 | ||
r_2100 water exchange 0.587434 4.270393 6.017051 5.279056 3.743349 11.968339 14.318412 29.461122 0.046558 0.484276 11.788595 12.588360 | ||
r_2104 xylitol exchange 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.001395 0.016699 0.687846 1.109349 | ||
r_2111 growth 0.007818 0.062000 0.086000 0.101764 0.031000 0.162000 0.178000 0.233000 0.000800 0.007000 0.053900 0.090000 | ||
EXC_OUT_s_3717 protein exchange (OUT) 0.000126 0.000981 0.001593 0.001480 0.001228 0.006778 0.006608 0.007255 0.000000 0.000000 0.000000 0.000000 | ||
r_4046 non-growth associated maintenance reaction 0.000002 0.023711 5.820578 0.000034 3.240072 6.293890 13.820184 33.490490 0.105586 1.237592 8.291982 7.669704 |
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Large diffs are not rendered by default.
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,12 @@ | ||
reaction acetate glycerol xylose xylose (no PK-PTA) | ||
acetate exchange -5.000000e-01 0 0 0 | ||
ammonium exchange -1.446275e-03 -1.999921e-03 -1.780501e-03 -1.709859e-03 | ||
carbon dioxide exchange 5.032288e-01 3.130608e-01 3.884278e-01 4.126924e-01 | ||
D-xylose exchange 0 0 -2.000000e-01 -2.000000e-01 | ||
glycerol exchange 0 -3.333333e-01 0 0 | ||
H+ exchange -4.895723e-01 -1.238315e-01 1.283747e-02 1.232813e-02 | ||
oxygen exchange -3.075425e-01 -2.091307e-01 -1.475195e-01 -1.813422e-01 | ||
phosphate exchange -1.678881e-03 -2.321571e-03 -2.066861e-03 -1.984857e-03 | ||
water exchange 5.529470e-01 7.151447e-01 4.496356e-01 4.714716e-01 | ||
lipid backbone exchange 6.035038e-02 8.345300e-02 7.429701e-02 7.134922e-02 | ||
lipid chain exchange 6.028702e-02 8.336538e-02 7.421900e-02 7.127431e-02 |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,13 @@ | ||
,growth,acetate,glycerol,xylose,co2,o2,xylitol,ammonium,protExch,lipidsRatio,proteinRatio | ||
A120,0.009,0.512667333,,,0.258,-2.446423378,,0.011897343,0.000125786,0.534,0.233423761 | ||
A100,0.062,3.172305443,,,4.243,-2.290039723,,0.102895526,0.000980834,0.4,0.238080093 | ||
A80,0.086,4.343820155,,,5.749,-4.576271055,,0.189400418,0.001593176,0.378,0.278774698 | ||
A60,0.114,4.020461624,,,5.12,-3.811746902,,0.241149653,0.001480052,0.272,0.277617088 | ||
G120,0.031,,1.003375652,,0.8077,-4.03836813,,0.051908645,0.001227674,0.439,0.213697615 | ||
G100,0.162,,5.226696569,,2.3119,-3.572956103,,0.419473453,0.006777971,0.523,0.294960672 | ||
G80,0.178,,6.285164983,,1.8693,-5.437861108,,0.464773989,0.006607892,0.39,0.31817701 | ||
G60,0.233,,13.2747394,,3.0004,-9.819610125,,1.304989817,0.007254707,0.366,0.488763919 | ||
X120,0.0008,,,0.017240889,0.103,-2.163,0.001395384,0.000733333,0,0.6,0.1551 | ||
X100,0.007,,,0.164335802,0.615,-2.321,0.01669892,0.008555556,3.367E-09,0.5,0.355409519 | ||
X80,0.0539,,,2.594055279,2.488,-10.343,0.687846423,0.18528125,2.16139E-08,0.412,0.339935599 | ||
X60,0.09,,,3.453142849,3.088,-9.067,1.109348557,0.452844198,0,0.334,0.324461679 |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,177 @@ | ||
clear;clc;if ~exist('scripts') | ~endsWith(scripts,'ComplementaryScripts'); run('../../init_rhtoGEM.m'); end | ||
%% Generate condition-specific models by modifying biomass composition and | ||
% constraining with measured fluxes | ||
model = importModel([root '\ModelFiles\xml\rhto.xml']); | ||
|
||
%% Readjust some reactions | ||
model = setParam(model,'lb','r_4046',0); % NGAM set to zero, maximized later on | ||
|
||
% Add protein excretion reation. | ||
model = addExchangeRxns(model,'out','s_3717'); | ||
lipidPos = find(contains(model.rxnNames,{'lipid backbone pseudoreaction','lipid chain pseudoreaction'})); | ||
[scaleLp,~] = find(model.S(:,lipidPos)<0); | ||
ngamIdx = getIndexes(model,'r_4046','rxns'); | ||
[~,allExIdx] = getExchangeRxns(model,'both'); | ||
allExIdx = [allExIdx; getIndexes(model,'r_4046','rxns')]; % Include NGAM | ||
|
||
%% Load experimental data | ||
expDat = readLopesData(); | ||
cd([scripts '/experimental/']) | ||
|
||
%% Make condition specific models, set rates and adjust biomass equation | ||
clear out | ||
for i=1:length(expDat.sample) | ||
models{i} = model; | ||
%% Set exchange fluxes | ||
models{i} = setParam(models{i},'eq',{'r_1634','r_1808','r_1718','r_1714'},[-expDat.rates(i,2:4),0]); | ||
%% Scale lipids in biomass to fit measured total lipid content | ||
fL = 1.01; | ||
while abs(fL-1) > 0 | ||
fL0 = fL; | ||
models{i}.S(scaleLp,lipidPos) = full(models{i}.S(scaleLp,lipidPos))*fL; | ||
[~,~,~,~,~,L] = sumBioMass(models{i}); | ||
fL = expDat.lipids(i)/L; | ||
end | ||
%% Scale protein & other biomass components | ||
[models{i},GAMpol] = changeOtherComp(models{i},expDat,i); | ||
models{i} = setGAM(models{i},GAMpol); | ||
%% Protein excretion | ||
if ~(expDat.rates(i,9) == 0) | ||
[~,P] = sumBioMass(models{i}); % Weight of protein pseudometabolite | ||
pIdx = getIndexes(models{i},'EXC_OUT_s_3717','rxns'); | ||
models{i}.S(find(models{1}.S(:,pIdx)),pIdx) = -1/P; % Make flux represent 1 g/gDCW/h | ||
models{i} = setParam(models{i},'eq','EXC_OUT_s_3717',expDat.rates(i,9)); | ||
end | ||
%% Xylitol excretion, and limit gas exchanges | ||
models{i} = setParam(models{i},'eq','r_2104',expDat.rates(i,7)); % Xylitol | ||
models{i} = setParam(models{i},'ub','r_1672',expDat.rates(i,5)); % CO2 | ||
models{i} = setParam(models{i},'lb','r_1992',expDat.rates(i,6)); % O2 | ||
models{i} = setParam(models{i},'ub','r_1654',-expDat.rates(i,8)); % Ammonium | ||
% Growth rate | ||
sol(1,i)=solveLP(models{i},1); | ||
sol(1,i)=solveLP(models{i},1); | ||
out(:,i)=sol(1,i).x; | ||
disp(['Growth rate in sample ' num2str(i) ': ' num2str(-sol(1,i).f)]); | ||
end | ||
|
||
%% Compare simulated growth with measured growth. Set the lower value as | ||
% lower bound and then maximize NGAM | ||
for i=1:length(models) | ||
sim = sol(1,i).x(getIndexes(models{i},'r_2111','rxns')); | ||
exp = expDat.rates(i,1); | ||
models{i} = setParam(models{i},'lb','r_2111',min([sim,exp])); | ||
models{i} = setParam(models{i},'obj','r_4046',1); | ||
sol(1,i) = solveLP(models{i},1); | ||
fba(:,i)=sol(1,i).x; | ||
end | ||
%% Attempt to fix CO2 and O2 close to measured value | ||
for i=1:length(models) | ||
disp(['Model: ' num2str(i)]) | ||
models{i} = setParam(models{i},'obj','r_1672',1); | ||
tmp2 = solveLP(models{i}); | ||
models{i} = setParam(models{i},'lb', 'r_1672', 0.999*-tmp2.f); | ||
models{i} = setParam(models{i},'obj','r_1992',-1); | ||
tmp2 = solveLP(models{i}); | ||
models{i} = setParam(models{i},'ub', 'r_1992', 0.999*tmp2.f); | ||
models{i} = setParam(models{i},'obj','r_4046',1); | ||
tmp2 = solveLP(models{i}); | ||
fba(:,i)=tmp2.x; | ||
disp(num2str(-tmp2.f)) | ||
end | ||
|
||
%% Random sampling | ||
% Fix measured exchange fluxes around 5% of the value from FBA | ||
exIdx = getIndexes(model,{'r_1634','r_1808','r_1718','r_1714',... | ||
'EXC_OUT_s_3717','r_2104','r_1672','r_1992','r_1654','r_2111','r_4046'},'rxns'); | ||
|
||
nsamples = 5000;% | ||
%load([root '/scrap/randsampl.mat'], 'goodRxns*') | ||
for i=1:numel(models) | ||
% Fix measured exchange fluxes + NGAM around 5% of the value from FBA. | ||
fluxes = sol(1,i).x(exIdx); | ||
models{i} = setParam(models{i},'var',exIdx,fluxes,5); | ||
if i==1 | ||
[~, goodRxns1] = randomSampling(models{i},1,true,true,true); | ||
elseif i==5 | ||
[~, goodRxns5] = randomSampling(models{i},1,true,true,true); | ||
elseif i==9 | ||
[~, goodRxns9] = randomSampling(models{i},1,true,true,true); | ||
end | ||
if i<5 | ||
rs{i} = randomSampling(models{i},nsamples,true,true,true,goodRxns1); | ||
elseif i<9 | ||
rs{i} = randomSampling(models{i},nsamples,true,true,true,goodRxns5); | ||
else | ||
rs{i} = randomSampling(models{i},nsamples,true,true,true,goodRxns9); | ||
end | ||
end | ||
|
||
for i=1:length(rs) | ||
fluxMean(:,i) = full(mean(rs{i},2)); | ||
fluxSD(:,i) = full(std(rs{i},0,2)); | ||
end | ||
|
||
save([root '/scrap/randsampl.mat'], 'fba', 'rs', 'fluxMean', 'models', 'goodRxns*','fluxSD') | ||
|
||
%% Write files | ||
%FBA results, only exchange fluxes. For internal fluxes, refer to mean from | ||
%random sampling. | ||
out = [models{1}.rxns models{1}.rxnNames num2cell(fba)]; | ||
out = out(allExIdx,:); | ||
rmIdx = find(sum(cell2mat(out(:,3:end)),2) == 0); | ||
out(rmIdx,:) =[]; | ||
fid = fopen([data '/Lopes2019/exp_pFBA_exchFluxes.tsv'],'w'); | ||
fprintf(fid,['%s' repmat('\t%s',1,13) '\n'],["reaction" "reactionName" string(expDat.sample')]); | ||
for j=1:size(out,1) | ||
fprintf(fid,['%s\t%s' repmat('\t%f',1,12) '\n'],out{j,:}); | ||
end | ||
fclose(fid); | ||
|
||
%% Write file with all mean fluxes | ||
clear out;for i=1:numel(models); out(:,i) = fluxMean(:,i); end | ||
clear out2;for i=1:numel(models); out2(:,i) = fluxSD(:,i); end | ||
|
||
out = [models{1}.rxns models{1}.rxnNames num2cell(out) num2cell(out2)]; | ||
fid = fopen([data '/Lopes2019/exp_randSampl_allFluxes.tsv'],'w'); | ||
fprintf(fid,['%s' repmat('\t%s',1,25) '\n'],["rxnID" "reactionName" ... | ||
strcat('mean_',string(expDat.sample')) strcat('std_',string(expDat.sample'))]); | ||
for j=1:length(out) | ||
fprintf(fid,['%s\t%s' repmat('\t%d',1,24) '\n'],out{j,:}); | ||
end | ||
fclose(fid); | ||
|
||
% Write exchange fluxes only | ||
out = out(allExIdx,:); | ||
rmIdx = find(sum(cell2mat(out(:,3:end)),2) == 0); | ||
out(rmIdx,:) =[]; | ||
fid = fopen([data '/Lopes2019/exp_randSampl_exchFluxes.tsv'],'w'); | ||
fprintf(fid,['%s' repmat('\t%s',1,25) '\n'],["rxnID" "reactionName" ... | ||
strcat('mean_',string(expDat.sample')) strcat('std_',string(expDat.sample'))]); | ||
for j=1:length(out) | ||
fprintf(fid,['%s\t%s' repmat('\t%d',1,24) '\n'],out{j,:}); | ||
end | ||
fclose(fid); | ||
|
||
%% NADH and NADPH reactions | ||
for i={'NADPH','NADH'} | ||
[fluxes, rxnIdx] = getMetProduction(models{1},i,fluxMean,true); | ||
clear out | ||
out.rxns = models{1}.rxns(rxnIdx); | ||
out.rxnNames= models{1}.rxnNames(rxnIdx); | ||
out.rxnEqns = constructEquations(models{1},rxnIdx); | ||
out.fluxes = num2cell(fluxes); | ||
out = [out.rxns out.rxnNames out.rxnEqns out.fluxes]; | ||
fid = fopen([data '/Lopes2019/exp_randSampl_' i{1} '_productionFluxes.tsv'],'w'); | ||
fprintf(fid,['%s' repmat('\t%s',1,14) '\n'],["rxnID" "rxnName" ... | ||
"rxnEqn" string(expDat.sample')]); | ||
for j=1:length(out) | ||
fprintf(fid,['%s\t%s\t%s' repmat('\t%f',1,12) '\n'],out{j,:}); | ||
end | ||
fclose(fid); | ||
end | ||
|
||
%% Export models | ||
for i=1:numel(models) | ||
fn = ['rhto_' expDat.sample{i} '.xml']; | ||
exportModel(models{i},[data '/Lopes2019/' fn]); | ||
end |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,12 @@ | ||
function data = readLopesData() | ||
fid = fopen('../../ComplementaryData/data/fluxData_Lopes2019.csv'); | ||
types = textscan(fid,repmat('%s',[1,12]),1,'Delimiter',','); | ||
dat = textscan(fid,['%s' repmat(' %f32 ',[1,11])],'Delimiter',',','HeaderLines',1); | ||
fclose(fid); | ||
data.types = horzcat(types{2:end}); | ||
data.sample = dat{1}; | ||
data.rates = cat(2,dat{[2:10]}); | ||
data.rates(isnan(data.rates)) = 0; | ||
data.lipids = dat{11}; | ||
data.protein = dat{12}; | ||
end |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,60 @@ | ||
clear;clc;if ~exist('scripts') | ~endsWith(scripts,'ComplementaryScripts'); run('../../init_rhtoGEM.m'); end | ||
%% Theoretical simulations | ||
% Take rhtoGEM model as it is | ||
cd([scripts '/experimental/']); | ||
model = importModel([root '\ModelFiles\xml\rhto.xml']); | ||
model.rev(getIndexes(model,{'r_1634','r_1718','r_1808'},'rxns')) = 1; | ||
|
||
% Allow excretion of lipid chains and backbones. This allows for optimizing | ||
% production of lipid pseudometabolite | ||
chainExIdx = getIndexes(model,'r_4064','rxns'); | ||
backbExIdx = getIndexes(model,'r_4062','rxns'); | ||
tmp = setParam(model,'ub',[chainExIdx,backbExIdx],1000); | ||
tmp = setParam(tmp,'obj','r_4062',1); % Maximize lipid pseudometabolite production | ||
tmp = setParam(tmp,'eq','r_4046',0); % Set NGAM to 0 | ||
|
||
% Get optimal solutions for each carbon source | ||
tmp = setParam(tmp,'lb',{'r_1634','r_1808','r_1718','r_1714'},[-1/2,0,0,0]); | ||
solAce = solveLP(tmp,1); | ||
|
||
tmp = setParam(tmp,'lb',{'r_1634','r_1808','r_1718','r_1714'},[0,-1/3,0,0]); | ||
solGly = solveLP(tmp,1); | ||
|
||
tmp = setParam(tmp,'lb',{'r_1634','r_1808','r_1718','r_1714'},[0,0,-1/5,0]); | ||
solXyl = solveLP(tmp,1); | ||
|
||
rxnIdx = find(ismember(model.rxnNames,'phosphoketolase')); | ||
tmp = setParam(tmp,'eq',rxnIdx,0); % Disable phosphoketolase | ||
solXylnoPK = solveLP(tmp,1); | ||
|
||
% Print maximum lipid yields | ||
fprintf(['Lipid pseudometabolite produced per C from each carbon source:\n' ... | ||
'Acetate: ' num2str(-solAce.f) ' Glycerol: ' num2str(-solGly.f) ... | ||
' Xylose: ' num2str(-solXyl.f) ' Xylose (no PK-PTA): ' num2str(-solXylnoPK.f) '\n']) | ||
|
||
% Write file with all fluxes | ||
out = [solAce.x, solGly.x, solXyl.x solXylnoPK.x]; | ||
out = [model.rxnNames num2cell(out)]; | ||
fid = fopen([data '/Lopes2019/theor_maxLipid_allFluxes.tsv'],'w'); | ||
fprintf(fid,'%s\t%s\t%s\t%s\t%s\n',["reaction" "Acetate" "Glycerol" "Xylose" ... | ||
"Xylose (no PK-PTA)"]); | ||
for j=1:length(out) | ||
fprintf(fid,'%s\t%d\t%d\t%d\t%d\n',out{j,:}); | ||
end | ||
fclose(fid); | ||
|
||
%% Write file with exchange fluxes | ||
[~,idx] = getExchangeRxns(model,'both'); | ||
out = out(idx,:); | ||
|
||
% Remove exchange reactions not carrying flux | ||
rmIdx = sum(cell2mat(out(:,2:end)),2) == 0; | ||
out = out(~rmIdx,:); | ||
|
||
fid = fopen([data '/Lopes2019/theor_maxLipid_exchFluxes.tsv'],'w'); | ||
fprintf(fid,['%s' repmat('\t%s',1,4) '\n'],["reaction" "acetate" ... | ||
"glycerol" "xylose" "xylose (no PK-PTA)"]); | ||
for j=1:length(out) | ||
fprintf(fid,['%s' repmat('\t%d',1,4) '\n'],out{j,:}); | ||
end | ||
fclose(fid); |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,18 @@ | ||
function [fluxes, rxnIdx] = getMetProduction(model,metName,solVector,nonZero) | ||
% Which fluxes produce/consume a metabolite (excluding transport | ||
% reactions). | ||
% nonZero only include non zero fluxes, default false | ||
|
||
metLoc = find(ismember(model.metNames,metName)); | ||
subS = full(model.S(metLoc,:)); | ||
[~,rxnIdx] = find(subS); | ||
coeff = sum(subS(:,rxnIdx),1); | ||
|
||
fluxes = solVector(rxnIdx,:) .* transpose(coeff); | ||
|
||
if nonZero | ||
nonZero = sum(fluxes,2)~=0; | ||
fluxes = fluxes(nonZero,:); | ||
rxnIdx = rxnIdx(nonZero); | ||
end | ||
end |