Skip to content

Commit 09a9e24

Browse files
edkerkedkerk
committed
fix: rxn balance curations and script refactor
1 parent 33cf2dc commit 09a9e24

File tree

4 files changed

+124
-68
lines changed

4 files changed

+124
-68
lines changed

‎code/modelCuration/v9_1_0.m‎

Lines changed: 87 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -61,9 +61,9 @@
6161
% HcytIdx = getIndexes(model,'s_0794','mets'); % H+[c]
6262
% HextIdx = getIndexes(model,'s_0796','mets'); % H+[e]
6363
%
64-
% symporterIDs = find(model.S(HcytIdx,:) & model.S(HextIdx,:));
64+
% symporterIDs = transpose(find(model.S(HcytIdx,:) & model.S(HextIdx,:)));
6565
% for i = 1:length(symporterIDs)
66-
% if strcmp(model.rxns(symporterIDs(i)), 'r_1258')
66+
% if ismember(model.rxns(symporterIDs(i)), {'r_1258'})
6767
% % Ignore the sodium transporter, without it, the model does not work
6868
% continue
6969
% end
@@ -125,20 +125,76 @@
125125
model.metCharges(strcmp(model.mets,'s_3906'))=-1;
126126
model.metCharges(strcmp(model.mets,'s_4263'))=-1;
127127

128+
% Balance the reactions 'r_0774' and 'r_0775', 'NAPRtase' by removing H+
129+
% consumption and adding a H2O as a reactant
130+
model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_0774'))=0; % Cytosolic
131+
model.S(find(strcmp(model.mets,'s_0803')),strcmp(model.rxns,'r_0774'))=-1; % Cytosolic
132+
model.S(find(strcmp(model.mets,'s_0799')),strcmp(model.rxns,'r_0775'))=0; % Mitochondrial
133+
model.S(find(strcmp(model.mets,'s_0807')),strcmp(model.rxns,'r_0775'))=-1; % Mitochondrial
134+
135+
% Balance the reaction r_0721, 'malonyl-CoA-ACP transacylase' by adding a
136+
% proton as reactant
137+
model.S(find(strcmp(model.mets,'s_0799')),strcmp(model.rxns,'r_0721'))=-1;
138+
139+
% Balance the reaction r_1603, '4-amino-5-hydroxymethyl-2-methylpyrimidine synthetase'
140+
% by adding a proton as reactant
141+
model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_1603'))=-1;
142+
143+
% Balance the reactions r_2142-2144, 'B-ketoacyl-ACP synthase' by removing
144+
% proton as reactant
145+
model.S(find(strcmp(model.mets,'s_0799')),ismember(model.rxns,{'r_2142','r_2143','r_2144'}))=0;
146+
147+
% Balance the reactions r_2232, 'peroxisomal acyl-CoA thioesterase (4:0)'
148+
% by correcting H+
149+
model.S(find(strcmp(model.mets,'s_0801')),strcmp(model.rxns,'r_2232')) = 1;
150+
151+
% Correct product and reactant of r_2236 and r_2254, part of peroxisomal
152+
% beta-oxidation, where the intermediate metabolite should be
153+
% trans-but-2-enoyl-CoA, not but-2-enoyl-CoA
154+
metsToAdd.metNames = 'trans-but-2-enoyl-CoA';
155+
metsToAdd.compartments = 'p';
156+
metsToAdd.metSmiles = 'C/C=C/C(=O)SCCNC(=O)CCNC(=O)[C@@H](C(C)(C)COP(=O)(O)OP(=O)(O)OC[C@@H]1[C@H]([C@H]([C@@H](O1)N2C=NC3=C(N=CN=C32)N)O)OP(=O)(O)O)O';
157+
metsToAdd.metFormulas = 'C25H36N7O17P3S';
158+
metsToAdd.metCharges = -4;
159+
metsToAdd.metMiriams{1} = struct('name',{{'chebi';'metanetx.chemical'}},...
160+
'value',{{'CHEBI:50998';'MNXM1364409'}});
161+
162+
model = addMets(model,metsToAdd,false,'s_');
163+
model = removeMets(model, {'but-2-enoyl-CoA'}, true, true, true, true);
164+
model.S(end,ismember(model.rxns,'r_2236')) = 1;
165+
model.S(end,ismember(model.rxns,'r_2254')) = -1;
166+
model.S(find(strcmp(model.mets,'s_0801')),ismember(model.rxns,'r_2236')) = 0;
167+
model.S(find(strcmp(model.mets,'s_0801')),ismember(model.rxns,'r_2254')) = 0;
168+
model.S(find(strcmp(model.mets,'s_0801')),ismember(model.rxns,'r_2284')) = +4;
169+
170+
% Balance the reaction r_4629, 'alcohol acyltransferase (hexanoyl-CoA)'
171+
% by adding a proton as product
172+
model.S(find(strcmp(model.mets,'s_0799')),strcmp(model.rxns,'r_4629')) = +4;
173+
174+
% Balance the reaction r_4322, 'GPI mannosyltransferase 4'
175+
% by removing proton
176+
model.S(find(strcmp(model.mets,'s_0795')),strcmp(model.rxns,'r_4322')) = 0;
177+
178+
% Balance the reaction r_4679, 'short-chain-fatty-acid-CoA ligase (propionate)'
179+
% by adding a proton as product
180+
model.S(find(strcmp(model.mets,'s_0801')),strcmp(model.rxns,'r_4679')) = 1;
181+
182+
% Balance the reaction r_4701, 'L-cysteine hydrogen-sulfide-lyase'
183+
% by adding a proton as product
184+
model.S(find(strcmp(model.mets,'s_0799')),strcmp(model.rxns,'r_4701')) = 1;
185+
128186
% Balance the reaction r_4702, 'L-cysteine:2-oxoglutarate aminotransferase'
129187
% by adding a proton as reactant
130188
model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4702'))=-1;
131189

