feat: include anaerobic curation figures

Author: edkerk

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_ | 26-Feb-2026 | develop | 4102 | 2748 | 1143 |
+| _Saccharomyces cerevisiae_ | 26-May-2026 | develop | 4102 | 2748 | 1143 |
 
 # Installation & usage
 

code/modelTests/anaerobic_flux_predictions.m

@@ -9,7 +9,6 @@
 
 colors = orderedcolors("glow12");
 
-figure;
 data_sets=unique(text_flux(:,6));
 merged_data=[];
 merged_sim=[];
@@ -55,8 +54,8 @@
 plot(x,y,'--','MarkerSize',6,'Color',[64,64,64]/256)
 ylim([0 threshold])
 xlim([0 threshold])
-text(5,threshold/2,['mean relative error: ' num2str(mean_relative_error)]);
-text(5,threshold/2-10,['R^2: ' num2str(R2)]);
+text(12,threshold/2-10,['mean relative error: ' num2str(mean_relative_error)]);
+text(12,threshold/2-5,['R^2: ' num2str(R2)]);
 
 legend(data_sets);
 xlabel('Experimental 100 \cdot v_i/v_{Glx}','FontSize',14,'FontName','Helvetica')

code/modelTests/anaerobiosis.m

@@ -1,51 +1,83 @@
 clear; close all
-% Load the model and apply the corrections for *all* models
+% Load the old and new model
 model902 = getEarlierModelVersion('9.0.2');
-cd('../modelCuration/')
-v9_1_0;
+model910 = loadYeastModel();
 
 %% Run growth tests
-R2new = growth(model);
-R2old = growthOld(model);
-% R increased from 0.8528 to 0.9085
+funcs   = {@growth, @growthOld, @growth, @growthOld};
+models  = {model910, model910, model902, model902};
+titles  = {'v9.1.0 model, new anaerobic script', ...
+           'v9.1.0 model, old anaerobic script', ...
+           'v9.0.2 model, new anaerobic script', ...
+           'v9.0.2 model, old anaerobic script'};
 
-% Repeat with the model *before* generic 9.1.0 curations were done
-R2new902 = growth(model902);
-R2old902 = growthOld(model902);
-% R decreased from 0.8256 to 0.9001
+fig1 = figure('Name', 'Growth Comparison', 'Position', [100 100 1000 800]);
+for k = 1:4
+    subplot(2,2,k);
+    funcs{k}(models{k});
+    title(titles{k});
+end
+sgtitle('Growth: v9.1.0 vs v9.0.2, new vs old anaerobic script');
+saveas(fig1, '..\..\data\testResults\v910_growth.png');
+
+% R2 of growth prediction increased from 0.8256 to 0.9085 as result of all
+% curations.
 
 %% Convert to anaerobic
 cd('../otherChanges/')
-modelAn = anaerobicModel(model);
-modelAnOld = anaerobicModelOld(model);
+modelAn910 = anaerobicModel(model910);
+modelAnOld910 = anaerobicModelOld(model910);
 modelAn902 = anaerobicModel(model902);
 modelAnOld902 = anaerobicModelOld(model902);
 
+%% Anaerobic flux predictions
 cd('../modelTests/');
-%% flux predictions
-R2fluxnew = anaerobic_flux_predictions(modelAn);
-R2fluxold = anaerobic_flux_predictions(modelAnOld);
-R2fluxnew902 = anaerobic_flux_predictions(modelAn902);
-R2fluxold902 = anaerobic_flux_predictions(modelAnOld902);
-% R 
+models = {modelAn910, modelAnOld910, modelAn902, modelAnOld902};
+titles  = {'v9.1.0 model, new anaerobic script', ...
+           'v9.1.0 model, old anaerobic script', ...
+           'v9.0.2 model, new anaerobic script', ...
+           'v9.0.2 model, old anaerobic script'};
+fig2 = figure('Name', 'Anaerobic Comparison', 'Position', [100 100 1000 800]);
+for k = 1:4
+    subplot(2,2,k);
+    anaerobic_flux_predictions(models{k});
+    title(titles{k});
+end
+sgtitle('Anaerobic flux predictions: v9.1.0 vs v9.0.2, new vs old anaerobic script');
+saveas(fig2, '..\..\data\testResults\v910_anaerobic_fluxes.png');
 
+% R2 of flux predictions increased from 0.7858 to 0.9175 as a response to
+% all curations
 
 %% Plot
