fix: unbalanced reactions and incorrect metFormulas (SysBioChalmers#360)

* feat: increaseVersion uses YAML model

* fix: correct metFormulas of dolichol-related mets

* feat: additional correcctions

* fix: r_4196 NAD missing

* chore: generate new model files

README.md

@@ -33,7 +33,7 @@ This repository contains the current consensus genome-scale metabolic model of _
 
 | Taxonomy | Latest update | Version | Reactions | Metabolites | Genes |
 |:-------|:--------------|:------|:------|:----------|:-----|
-| _Saccharomyces cerevisiae_ | 16-Oct-2023 | develop | 4131 | 2806 | 1163 |
+| _Saccharomyces cerevisiae_ | 02-Dec-2023 | develop | 4130 | 2805 | 1162 |
 
 # Installation & usage
 

code/increaseVersion.m

@@ -47,7 +47,7 @@ function increaseVersion(bumpType)
 
 %Load model:
 disp('Loading model file')
-model = importModel('../model/yeast-GEM.xml');
+model = readYAMLmodel('../model/yeast-GEM.yml');
 
 %Run tests
 cd modelTests
@@ -109,6 +109,7 @@ function increaseVersion(bumpType)
 %Include tag and save model:
 disp('Write model files')
 model.id = ['yeastGEM_v' newVersion];
+model.version = newVersion;
 saveYeastModel(model,true,true,true)   %only save if model can grow
 
 %Check for any unexpected file changes
@@ -127,9 +128,9 @@ function increaseVersion(bumpType)
                     change = true;
                 end
             case 'model/yeast-GEM.yml'
-                %.yml file: 2 lines should be added & 2 lines should be deleted
-                %(1 with version information, 1 with current date)
-                if eval([diff_i{1} ' > 2']) || eval([diff_i{2} ' > 2'])
+                %.yml file: 3 lines should be added & 3 lines should be deleted
+                %(2 with version information, 1 with current date)
+                if eval([diff_i{1} ' > 3']) || eval([diff_i{2} ' > 3'])
                     disp(['NOTE: File ' diff_i{3} ' is changing more than expected'])
                     change = true;
                 end                

code/modelCuration/v9_0_0.m

@@ -0,0 +1,67 @@
+% This scripts applies curations to be applied on yeast-GEM release 8.7.1, to
+% get to yeast-GEM release 9.0.0.
+% Indicate which Issue/PR are addressed. If multiple curations are performed
+% before a new release is made, just add the required code to this script. If
+% more extensive coding is required, you can write a separate (generic) function
+% that can be kept in the /code/modelCuration folder. Otherwise, try to use
+% existing functions whenever possible. In particular /code/curateMetsRxnsGenes
+% can do many types of curation.
+
+%% Load yeast-GEM 8.7.1 (requires local yeast-GEM git repository)
+cd ..
+codeDir=pwd();
+model = getEarlierModelVersion('8.7.1');
+model.id='yeastGEM_develop';
+model.version='';
+%dataDir=fullfile(pwd(),'..','data','modelCuration','v9.0.0'); % No dataDir required for these curations
+cd modelCuration
+
+%% Correct inbalanced reactions, based on metFormulas
+% While dolichol can have any number of isoprenoid units, in yeast-GEM it is
+% defined as 4 units. This means that there is no need to keep R-subgroups as
+% part of dolichol-derived metabolites to indicate that unspecified length.
+% Define which metabolites are dolichol-derived and remove the R from their
+% metabolite formula.
+dolMets = getIndexes(model,{'s_3765','s_3767','s_3888','s_3911'},'mets');
+model.metFormulas(dolMets) = regexprep(model.metFormulas(dolMets),'R','');
+
+% r_4722 (polyphosphate polymerase) is unbalanced, 2 ADP is missing as product
+model = changeRxns(model,'r_4722','2 ATP[c] + H2O[c] => H+[c] + polyphosphate[v] + 2 ADP[c]',3);
+
+% r_4240 is unbalanced, but also has a generic reactant (protein asparagine)
+% and describes a process (protein modification) that is out-of-scope of a
+% metabolic model. Remove the reactions to resolve all 3 issues at once.
+model = removeReactions(model,'r_4240',true,true,true);
+
+% Some glycan metFormulas were incorrect absent. They were manually curated
+% by drawing out the structures.
+glycMets = getIndexes(model,{'s_3932','s_4003','s_4002'},'mets');
+model.metFormulas(glycMets) = {'C50H82N4O37R','C38H70N2O36P2R2','C32H60N2O31P2R2'};
+
+% r_0774 and r_0775 are unbalanced due to a missing H2O.
+model = changeRxns(model,{'r_0774','r_0775'},...
+    {'ATP[c] + H+[c] + nicotinate[c] + PRPP[c] => ADP[c] + diphosphate[c] + nicotinic acid D-ribonucleotide[c] + phosphate[c]',...
+     'ATP[m] + H+[m] + nicotinate[m] + PRPP[m] => ADP[m] + diphosphate[m] + nicotinic acid D-ribonucleotide[m] + phosphate[m]'},3);
+
+% r_4196 (NADH:ferricytochrome-b5 oxidoreductase) is unbalanced, NAD is missing as product
+model = changeRxns(model,'r_4196','NADH[erm] + 2 Ferricytochrome b5[erm] <=> H+[erm] + 2 Ferrocytochrome b5[erm] + NAD[erm]',3);
+
+%% DO NOT CHANGE OR REMOVE THE CODE BELOW THIS LINE.
+% Show some metrics:
+cd(fullfile(codeDir,'modelTests'))
+disp('Run gene essentiality analysis')
+[new.accuracy,new.tp,new.tn,new.fn,new.fp] = essentialGenes(model);
+fprintf('Genes in model: %d\n',numel(model.genes));
+fprintf('Gene essentiality accuracy: %.4f\n', new.accuracy);
+fprintf('True non-essential genes: %d\n', numel(new.tp));
+fprintf('True essential genes: %d\n', numel(new.tn));
+fprintf('False non-essential genes: %d\n', numel(new.fp));
+fprintf('False essential genes: %d\n', numel(new.fn));
+fprintf('\nRun growth analysis\n')
+R2=growth(model);
+fprintf('R2 of growth prediction: %.4f\n', R2);
+
+% Save model:
+cd ..
+saveYeastModel(model)
+cd modelCuration

data/databases/model_metDeltaG.csv

@@ -2372,7 +2372,6 @@ s_3882,10000000
 s_3883,-356.59
 s_3884,10000000
 s_3885,10000000
-s_3886,10000000
 s_3887,10000000
 s_3888,83.01
 s_3889,10000000

data/databases/model_rxnDeltaG.csv

@@ -3602,7 +3602,6 @@ r_4236,16.42
 r_4237,10000000
 r_4238,10000000
 r_4239,10000000
-r_4240,10000000
 r_4241,10000000
 r_4242,10000000
 r_4243,10000000

model/dependencies.txt

@@ -1,4 +1,4 @@
 MATLAB	9.12.0.2039608 (R2022a) Update 5
 libSBML	5.19.0
-RAVEN_toolbox	commit 2028181
+RAVEN_toolbox	2.8.6
 COBRA_toolbox	commit 28acb94

model/yeast-GEM.txt

@@ -3561,7 +3561,7 @@ r_4192	S-adenosyl-L-methionine[c] + [cytochrome c]-L-lysine[c] <=> H+[c] + S-ade
 r_4193	S-adenosyl-L-methionine[n] + L-lysine-[histone][n] <=> H+[n] + S-adenosyl-L-homocysteine[n] + N6-methyl-L-lysine-[histone][n] 	YHR119W or YDR440W or YJL168C	-1000.00	1000.00	  0.00
 r_4194	Zn(2+)[c] <=> Zn(2+)[er] 	YIL023C	-1000.00	1000.00	  0.00
 r_4195	Zn(2+)[m] <=> Zn(2+)[v] 	YMR243C or YOR316C	-1000.00	1000.00	  0.00
-r_4196	NADH[erm] + 2 Ferricytochrome b5[erm] <=> H+[erm] + 2 Ferrocytochrome b5[erm] 	YIL043C or YML125C	-1000.00	1000.00	  0.00
+r_4196	NADH[erm] + 2 Ferricytochrome b5[erm] <=> H+[erm] + NAD[erm] + 2 Ferrocytochrome b5[erm] 	YIL043C or YML125C	-1000.00	1000.00	  0.00
 r_4197	NADH[m] + 2 Ferricytochrome b5[m] <=> H+[m] + NAD[m] + 2 Ferrocytochrome b5[m] 	YIL043C or YKL150W or YML125C	-1000.00	1000.00	  0.00
 r_4198	NADP(+)[c] + (S)-benzoin[c] <=> H+[c] + NADPH[c] + benzil[c] 	YIR036C	-1000.00	1000.00	  0.00
 r_4199	RX[er] + glutathione[er] <=> HX[er] + R-S-glutathione[er] 	YIR038C	-1000.00	1000.00	  0.00
@@ -3602,7 +3602,6 @@ r_4236	(R)-lactate[c] + H+[c] <=> H2O[c] + methylglyoxal[c] 	YDR533C or YMR322C
 r_4237	H2O[v] + Ca(2+)[c] + ATP[v] => H+[v] + phosphate[v] + Ca(2+)[v] + ADP[v] 	YGL006W	  0.00	1000.00	  0.00
 r_4238	H2O[g] + Ca(2+)[c] + ATP[g] => H+[g] + phosphate[g] + Ca(2+)[g] + ADP[g] 	YGL167C	  0.00	1000.00	  0.00
 r_4239	Arg-tRNA(Arg)[c] + generic protein[c] <=> H+[c] + tRNA(Arg)[c] + L-Arginyl-protein[c] 	YGL017W	-1000.00	1000.00	  0.00
-r_4240	Protein asparagine[er] + G00008[er] <=> H+[er] + Dolichyl diphosphate[er] + G00009[er] 	YGL022W	-1000.00	1000.00	  0.00
 r_4241	H2O[er] + G00009[er] <=> G00171[er] + D-glucose[er] 	YGL027C	-1000.00	1000.00	  0.00
 r_4242	GDP-alpha-D-mannose[er] + G10694[er] <=> H+[er] + GDP[er] + G01813[er] 	YGL038C	-1000.00	1000.00	  0.00
 r_4243	GDP-alpha-D-mannose[er] + G00003[er] <=> H+[er] + GDP[er] + G00004[er] 	YGL065C	-1000.00	1000.00	  0.00
@@ -4073,7 +4072,7 @@ r_4718	H2O[c] + isobutyraldehyde[c] + NAD[c] <=> 2 H+[c] + NADH[c] + isobutyrate
 r_4719	2-methylbutanal[c] + H2O[c] + NAD[c] <=> 2 H+[c] + NADH[c] + 2-methylbutyrate[c] 	YOR374W or YPL061W	-1000.00	1000.00	  0.00
 r_4720	ethanol[c] + H+[c] + isobutyrate[c] <=> H2O[c] + ethyl isobutyrate[c] 	YOR126C	-1000.00	1000.00	  0.00
 r_4721	ethanol[c] + H+[c] + 2-methylbutyrate[c] <=> H2O[c] + ethyl 2-methylbutyrate[c] 	YOR126C	-1000.00	1000.00	  0.00
-r_4722	2 ATP[c] + H2O[c] => H+[c] + polyphosphate[v] 	YJL012C	  0.00	1000.00	  0.00
+r_4722	2 ATP[c] + H2O[c] => 2 ADP[c] + H+[c] + polyphosphate[v] 	YJL012C	  0.00	1000.00	  0.00
 r_4723	H2O[v] + polyphosphate[v] <=> H+[v] + 2 phosphate[v] 	YHR201C or YDR452W	-1000.00	1000.00	  0.00
 r_4724	H2O[m] + polyphosphate[m] <=> H+[m] + 2 phosphate[m] 	YHR201C or YDR452W	-1000.00	1000.00	  0.00
 r_4725	H2O[n] + polyphosphate[n] <=> H+[n] + 2 phosphate[n] 	YHR201C or YDR452W	-1000.00	1000.00	  0.00