132-
% Balance the reaction r_4703, 'L-cysteine:2-oxoglutarate aminotransferase'
190+
% Balance the reaction r_4703, '3-mercaptopyruvate sulfurtransferase'
133191
% by adding a proton as product
134192
model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4703'))=1;
135193

136-
% Balance the reactions 'r_0774' and 'r_0775', 'NAPRtase' by removing H+
137-
% consumption and adding a H2O as a reactant
138-
model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_0774'))=0;
139-
model.S(find(strcmp(model.mets,'s_0803')),strcmp(model.rxns,'r_0774'))=-1;
140-
model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_0775'))=0;
141-
model.S(find(strcmp(model.mets,'s_0803')),strcmp(model.rxns,'r_0775'))=-1;
194+
% Balance the reaction r_4707, 'trithionate thiosulfohydrolase'
195+
% by adding a proton as product
196+
model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4707'))=1;
197+
142198

143199
%% ========================================================================
144200
% This section focuses on individual reactions that have the wrong
@@ -147,25 +203,33 @@
147203
% Make GCY1 irreversible. Has a positive DeltaGo' (+20.9) and is part of a
148204
% transhydrogenase cycle (NADH -> NADPH) at the cost of one ATP. High
149205
% cytosolic NADPH/NADP ratio makes it thermodynamically infeasible that it
150-
% runs in reverse direction.
151-
model = setParam(model,'ub','r_0487',0);
206+
% runs in reverse direction. First swap direction, then make irreversible
207+
model.S(:,strcmp(model.rxns,'r_0487')) = -model.S(:,strcmp(model.rxns,'r_0487'));
208+
model = setParam(model,'lb','r_0487',0);
209+
model = setParam(model,'rev','r_0226',0);
152210

153211
% The mitochondrial ATP synthase is able to run in reverse, which occurs
154212
% in anaerobic conditions
155213
model = setParam(model,'lb','r_0226',-1000);
214+
model = setParam(model,'rev','r_0226',1);
156215

157216
% Rename r_0227, it is the plasma membrane ATPase, not a cytosolic ATPase
158217
model.rxnNames(strcmp(model.rxns,'r_0227')) = {'ATPase, plasma membrane'};
159218

