fix: curate auxotrophy and correct

README.md

@@ -39,7 +39,7 @@ This repository contains the current consensus genome-scale metabolic model of _
 
 | Taxonomy | Latest update | Version | Reactions | Metabolites | Genes |
 |:-------|:--------------|:------|:------|:----------|:-----|
-| _Saccharomyces cerevisiae_ | 27-Oct-2024 | develop | 4130 | 2805 | 1162 |
+| _Saccharomyces cerevisiae_ | 14-Nov-2024 | develop | 4131 | 2806 | 1161 |
 
 # Installation & usage
 

code/modelCuration/v9_0_2.m

@@ -0,0 +1,161 @@
+% This scripts applies curations to be applied on yeast-GEM release 9.0.1, to
+% get to yeast-GEM release 9.0.2.
+
+%% Load yeast-GEM 9.0.1 (requires local yeast-GEM git repository)
+cd ..
+codeDir=pwd();
+model = getEarlierModelVersion('9.0.1');
+model.id='yeastGEM_develop';
+model.version='';
+%dataDir=fullfile(pwd(),'..','data','modelCuration','v9.0.0'); % No dataDir required for these curations
+cd modelCuration
+
+%% Issue #371. Implement auxotrophy-based curation, as proposed by Han
+% et al., 2024 doi:10.1016/j.synbio.2024.07.006. Not all changes are
+% implemented as proposed by the authors.
+
+% To resolve AAT2 causing asparate auxotrophy
+model                       = setParam(model,'eq','r_0217',0);
+rxnIdx                      = getIndexes(model,'r_0217','rxns');
+model.rxnNotes{rxnIdx}      = 'No significant mitochondrial flux';
+model.rxnReferences{rxnIdx} = '10.1016/j.femsyr.2004.09.008';
+
+% To resolve ALD2 and ALD3 causing pantothenic acid auxotrophy
+model                       = changeGrRules(model,'r_0172','YMR170C or YMR169C',true);
+rxnIdx                      = getIndexes(model,'r_0172','rxns');
+model.rxnReferences{rxnIdx} = '10.1093/genetics/163.1.69';
+
+% To resolve CHO2 causing choline auxotrophy
+CHO2rxns                    = {'r_2488','r_2489','r_2490','r_2491','r_2492','r_2493','r_2494','r_2495'};
+model                       = changeGrRules(model,CHO2rxns,'YGR157W',true);
+rxnIdx                      = getIndexes(model,CHO2rxns,'rxns');
+model.rxnReferences(rxnIdx) = {'10.1016/0005-2760(90)90145-N'};
+% While the curation is correct, CHO2 and OPI3 are not strictly essential
+% genes, choline auxotrophy is only in relation to the production of
+% phosphocholines. If no cholines are available, the phospholipid
+% composition would be different. But because of the fixed biomass
+% composition in yeast-GEM, CHO2 and OPI3 act as essential genes. This
+% causes false positives in the essentialGenes function, as the genes not
+% truly essential.
+
+% To resolve CYS3 causing cysteine auxotrophy
+model                       = setParam(model,'eq','r_0312',0);
+rxnIdx                      = getIndexes(model,'r_0312','rxns');
+model.rxnNotes{rxnIdx}      = 'Regulatory role only';
+model.rxnReferences{rxnIdx} = '10.1111/j.1574-6968.2003.tb11531.x';
+
+% Han et al. suggest to block TUM1 reaction r_4703, to resolve CYS4 as
+% causing cysteine auxotrophy. However, there is good evidence that this
+% reaction is taking place (doi:10.1093/femsyr/fow100). Instead, make the
+% reaction irreversible.
+model                       = setParam(model,'lb','r_4703',0);
+
+% To resolve ERG10 causing ergosterol auxotrophy
+model = removeReactions(model,'r_0559',true,true); % ERG13, mitochondrial
+% only, see e.g. https://www.uniprot.org/uniprotkb/P54839/entry
+
+% To resolve GFA1 causing D-glucosamine auxotrophy: add chitin to
+% carbohydrate pseudoreaction and adjust other reactions proportionally
+rxn  = getIndexes(model,'r_4048','rxns');
+met1 = getIndexes(model,'s_0509','mets');%chitin
+met2 = getIndexes(model,'s_0001','mets');%(1->3)-beta-D-glucan
+met3 = getIndexes(model,'s_0004','mets');%(1->6)-beta-D-glucan
+met4 = getIndexes(model,'s_1107','mets');%mannan
+met5 = getIndexes(model,'s_0773','mets');%glycogen
+met6 = getIndexes(model,'s_1520','mets');%trehalose
+model.S(met1, rxn) = -0.02361;
+model.S(met2, rxn) = -0.73914;
+model.S(met3, rxn) = -0.24696;
+model.S(met4, rxn) = -0.70204;
+model.S(met5, rxn) = -0.35689;
+model.S(met6, rxn) = -0.13655;
+% And change gene association
+model = changeGrRules(model,'r_0477','YKL104C',true);
+
+% To resolve GSH1 causing glutatione auxotrophy: add glutathione to the
+% cofactor pseudoreaction
+rxn = getIndexes(model,'r_4598','rxns');
+met = getIndexes(model,'s_0750','mets');
+model.S(met, rxn) = -1e-06;
+% The above change will make GSH1 and GSH2 essential genes. This will be
+% reported in essentialGenes as false positives, as their essentiality is
+% debated (see for instance doi:10.1016/S1567-1356(02)00081-8). We here
+% argue that absence of glutathione very strongly affects cell growth and
+% stress tolerance.
+
+% To resolve ERG13/20 causing ergosterol auxotrophy; and HEM1/2/3/4/12 
+% causing heme auxotrophy. Heme transport and exchange should be present to
+% allow uptake if supplemented
+metsToAdd.metNames      = 'heme a';
+metsToAdd.compartments  = 'e';
+model                   = addMets(model,metsToAdd,true,'s_');
+model                   = addTransport(model,'c','e','heme a',true,false,'r_');
+model                   = addExchangeRxns(model,'out',model.mets(end));
+model.rxns(end)         = generateNewIds(model,'rxns','r_',1);
+model.rxnNames(end)     = {'heme a exchange'};
+model.subSystems{end}   = {'Exchange reaction'};
+
+% To resolve MET3 causing methionine auxotrophy
+rxnIdx = getIndexes(model,{'r_0223','r_0224','r_1026'},'rxns');
+model.eccodes(rxnIdx) = {'2.7.7.53','2.7.7.53','2.7.7.5'};
+% Han et al. suggest to block APA1 reaction r_1026, but their cited
+% reference (doi:10.1007/BF00264692) contains no reasoning for this.
+% Instead, there is good evidence for the function of this enzyme:
+% (doi:10.1128/jb.171.12.6437-6445.1989), although the reaction is assayed
+% in the reverse direction. So instead, invert the reaction (and keep it
+% irreversible), this has the same effect for MET3 as blocking it.
+model.S(:,rxnIdx(3)) = - model.S(:,rxnIdx(3));
+
+% To resolve MET13 causing methionine auxotrophy
+model = changeGrRules(model,'r_0080','(YGL125W and YPL023C) or YGL125W',true);
+
+% To resolve MET17 causing methionine auxotrophy
+model = setParam(model,'lb','r_0815',0);
+
+% To resolve THI4 causing thiamine auxotrophy and HOM2 causing methionine
+% and threonine auxotrophy
+model = removeReactions(model,'r_2070',true,true);
+model = removeReactions(model,'r_2071',true,true);
+rxnsToAdd.rxns      = generateNewIds(model,'rxns','r_',2);
+rxnsToAdd.equations = {'s_3910 + s_0803 => s_0423 + s_0293 + s_0456 + s_0794', ...
+                       's_1003 + s_1198 + s_0841 => s_1216 + s_3910 + 3 s_0803 + s_0794'};
+rxnsToAdd.rxnNames  = {'HET-P synthase (thiazole synthase)','adenylated thiazole synthase'};
+rxnsToAdd.lb        = [0,0];
+rxnsToAdd.ub        = [1000,1000];
+rxnsToAdd.eccodes   = {'','2.4.2.60'};
+rxnsToAdd.subSystems= {'Thiamine metabolism','Thiamine metabolism'};
+rxnsToAdd.grRules   = {'YGR144W','YGR144W'};
+rxnsToAdd.rxnMiriams{1} = struct('name',{{'kegg.reaction';'kegg.pathway';'metanetx.reaction'}},...
+                                'value',{{'R10685';'sce00730';'MNXR139808'}});
+rxnsToAdd.rxnMiriams{2} = struct('name',{{'kegg.reaction';'kegg.pathway';'metanetx.reaction'}},...
+                                'value',{{'R10711';'sce00730';'MNXR139812'}});
+rxnsToAdd.rxnReferences = {'10.1021/ja067606t','10.1021/ja067606t'};
+rxnsToAdd.rxnConfidenceScores = [3,3];
+model = addRxns(model,rxnsToAdd);
+
+% To resolve URA2 causing uracil auxotrophy
+model = changeGrRules(model,'r_0250','YJL130C or (YJR109C and YOR303W)',true);
+
+%% Issue #372. Correct UniProt identifier of YCR024C and YCR024C-A
+geneIdx = getIndexes(model,'YCR024C','genes');
+model.geneMiriams{geneIdx}.value = {'P25345'};
+
+%% 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,true);
+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,true);
+fprintf('R2 of growth prediction: %.4f\n', R2);
+
+% Save model:
+cd ..
+saveYeastModel(model)
+cd modelCuration

