diff --git a/analysis/test/README.md b/analysis/test/README.md new file mode 100644 index 00000000..f5d149cc --- /dev/null +++ b/analysis/test/README.md @@ -0,0 +1,7 @@ +

+ pdf   MLX   M file   viewtutorials +

+

+

+ +

diff --git a/analysis/test/tutorial_pFBA.m b/analysis/test/tutorial_pFBA.m new file mode 100644 index 00000000..5e93eb8f --- /dev/null +++ b/analysis/test/tutorial_pFBA.m @@ -0,0 +1,300 @@ +%% A step-by-step guide to parsimoneous enzyme usage Flux Balance Analysis - pFBA +% *Author(s): Francisco José Pardo Palacios, Ines Thiele, LCSB, University of +% Luxembourg.* +% +% *Reviewer(s): Sebastián Mendoza, Center for Mathematical Modeling, University +% of Chile. Catherine Clancy, LCSB, University of Luxembourg.* +% +% *Lin Wang (Costas D. Maranas Lab), Joshua Chan (Costas D. Maranas Lab), +% Chiam Yu Ng (Costas D. Maranas Lab)* +%% INTRODUCTION +% This tutorial shows how the parsimoneous enzyme usage Flux Balance Analysis +% (pFBA), as described in Lewis et al.$$^1$, has been implemented in The COBRA +% Toolbox as the function |pFBA()|. +% +% The main aim of the tutorial is to explain how the calculations are carried +% out in order to understand the pFBA analysis, and to be able to classify, under +% certain conditions, the genes of a model as: essential, pFBA optima, Enzymatically +% Less Efficient (ELE), Metabolically Less Efficient (MLE) or pFBA no-flux genes +% (Figure 1). +% +% # _*Essential genes:_* metabolic genes necessary for growth in the given media. +% # _*pFBA optima: _*non-essential genes contributing to the optimal growth +% rate and minimum gene-associated flux. +% # _*Enzymatically less efficient (ELE): _*genes requiring more flux through +% enzymatic steps than alternative pathways that meet the same predicted growth +% rate. +% # _*Metabolically less efficient (MLE):_* genes requiring a growth rate reduction +% if used. +% # _*pFBA no-flux:_* genes that are unable to carry flux in the experimental +% conditions. +% +% +% +% +% Figure 1: Gene/enzyme classification scheme used by pFBA +% +% This tutorial will use the_ E. coli core _reconstruction$$^2$ as the model +% of choice, and will be called herein as |modelEcore|. The results obtained could +% then be compared to data from evolved E. coli and observe if the |modelEcore| +% can predict its evolution. In order to investigate this, all the steps described +% in the pFBA flowchart (Figure 2) should be followed, and are demonstrated in +% this tutorial. +% +% +% +% +% Figure 2: pFBA flowchart +%% EQUIPMENT SETUP +%% *Initialize the COBRA Toolbox.* +% If necessary, initialize The Cobra Toolbox using the |initCobraToolbox| function. +%% +initCobraToolbox(false) % false, as we don't want to update +%% *Setting the *optimization* solver.* +% This tutorial will be run with a |'glpk'| package, which is a linear programming +% ('|LP'|) solver. The |'glpk'| package does not require additional instalation +% and configuration. + +solverName = 'glpk'; +solverType = 'LP'; +changeCobraSolver(solverName, solverType); +% changeCobraSolver('ibm_cplex'); +% changeCobraSolver('gurobi'); +%% +% However, for the analysis of larger models, such as Recon 2.04$$^3$, it +% is not recommended to use the |'glpk'| package but rather an industrial strength +% solver, such as the |'gurobi' or 'ibm_cplex'| package. +% +% A solver package may offer different types of optimization programmes to +% solve a problem. The above example used a LP optimization, other types of optimization +% programmes include; mixed-integer linear programming ('|MILP|'), quadratic programming +% ('|QP|'), and mixed-integer quadratic programming ('|MIQP|'). +%% *Model setup.* +% Load the modelEcore, and define the uptake of nutrients by the modelEcore. +% The substrate used is glucose, and for this tutorial limit its uptake up to +% 18 mmol/(gDW·h). + +modelFileName = 'ecoli_core_model.mat'; +modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models. +modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded +model = readCbModel(modelFileName); +model = changeRxnBounds(model, 'EX_glc(e)', -18, 'l'); +%% PROCEDURE +%% Identify essentail reactions: perform a gene knocked-out analysis. +% If the modelEcore is not able to grow when a certain gene is knocked-out, +% the name of that gene will be saved as an essential gene. Even if a very small +% growth is calculated, the model will be considered as not growing and the gene +% will be recorded in an 'essential_genes' vector. Here no growth is defined as +% growth lower than 0.000001. The remaining non-essentail genes will be stored +% in a 'non_EG' vector. +%% +[grRatio, grRateKO, grRateWT, delRxns,... + hasEffect] = singleGeneDeletion(model, 'FBA', model.genes); +essential_genes = []; +non_EG = []; +tol = 1e-6; +for n = 1:length(grRateKO) + if (grRateKO(n)=tol); +pFBAopt_Rxns = unique(regexprep(pFBAopt_Rxns, '_[f|b]$','')); +pFBAopt_Rxns = setdiff(pFBAopt_Rxns, pFBAEssentialRxns) +ELE_Rxns = model.rxns((abs(minFlux3)+abs(maxFlux3))<=tol); +ELE_Rxns = setdiff(ELE_Rxns, RxnMLEname); +ELE_Rxns = setdiff(ELE_Rxns, ZeroFluxRxns) +RxnELE = findRxnIDs(model, ELE_Rxns); +RxnOptima = findRxnIDs(model, pFBAopt_Rxns); +%% Classify the genes: +% The last step is to associate the genes that are related with each reaction. +% The main point of this is to classify the genes into the 5 different groups +% (Figure 3) and store them into different vectors: +% +% # _*Essential genes:_* metabolic genes necessary for growth in the given media +% ('essential_genes'). +% # _*pFBA optima: _*non-essential genes contributing to the optimal growth +% rate and minimum gene-associated flux ('OptimaGenes'). +% # _*Enzymatically less efficient (ELE): _*genes requiring more flux through +% enzymatic steps than alternative pathways that meet the same predicted growth +% rate ('ELEGenes'). +% # _*Metabolically less efficient (MLE):_* genes requiring a growth rate reduction +% if used ('MLEGenes'). +% # _*pFBA no-flux:_* genes that are unable to carry flux in the experimental +% conditions ('pFBAnoFluxGenes'). +% +% +% +% +% Figure 3: Gene classes retrieved through pFBA. +% +% Some genes may not fit in any of this 5 categories. These genes will be +% saved in a vetor called 'remainingGenes'. +%% +OptimaGenes = []; +restGenes = pFBAfluxGenes; +for i = 1:length(RxnOptima) + listGenes = find(RxnGMat(RxnOptima(i), :)); + for n = 1:length(listGenes) + pos = find(pFBAfluxGenes==listGenes(n)); + if pos + OptimaGenes = [OptimaGenes; model.genes(pFBAfluxGenes(pos))]; + restGenes(pos,1) = 0; + end + end +end +OptimaGenes = unique(OptimaGenes); +restGenes(restGenes==0) = []; + +ELEGenes = []; +restGenes2 = restGenes; +for i = 1:length(RxnELE) + listGenes = find(RxnGMat(RxnELE(i), :)); + for n = 1:length(listGenes) + pos = find(restGenes==listGenes(n)); + if pos + ELEGenes = [ELEGenes; model.genes(restGenes(pos))]; + restGenes2(pos, 1) = 0; + end + end +end +ELEGenes = unique(ELEGenes); +restGenes2(restGenes2==0) = []; + +MLEGenes = []; +finalRemainingGenes = restGenes2; +for i = 1:length(RxnMLE) + listGenes = find(RxnGMat(RxnMLE(i),:)); + for n = 1:length(listGenes) + pos = find(restGenes2==listGenes(n)); + if pos + MLEGenes = [MLEGenes; model.genes(restGenes2(pos))]; + finalRemainingGenes(pos, 1) = 0; + end + end +end +MLEGenes = unique(MLEGenes); +finalRemainingGenes(finalRemainingGenes==0) = []; + +remainingGenes = []; +for n = 1:length(finalRemainingGenes) + remainingGenes = [remainingGenes; model.genes(finalRemainingGenes(n))]; +end +%% *Print results:* + +essential_genes +OptimaGenes +ELEGenes +MLEGenes +pFBAnoFluxGenes +%% TIMING +% The tutorial runs in a few minutes. +%% References +% [1] Lewis et al. Omic data from evolved E. coli are consistent with computed +% optimal growth from genome-scale models. _Mol Syst Biol. _6:390 (2010). +% +% [2] Orth, J., Fleming, R.M., Palsson B. Ø. Reconstruction and Use of Microbial +% Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational +% Guide. _EcoSal Plus._ 4(1) (2010). +% +% [3] Thiele, I., et al. A community-driven global reconstruction of human +% metabolism. _Nat Biotechnol_. 31(5):419-425 (2013). \ No newline at end of file diff --git a/analysis/test/tutorial_pFBA.mlx b/analysis/test/tutorial_pFBA.mlx new file mode 100644 index 00000000..5f04a5c1 Binary files /dev/null and b/analysis/test/tutorial_pFBA.mlx differ