160219
% Make sure both formate-THF ligases are reversible (r_0447 already is).
161220
model = setParam(model,'lb','r_0446',-1000);
221+
model = setParam(model,'rev','r_0446',1);
162222

163223
% Make methylenetetrahydrofolate dehydrogenases ADE3 and MIS1 irreversible
164-
model = setParam(model,'ub',{'r_0732','r_0733'},0);
224+
% First swap direction, then make irreversible
225+
model.S(:,strcmp(model.rxns,'r_0732')) = -model.S(:,strcmp(model.rxns,'r_0732'));
226+
model.S(:,strcmp(model.rxns,'r_0733')) = -model.S(:,strcmp(model.rxns,'r_0733'));
227+
model = setParam(model,'lb',{'r_0732','r_0733'},0);
228+
model = setParam(model,'rev',{'r_0732','r_0733'},1);
229+
%model = setParam(model,'lb','r_0731',-1000);
165230

166-
% There is no evidence for this PFK1 side reaction in yeast. Consider
167-
% removing the reaction completley
168-
model = setParam(model,'eq','r_0887',0);
231+
% There is no evidence for this PFK1 side reaction in yeast
232+
model = removeReactions(model,'r_0887',true,true,true);
169233

170234
% TYR1 incorrectly annotated as using NAD, should be NADP
171235
model.S(find(strcmp(model.mets,'s_1212')),strcmp(model.rxns,'r_0939'))=0; %NADPH
@@ -180,7 +244,6 @@
180244
% Make polyphosphate hydrolase and diphosphate transport over cell membrane
181245
% both irreversible
182246
model = setParam(model,'lb',{'r_4723','r_4724','r_4725'},0);
183-
model = setParam(model,'lb','r_4460',0);
184247

185248
% While r_0013 was elementary balanced, it was not charged balanced. The
186249
% reaction mechanism was incorrect. Corrected to mimic a combination of
@@ -203,16 +266,6 @@
203266
model.metFormulas(idx) = {'CH2O3S'};
204267
model.metCharges(idx) = -1;
205268

206-
207-
% [~,metFormulae] = computeMetFormulae(model,'metMwRange','s_0338','fillMets','none')
208-
% model.metFormulas(getIndexes(model,'s_0329','mets')) = {'C17H28NO17P'};
209-
% model.metFormulas(getIndexes(model,'s_0330','mets')) = {'C33H58NO18P'};
210-
% model.metFormulas(getIndexes(model,'s_0331','mets')) = {'C19H30NO18P'};
211-
% model.metFormulas(getIndexes(model,'s_0334','mets')) = {'C39H68NO23P'};
212-
% model.metFormulas(getIndexes(model,'s_0337','mets')) = {'C44H78N2O27P2'};
213-
% model.metFormulas(getIndexes(model,'s_0338','mets')) = {'C41H38N2O41P3'};
214-
% model.metFormulas(getIndexes(model,'s_0339','mets')) = {'C50H88N2O32P2'};
215-
216269
% Copy annotations between same metabolites in separate compartments
217270
model.metFormulas(getIndexes(model,'s_4211','mets')) = model.metFormulas(getIndexes(model,'s_2885','mets'));
218271
model.metCharges(getIndexes(model,'s_4211','mets')) = model.metCharges(getIndexes(model,'s_2885','mets'));
@@ -227,11 +280,18 @@
227280
% metabolic process, and will therefore be removed
228281
model = removeReactions(model,'r_4325',true,true,true);
229282

283+
% r_4568 is an unbalanced, dead-end, non-gene-associated reaction without
284+
% known reaction mechanism
285+
model = removeReactions(model,'r_4568',true,true,true);
286+
287+
% r_4704, r_4706 and r_4709 revolve around the unspecific metabolite
288+
% alkanesulfonate, with unbalanced r_4706
289+
model = removeReactions(model,{'r_4704','r_4706','r_4709'},true,true,true);
290+
230291
% r_0229
231292
%model = changeRxns(model,'r_0229','dethiobiotin[c] + polysulphur[c] <=> biotin[c] + 2 H+[c]',2)
232293
%dethiobiotin[c] + hydrogen sulfide[c] + 2 S-adenosyl-L-methionine[c] + 2 H+[c] <=> biotin[c] + 2 L-methionine[c] + 2 5'-Deoxyadenosine
233294