-plotAnaerobic(modelAn)
+models = {modelAn910, modelAnOld910, modelAn902, modelAnOld902};
+titles  = {'v9.1.0 model, new anaerobic script', ...
+           'v9.1.0 model, old anaerobic script', ...
+           'v9.0.2 model, new anaerobic script', ...
+           'v9.0.2 model, old anaerobic script'};
+fig3 = figure('Name', 'Anaerobic Comparison', 'Position', [100 100 1000 800]);
+for k = 1:4
+    subplot(2,2,k);
+    plotAnaerobic(models{k});
+    title(titles{k});
+end
+sgtitle('Anaerobic fluxes (Sjöberg data): v9.1.0 vs v9.0.2, new vs old anaerobic script');
+saveas(fig3, '..\..\data\testResults\v910_anaerobic_sjoberg.png');
+
 %% Set glucose uptake rate and solve pFBA
-modelAn = setParam(modelAn,'eq','r_1714',-23);
-res=solveLP(modelAn,1);
-FLUX = res.x;
-v_AStr = res.x(getIndexes(modelAn,'r_1115','rxns'));
-v_ATPase = res.x(getIndexes(modelAn,'r_0227','rxns'));
-v_glc = res.x(getIndexes(modelAn,'r_1714','rxns'));
+temp_model = setParam(modelAn910,'eq','r_1714',-23);
+res=solveLP(temp_model,1); FLUX = res.x;
+
+v_AStr = res.x(getIndexes(temp_model,'r_1115','rxns')); % Ammonium exchange
+v_ATPase = res.x(getIndexes(temp_model,'r_0227','rxns')); % ATPase
+v_glc = res.x(getIndexes(temp_model,'r_1714','rxns')); % Glucose uptake
+[v_AStr, v_ATPase, v_glc]
 
 %% Pack flux results into table
-temp_model=modelAn;
-rxns_reacs=constructEquations(temp_model,modelAn.rxns);
-tab=table(modelAn.rxns,modelAn.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,modelAn.grRules);
+rxns_reacs=constructEquations(temp_model,temp_model.rxns);
+tab=table(temp_model.rxns,temp_model.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,temp_model.grRules);
 
-[massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(modelAn);
+[massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(temp_model);
 
 tabImbalance = [tab,table(imBalancedRxnBool,imBalancedCharge,imBalancedMass)];
 tabImbalance(~imBalancedRxnBool,:) = [];

code/modelTests/anaerobiosis.mlx

@@ -36,7 +36,6 @@
 mod_data(27:32,:) = simulateChemostat(model_origin,exp_data(27:32,:),2,'N');
 
 % plot the figure
-figure
 hold on
 cols = [215,25,28;253,174,97;171,217,233;44,123,182]/256;
 b(1) = plot(exp_data(1:9,4),mod_data(1:9,4),'o','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(2,:));
@@ -60,6 +59,7 @@
 hold off
 R2=corrcoef(exp_data(:,4),mod_data(:,4));
 R2=R2(2)^2;
+text(0.25,0.1,['R2:',num2str(R2)])
 
 if writeOutput
     saveas(gcf,fullfile(funcDir,'..','..','data','testResults','growth.png'));

code/modelTests/growth.m

@@ -36,7 +36,6 @@
 mod_data(27:32,:) = simulateChemostat(model_origin,exp_data(27:32,:),2,'N');
 
 % plot the figure
-figure
 hold on
 cols = [215,25,28;253,174,97;171,217,233;44,123,182]/256;
 b(1) = plot(exp_data(1:9,4),mod_data(1:9,4),'o','MarkerSize',10,'MarkerEdgeColor','k','MarkerFaceColor',cols(2,:));
@@ -60,6 +59,7 @@
 hold off
 R2=corrcoef(exp_data(:,4),mod_data(:,4));
 R2=R2(2)^2;
+text(0.25,0.1,['R2:',num2str(R2)])
 
 if writeOutput
     saveas(gcf,fullfile(funcDir,'..','..','data','testResults','growth.png'));

code/modelTests/growthOld.m

@@ -12,7 +12,7 @@ function plotAnaerobic(modelAn)
 v_AStr = FLUX(getIndexes(modelAn,'r_1115','rxns'));
 v_ATPase = FLUX(getIndexes(modelAn,'r_0227','rxns'));
 %% Show relative accuracy of main extracellular products
-figure;
+%figure;
 %glycerol ethanol Co2
 %4.5 ± 0.4  31 ± 2  38 ± 10
 data=[4.5 31 38 0.36];