Tutorial 3 — Knockouts, MOMA and omics data¶

This exercise shows how to run FBA and minimization of metabolic adjustment (MOMA) simulations, and how a GEM can serve as a scaffold for interpreting microarray data. It uses a simplified model of yeast metabolism (smallYeast.yml), imported with readYAMLmodel.

It is assumed you have completed Tutorial 2.

Skeleton and solutions

tutorial3.m is a skeleton script for you to work through. The completed answers are in the companion script tutorial3_solutions.m. The code excerpts below follow the solutions; try to derive them yourself before looking.

Step by step¶

1. Import and check ATP yield¶

Import the model, restrict it to anaerobic glucose growth, and maximize ATP hydrolysis (ATPX):

model = readYAMLmodel('smallYeast.yml');
model = setParam(model, 'ub', {'glcIN' 'o2IN'}, [1 0]);
model = setParam(model, 'obj', {'ATPX'}, 1);
sol = solveLP(model);
printFluxes(model, sol.x, false);

The ATP production rate should be 2.0. It was 4.0 in Tutorial 2, but there sucrose was used instead of glucose.

2. Check product yields¶

Switch to fully aerobic conditions and check the yield of different products, maximizing each in turn:

model = setParam(model, 'ub', {'glcIN' 'o2IN'}, [1 1000]);
model = setParam(model, 'obj', {'ethOUT'}, 1);
sol = solveLP(model);
fprintf(['Yield of ethanol is ' num2str(sol.f) ' mol/mol\n']);

The same is repeated for acetate (acOUT), glycerol (glyOUT) and biomass (biomassOUT).

3. Aerobic vs. limited oxygen¶

Solve for biomass production under full aeration, then restrict the oxygen uptake upper bound to 0.5 and solve again:

solA = solveLP(model);
model = setParam(model, 'ub', {'o2IN'}, 0.5);
solB = solveLP(model);

4. Single gene deletions with FBA¶

Set the model to anaerobic growth maximizing biomass, then scan single gene deletions. After running the deletion scan, keep only solutions that still grow and look for the one giving the highest glycerol production:

model = setParam(model, 'eq', {'o2IN'}, 0);
model = setParam(model, 'obj', {'biomassOUT'}, 1);
sol = solveLP(model);
printFluxes(model, sol.x, true);

[genes, fluxes, originalGenes, details] = findGeneDeletions(model, 'sgd', 'fba');

I = getIndexes(model, {'biomassOUT'}, 'rxns');
J = getIndexes(model, {'glyOUT'}, 'rxns');
okSolutions = find(fluxes(I,:) > 10^-2);          % still growing
[maxGlycerol, J] = max(fluxes(J, okSolutions));
fprintf(['Glycerol production is ' num2str(maxGlycerol) ...
    ' after deletion of ' originalGenes{genes(okSolutions(J),:)} '\n']);

The model predicts a glycerol production of 0.23 mmol/gDW/h for the unperturbed anaerobic case. The best gene deletion corresponds to turning off the ZWF reaction (ZWF1, YNL241C). Compare the deletion strain to wild type and follow the changed fluxes around the redox metabolites:

model2 = setParam(model, 'eq', {'ZWF'}, 0);
sol2 = solveLP(model2);
followChanged(model, sol2.x, sol.x, 10, 10^-2, 0, {'NADPH' 'NADH' 'NAD' 'NADP'});

5. MOMA¶

FBA assumes the cell re-optimises after a perturbation. MOMA instead assumes the perturbed cell changes its metabolism as little as possible — useful when you have wild-type data and want to predict a mutant. Constrain the wild-type model to the recorded batch exchange rates, then define an unconstrained model with ZWF knocked out:

model = setParam(model, 'ub', {'acOUT' 'biomassOUT' 'co2OUT' 'ethOUT' 'glyOUT' 'glcIN' 'o2IN' 'ethIN'}, [0 0.67706 22.4122 19.0946 1.4717 15 1.6 0]*1.0001);
model = setParam(model, 'lb', {'acOUT' 'biomassOUT' 'co2OUT' 'ethOUT' 'glyOUT' 'glcIN' 'o2IN' 'ethIN'}, [0 0.67706 22.4122 19.0946 1.4717 15 1.6 0]*0.9999);

model2 = model;
I = getIndexes(model, getExchangeRxns(model), 'rxns');
model2.lb(I) = 0;  model2.ub(I) = 1000;
model2 = setParam(model2, 'eq', {'ZWF'}, 0);

[fluxA, fluxB, flag] = qMOMA(model, model2);

The glycerol production is higher in the deletion strain. Note that this is without any objectives, just by trying to maintain the cell's original flux distribution.

6. Reporter metabolites from microarray data¶

A GEM can highlight the metabolites around which significant transcriptional changes cluster. Load expression data and run the reporter-metabolites test:

[orfs, pvalues] = textread('expression.txt', '%s%f');
repMets = reporterMetabolites(model, orfs, pvalues);
[I, J] = sort(repMets.metPValues);

fprintf('TOP 10 REPORTER METABOLITES:\n');
for i = 1:min(numel(J), 10)
    fprintf([repMets.mets{J(i)} '\t' num2str(I(i)) '\n']);
end

Full script¶

% tutorial3
%   This exercise shows how to run FBA and minimization of metabolic
%   adjustment (MOMA) simulations and how one can use GEMs as a scaffold
%   for interpreting microarray data. A simplified model of yeast
%   metabolism is used in this approach as an example.
%   See Tutorial 3 in "RAVEN tutorials.docx" for more details.
%
%   It is assumed that the user has already completed Tutorial 2

%Import the model
model=readYAMLmodel('smallYeast.yml');

%Set the upper bound of glucose uptake to 1 and O2 uptake to unlimited
model=setParam(model,'ub',{'glcIN' 'o2IN'},[1 1000]);

%Set the objective to be ethanol production
model=setParam(model,'obj',{'ethOUT'},1);

%Solve the model
sol=solveLP(model);

%Print the resulting exchange fluxes
printFluxes(model,sol.x,true);

%Run a single gene deletion
[genes, fluxes, originalGenes, details]=findGeneDeletions(model,'sgd','fba');

%Get the indexes of these reactions
I=getIndexes(model,{'biomassOUT'},'rxns');
J=getIndexes(model,{'glyOUT'},'rxns');

okSolutions=find(fluxes(I,:)>10^-2); %Only look at solutions which are still growing
[maxGlycerol, J]=max(fluxes(J,okSolutions));
disp(maxGlycerol);
disp(originalGenes(genes(okSolutions(J),:)));

%Compare the ZWF1 deletion strain to wild type
model2=setParam(model,'eq',{'ZWF'},0);
sol2=solveLP(model2);
followChanged(model,sol2.x,sol.x, 10, 10^-2, 0,{'NADPH' 'NADH' 'NAD' 'NADP'});

%Reload the model with exchange metabolites removed for simulation
model=readYAMLmodel('smallYeast.yml');
sol=solveLP(model);

%Define another model where all exchange reactions are open
model2=model;
I=getIndexes(model,getExchangeRxns(model),'rxns');
model2.lb(I)=0;
model2.ub(I)=1000;

%Delete ZWF gene
model2=setParam(model2,'eq',{'ZWF'},0);

%Run MOMA
[fluxA, fluxB, flag]=qMOMA(model,model2);

%Read microarray results and calculate reporter metabolites (metabolites
%around which there are significant transcriptional changes)
[orfs, pvalues]=textread('expression.txt','%s%f');
repMets=reporterMetabolites(model,orfs,pvalues);
[I, J]=sort(repMets.metPValues);

fprintf('TOP 10 REPORTER METABOLITES:\n');
for i=1:min(numel(J),10)
    fprintf([repMets.mets{J(i)} '\t' num2str(I(i)) '\n']);
end