234-
235295
%% ========================================================================
236296
% Condition-specific gene expression. These can be enabled with scripts
237297
% Glycine cleavage only active when glycine is used as nitrogen source

‎code/modelTests/anaerobic_flux_predictions.m‎

Lines changed: 16 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@
77

88
sim_vals=[];
99

10-
1110
colors = orderedcolors("glow12");
1211

1312
figure;
@@ -17,31 +16,31 @@
1716
merged_names=[];
1817

1918
for i=1:length(data_sets)
20-
21-
index_data_set=find(strcmp(text_flux(:,6),data_sets(i,1)));
22-
index_glx=strcmp(text_flux(index_data_set,7),'r_1714');
23-
24-
model = setParam(model,'eq','r_1714',-cell2mat(vals_flux(index_data_set(find(index_glx)),4) ));
25-
26-
%% Solve the LP problem
19+
% Gather single data set
20+
idx_data_set = find(strcmp(text_flux(:,6),data_sets(i,1)));
21+
idx_glc = idx_data_set(strcmp(text_flux(idx_data_set,7),'r_1714'));
22+
% Set glucose uptake
23+
model = setParam(model,'eq','r_1714',-cell2mat(vals_flux(idx_glc,4)));
24+
% Solve the LP problem
2725
res=solveLP(model,1);
28-
rxns = text_flux(index_data_set,7);
29-
rxns(cellfun(@isempty,rxns)) = [];
30-
index_model=getIndexes(model,rxns,'rxns');
31-
exclude_data=(index_model==0);
32-
index_model(exclude_data)=[];
33-
index_data_set(exclude_data)=[];
26+
% Organize data and output
27+
rxns = text_flux(idx_data_set,7);
28+
include_data = ismember(rxns,model.rxns);
3429

35-
scaled_sim=abs(-100.*res.x(index_model)./res.x(getIndexes(model,'r_1714','rxns')));
30+
rxns(~include_data) = [];
31+
idx_data_set(~include_data) = [];
32+
33+
idx_model = getIndexes(model,rxns,'rxns');
34+
35+
scaled_sim=abs(-100.*res.x(idx_model)./res.x(getIndexes(model,'r_1714','rxns')));
3636

