tutorial3_solutions

 tutorial3_solutions
   This script contains the solutions for Tutorial 3, see Tutorial 3 in
   "RAVEN tutorials.docx" for more details. All the parameters are set in
   this script, rather than modifying the Excel model file.

0001 % tutorial3_solutions
0002 %   This script contains the solutions for Tutorial 3, see Tutorial 3 in
0003 %   "RAVEN tutorials.docx" for more details. All the parameters are set in
0004 %   this script, rather than modifying the Excel model file.
0005 
0006 %Import the Excel model
0007 model=importExcelModel('smallYeast.xlsx',true);
0008 
0009 %Step 1
0010 %Set the upper bound of glucose uptake to 1 and O2 uptake to zero
0011 model=setParam(model,'ub',{'glcIN' 'o2IN'},[1 0]);
0012 
0013 %Set the objective to be ATP hydrolysis
0014 model=setParam(model,'obj',{'ATPX'},1);
0015 
0016 %Solve the model
0017 sol=solveLP(model);
0018 
0019 %Print the resulting fluxes. The ATP production rate should be 2.0. It was
0020 %4.0 in Tutorial 2, but there sucrose was used instead of glucose.
0021 printFluxes(model,sol.x,false);
0022 
0023 %Step 2
0024 %Check the yield of different products and print the results
0025 %Change to fully aerobic
0026 model=setParam(model,'ub',{'glcIN' 'o2IN'},[1 1000]);
0027 model=setParam(model,'obj',{'ethOUT'},1);
0028 sol=solveLP(model);
0029 fprintf(['Yield of ethanol is ' num2str(sol.f) ' mol/mol\n']);
0030 model=setParam(model,'obj',{'acOUT'},1);
0031 sol=solveLP(model);
0032 fprintf(['Yield of acetate is ' num2str(sol.f) ' mol/mol\n']);
0033 model=setParam(model,'obj',{'glyOUT'},1);
0034 sol=solveLP(model);
0035 fprintf(['Yield of glycerol is ' num2str(sol.f) ' mol/mol\n']);
0036 model=setParam(model,'obj',{'biomassOUT'},1);
0037 sol=solveLP(model);
0038 fprintf(['Yield of biomass is ' num2str(sol.f) '/h\n']);
0039 
0040 %Step 3
0041 %Solve for both aerobic and anaerobic growth
0042 solA=solveLP(model);
0043 model=setParam(model,'ub',{'o2IN'},0.5);
0044 solB=solveLP(model);
0045 
0046 %Plot the differences
0047 %Load the map
0048 load 'pathway.mat' pathway;
0049 drawMap('Aerobic vs Anaerobic',pathway,model,solA.x,solB.x,[],'mapFBA.pdf',10^-5);
0050 
0051 %Step 4
0052 %Change to anaerobic growth and maximize for biomass
0053 model=setParam(model,'eq',{'o2IN'},0);
0054 model=setParam(model,'obj',{'biomassOUT'},1);
0055 sol=solveLP(model);
0056 printFluxes(model,sol.x,true);
0057 %One can see that the model predicts a glycerol production of 0.23
0058 %mmol/gDW/h
0059 
0060 %Run a single gene deletion
0061 [genes, fluxes, originalGenes, details]=findGeneDeletions(model,'sgd','fba');
0062 
0063 %Get the indexes of these reactions
0064 I=getIndexes(model,{'biomassOUT'},'rxns');
0065 J=getIndexes(model,{'glyOUT'},'rxns');
0066 
0067 okSolutions=find(fluxes(I,:)>10^-2); %Only look at solutions which are still growing
0068 [maxGlycerol, J]=max(fluxes(J,okSolutions));
0069 fprintf(['Glycerol production is ' num2str(maxGlycerol) ' after deletion of ' originalGenes{genes(okSolutions(J),:)} '\n']);
0070 
0071 %The best gene deletion corresponds to turning off the ZWF1 reaction
0072 %(YNL241C)
0073 model2=setParam(model,'eq',{'ZWF'},0);
0074 sol2=solveLP(model2);
0075 drawMap('ZWF1 deletion vs WT',pathway,model,sol2.x,sol.x,[],'mapZWF.pdf',10^-5);
0076 followChanged(model,sol2.x,sol.x, 10, 10^-2, 0,{'NADPH' 'NADH' 'NAD' 'NADP'});
0077 
0078 %Step 5
0079 %Set the exchange rates to the recorded batch values
0080 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);
0081 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);
0082 
0083 %Define another model where all exchange reactions are open.
0084 model2=model;
0085 I=getIndexes(model,getExchangeRxns(model),'rxns');
0086 model2.lb(I)=0;
0087 model2.ub(I)=1000;
0088 
0089 %Delete ZWF gene
0090 model2=setParam(model2,'eq',{'ZWF'},0);
0091 
0092 %Run MOMA
0093 [fluxA, fluxB, flag]=qMOMA(model,model2);
0094 drawMap('ZWF deletion vs wild type',pathway,model,fluxB,fluxA,[],'mapMOMA.pdf',10^-5);
0095 
0096 %As one can see, the glycerol production is higher in the deletion strain.
0097 %Note that this is without any objectives, just by trying to maintain the
0098 %cells original flux distribution.
0099 
0100 %Step 6
0101 %Read microarray results and calculate reporter metabolites (metabolites
0102 %around which there are significant transcriptional changes)
0103 [orfs, pvalues]=textread('expression.txt','%s%f');
0104 repMets=reporterMetabolites(model,orfs,pvalues);
0105 [I, J]=sort(repMets.metPValues);
0106 
0107 fprintf('TOP 10 REPORTER METABOLITES:\n');
0108 for i=1:min(numel(J),10)
0109     fprintf([repMets.mets{J(i)} '\t' num2str(I(i)) '\n']);
0110 end
0111 
0112 %Get all reactions involving those metabolites and display them on a map
0113 mets=ismember(model.mets,repMets.mets(J(1:10)));
0114 [~, I]=find(model.S(mets,:));
0115 pathway=trimPathway(pathway, model.rxns(I), true);
0116 drawMap('Reactions involving the top 10 Reporter Metabolites',pathway,model,ones(numel(model.rxns),1),zeros(numel(model.rxns),1),[],'mapRM.pdf',10^-5);