data/databases/model_metDeltaG.csv

@@ -2804,3 +2804,4 @@ s_4328,-593.58
 s_4329,-450.64
 s_4330,-659.91
 s_4331,-63.35
+s_4332,10000000

data/databases/model_rxnDeltaG.csv

@@ -442,7 +442,6 @@ r_0555,-27.78
 r_0556,10000000
 r_0557,15.51
 r_0558,-9.81
-r_0559,-8.4
 r_0560,-9.08
 r_0561,-17.09
 r_0562,-24.94
@@ -1467,8 +1466,6 @@ r_2066,10000000
 r_2067,10000000
 r_2068,10000000
 r_2069,10000000
-r_2070,10000000
-r_2071,10000000
 r_2072,0
 r_2073,10000000
 r_2074,-15.23
@@ -4129,3 +4126,7 @@ r_4775,NaN
 r_4776,NaN
 r_4777,NaN
 r_4778,NaN
+r_4779,0
+r_4780,NaN
+r_4781,NaN
+r_4782,NaN

data/testResults/essentialGenes.md

@@ -19,7 +19,6 @@ YCL052C
 YCR012W
 YDL004W
 YDL015C
-YDL103C
 YDL120W
 YDL166C
 YDR023W
@@ -34,9 +33,7 @@ YDR373W
 YDR376W
 YDR437W
 YDR487C