37-
data_vals=abs(cell2mat(vals_flux(index_data_set,5)));
37+
data_vals=abs(cell2mat(vals_flux(idx_data_set,5)));
3838
merged_data=[merged_data data_vals'];
3939
merged_sim=[merged_sim scaled_sim'];
4040
merged_names=[merged_names rxns'];
4141
plot(data_vals,scaled_sim,'^','MarkerFaceColor',colors(i,:),'MarkerEdgeColor',colors(i,:));
4242
hold on;
4343

44-
% res.x()
4544
end
4645

4746

‎code/modelTests/anaerobiosis.m‎

Lines changed: 14 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -21,39 +21,32 @@
2121
modelAn = setParam(modelAn,'eq','r_1714',-23);
2222
res=solveLP(modelAn,1);
2323
FLUX = res.x;
24-
v_AStr = res.x(getIndexes(modelAn,'r_1115','rxns'))
25-
v_ATPase = res.x(getIndexes(modelAn,'r_0227','rxns'))
26-
v_glc = res.x(getIndexes(modelAn,'r_1714','rxns'))
24+
v_AStr = res.x(getIndexes(modelAn,'r_1115','rxns'));
25+
v_ATPase = res.x(getIndexes(modelAn,'r_0227','rxns'));
26+
v_glc = res.x(getIndexes(modelAn,'r_1714','rxns'));
2727

2828
%% Pack flux results into table
2929
temp_model=modelAn;
30-
temp_model.metNames=strcat( modelAn.metNames, repmat('[',length(modelAn.mets),1),modelAn.compNames(temp_model.metComps),repmat(']',length(modelAn.mets),1) );
31-
rxns_reacs=printRxnFormula(temp_model,'rxnAbbrList',modelAn.rxns,'metNameFlag',1,'printFlag',0);
30+
rxns_reacs=constructEquations(temp_model,modelAn.rxns);
3231
tab=table(modelAn.rxns,modelAn.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,modelAn.grRules);
33-
%tab = sortrows(tab,"Var4","descend");
3432

3533
[massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(modelAn);
3634

3735
tabImbalance = [tab,table(imBalancedRxnBool,imBalancedCharge,imBalancedMass)];
3836
tabImbalance(~imBalancedRxnBool,:) = [];
37+
% Filter out exchange and SLIME reactions, these are known to be unbalanced
3938
tabImbalance(contains(tabImbalance.Var2,'exchange'),:) = [];
4039
tabImbalance(contains(tabImbalance.Var2,'SLIME'),:) = [];
41-
%tab=table(model.rxns,model.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,model.grRules);
42-
43-
44-
45-
out
46-
4740

4841

4942
%% Calculate formula and degree of reduciton of biomass
50-
[mwRange,metFormulae,elements,metEle]=computeMetFormulae(modelAn,'metMwRange','s_0450','fillMets','none','printLevel',0);
51-
Biomass_index = find(strcmp(modelAn.metNames,'biomass'));
52-
%Degree of reduction per element. Order of the elements 'C', 'H', 'O', 'N'
53-
DR_per_ele = [4, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
54-
DR_per_Cmol=sum(metEle(Biomass_index,:).*DR_per_ele)/metEle(Biomass_index,1)
55-
Biomass_formula_Cmol=metEle(Biomass_index,:)/metEle(Biomass_index,1);
56-
MW_per_Cmol_min = mwRange/metEle(Biomass_index,1)
57-
58-
43+
% [mwRange,metFormulae,elements,metEle]=computeMetFormulae(modelAn,'metMwRange','s_0450','fillMets','none','printLevel',0);
44+
% Biomass_index = find(strcmp(modelAn.metNames,'biomass'));
45+
% %Degree of reduction per element. Order of the elements 'C', 'H', 'O', 'N'
46+
% DR_per_ele = [4, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
47+
% DR_per_Cmol=sum(metEle(Biomass_index,:).*DR_per_ele)/metEle(Biomass_index,1)
48+
% Biomass_formula_Cmol=metEle(Biomass_index,:)/metEle(Biomass_index,1);
49+
% MW_per_Cmol_min = mwRange/metEle(Biomass_index,1)
50+
%
51+
%
5952

‎code/modelTests/growth.m‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -82,15 +82,19 @@
8282
%Simulate chemostats:
8383
mod_data = zeros(size(exp_data));
8484
solresult = zeros(length(model_origin.rxns),length(exp_data(:,1)));
85+
8586
if mode1 == 2
8687
model_origin = anaerobicModel(model_origin);
8788
end
89+
8890
if strcmp(mode2,'N')
89-
model_origin = scaleBioMass(model_origin,'protein',0.289,'',false);
90-
model_origin = scaleBioMass(model_origin,'lipid',0.048,'',false);
91-
model_origin = scaleBioMass(model_origin,'RNA',0.077,'carbohydrate',false);
91+
% P content in NH3-lim 0.1/h chemostat, per 10.1016/j.femsyr.2005.04.003
92+
model_origin = scaleBioMass(model_origin,'protein',0.28,'',false);
93+
% Assume that RNA decreased by the same amount (40%)
94+
model_origin = scaleBioMass(model_origin,'RNA',0.0329,'carbohydrate',false);
9295
model_origin = setParam(model_origin,'ub','r_0472',1000); %Glutamate synthase repressed in excess nitrogen
9396
end
97+
9498
for i = 1:length(exp_data(:,1))
9599
model_test= model_origin;
96100
%Fix glucose uptake rate and maximize growth:

0 commit comments

Comments
 (0)