-YEL058W
 YER023W
-YFL017C
 YGL001C
 YGL008C
 YGL047W
@@ -62,15 +59,13 @@ YJL005W
 YJL026W
 YJL091C
 YJL097W
-YJL101C
 YJR013W
 YJR016C
 YJR057W
 YKL004W
 YKL019W
 YKL035W
 YKL060C
-YKL104C
 YKL141W
 YKL152C
 YKL165C
@@ -91,11 +86,9 @@ YMR301C
 YNL036W
 YNL038W
 YNL267W
-YOL049W
 YOL066C
 YOR143C
 YOR149C
-YPL028W
 YPL076W
 YPL175W
 YPL231W
@@ -106,6 +99,7 @@ YBR126C
 YDR074W
 YER014W
 YER026C
+YGR157W
 YJR073C
 YKL067W
 YKR039W
@@ -483,7 +477,6 @@ YGR143W
 YGR144W
 YGR154C
 YGR155W
-YGR157W
 YGR170W
 YGR171C
 YGR177C
@@ -1072,6 +1065,7 @@ YBL076C
 YBR029C
 YDL045C
 YDL055C
+YDL103C
 YDL141W
 YDL205C
 YDR037W
@@ -1082,8 +1076,10 @@ YDR232W
 YDR341C
 YDR454C
 YDR531W
+YEL058W
 YER003C
 YER043C
+YFL017C
 YFL022C
 YFL045C
 YGL040C
@@ -1098,9 +1094,11 @@ YHR020W
 YHR074W
 YIL078W
 YIL083C
+YJL101C
 YJL167W
 YKL024C
 YKL088W
+YKL104C
 YLL018C
 YLR060W
 YML126C
@@ -1111,6 +1109,7 @@ YNL247W
 YNL256W
 YNR016C
 YNR043W
+YOL049W
 YOL097C
 YOL130W
 YOR074C
@@ -1120,6 +1119,7 @@ YOR176W
 YOR236W
 YOR278W
 YOR335C
+YPL028W
 YPL117C
 YPL160W
 YPR033C

data/testResults/growth.md

@@ -1,4 +1,4 @@
 ## R2 of growth rate prediction
-0.8804
+0.8798
 
 ![Growth curve](growth.png)

model/dependencies.txt

@@ -1,4 +1,4 @@
 MATLAB	24.1.0.2568132 (R2024a) Update 1
 libSBML	5.20.2
-RAVEN_toolbox	2.10.1
+RAVEN_toolbox	commit c90390e
 COBRA_toolbox	3.4