pFBA

PURPOSE ^

Parsimoneous enzyme usage Flux Balance Analysis - method that optimizes

SYNOPSIS ^

function [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, varargin)

DESCRIPTION ^

 Parsimoneous enzyme usage Flux Balance Analysis - method that optimizes
 the user's objective function and then minimizes the flux through the
 model and subsequently classifies each gene by how it contributes to the
 optimal solution. See Lewis, et al. Mol Syst Bio doi:10.1038/msb.2010.47

 Inputs:
 model                 COBRA model

 varargin:
 'geneoption'          1 = only minimize the sum of the flux through
                       gene-associated fluxes (default), 0 = minimize the
                       sum of all fluxes in the network
 'tol'                 tolerance (default = 1e-6)
 'map'                 map structure from readCbMap.m (no map written if empty
 'mapoutname'          File Name for map 
 'skipclass'           1 = Don't classify genes and reactions. Only return
                       model with the minimium flux set as upper bound.
                       0 = classify genes and reactions (default). 

 Output:
 GeneClasses           Structure with fields for each gene class
 RxnsClasses           Structure with fields for each reaction class
 modelIrrevFM          Irreversible model used for minimizing flux with
                       the minimum flux set as a flux upper bound

 
 
 ** note on maps: Red (6) = Essential reactions, Orange (5) = pFBA optima
       reaction, Yellow (4) = ELE reactions, Green (3) = MLE reactions,
       blue (2) = zero flux reactions, purple (1) = blocked reactions,
       black (0) = not classified
 
 Example:
 [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, 'geneoption',0, 'tol',1e-7)

 by Nathan Lewis Aug 25, 2010

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, varargin)
0002 % Parsimoneous enzyme usage Flux Balance Analysis - method that optimizes
0003 % the user's objective function and then minimizes the flux through the
0004 % model and subsequently classifies each gene by how it contributes to the
0005 % optimal solution. See Lewis, et al. Mol Syst Bio doi:10.1038/msb.2010.47
0006 %
0007 % Inputs:
0008 % model                 COBRA model
0009 %
0010 % varargin:
0011 % 'geneoption'          1 = only minimize the sum of the flux through
0012 %                       gene-associated fluxes (default), 0 = minimize the
0013 %                       sum of all fluxes in the network
0014 % 'tol'                 tolerance (default = 1e-6)
0015 % 'map'                 map structure from readCbMap.m (no map written if empty
0016 % 'mapoutname'          File Name for map
0017 % 'skipclass'           1 = Don't classify genes and reactions. Only return
0018 %                       model with the minimium flux set as upper bound.
0019 %                       0 = classify genes and reactions (default).
0020 %
0021 % Output:
0022 % GeneClasses           Structure with fields for each gene class
0023 % RxnsClasses           Structure with fields for each reaction class
0024 % modelIrrevFM          Irreversible model used for minimizing flux with
0025 %                       the minimum flux set as a flux upper bound
0026 %
0027 %
0028 %
0029 % ** note on maps: Red (6) = Essential reactions, Orange (5) = pFBA optima
0030 %       reaction, Yellow (4) = ELE reactions, Green (3) = MLE reactions,
0031 %       blue (2) = zero flux reactions, purple (1) = blocked reactions,
0032 %       black (0) = not classified
0033 %
0034 % Example:
0035 % [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, 'geneoption',0, 'tol',1e-7)
0036 %
0037 % by Nathan Lewis Aug 25, 2010
0038 
0039 
0040 
0041 %
0042 
0043 if nargin < 2
0044     tol = 1e-6;
0045     GeneOption = 1;
0046     map = []; % no map
0047     mapOutName = 'pFBA_map.svg';
0048     skipclass = 0; 
0049 end
0050 
0051 if mod(length(varargin),2)==0
0052     for i=1:2:length(varargin)-1
0053         switch lower(varargin{i})
0054             case 'geneoption', GeneOption = varargin{i+1};
0055             case 'tol', tol = varargin{i+1};
0056             case 'map', map = varargin{i+1};
0057             case 'mapoutname', mapOutName = varargin{i+1};
0058             case 'skipclass', skipclass = varargin{i+1};
0059             otherwise, options.(varargin{i}) = varargin{i+1};
0060         end
0061     end
0062 else
0063     error('Invalid number of parameters/values');
0064 end
0065 
0066 if ~exist('GeneOption','var'), GeneOption = 1;              end
0067 if ~exist('tol','var'), tol = 1e-6;                         end
0068 if ~exist('map','var'), map = [];                           end
0069 if ~exist('mapOutName','var'), mapOutName = 'pFBA_map.svg'; end
0070 if ~exist('skipclass','var'), skipclass = 0;                end
0071 
0072 if skipclass % skip the model reduction and gene/rxn classification
0073     % minimize the network flux
0074 FBAsoln = optimizeCbModel(model);
0075 model.lb(model.c==1) = FBAsoln.f;
0076 [ MinimizedFlux modelIrrevFM]= minimizeModelFlux_local(model,GeneOption);
0077 modelIrrevFM = changeRxnBounds(modelIrrevFM,'netFlux',MinimizedFlux.f,'b');
0078 GeneClasses = [];
0079 RxnClasses = [];
0080 else
0081 
0082 % save a copy of the inputted model
0083 model_sav = model;
0084 
0085 % Remove all blocked reactions
0086 [selExc,selUpt] = findExcRxns(model,0,0); % find and open up all exchanges
0087 tempmodel = changeRxnBounds(model,model.rxns(selExc),-1000,'l');
0088 tempmodel = changeRxnBounds(tempmodel,model.rxns(selExc),1000,'u');
0089 tempmodel = reduceModel(tempmodel,tol); % reduce the model to find blocked reactions
0090 Blocked_Rxns = setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''));
0091 model = removeRxns(model,setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''))); % remove blocked reactions
0092 Ind2Remove = find(~and(sum(full(model.rxnGeneMat),1),1));
0093 Blocked_genes = model.genes(Ind2Remove);
0094 model.genes(Ind2Remove)={'dead_end'}; % make sure genes that are unique to blocked reactions are tagged for removal
0095 
0096 % find essential genes
0097 grRatio = singleGeneDeletion(model);
0098 grRatio(isnan(grRatio))=0;
0099 pFBAEssential = model.genes(grRatio<tol);
0100 
0101 if nargout > 1
0102     % find essential reactions
0103 RxnRatio = singleRxnDeletion(model);
0104 RxnRatio(isnan(RxnRatio))=0;
0105 pFBAEssentialRxns = model.rxns(RxnRatio<tol);
0106 end
0107     
0108 
0109 % remove zero flux rxns
0110 tempmodel = reduceModel(model,tol);
0111 ZeroFluxRxns = setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''));
0112 model = removeRxns(model,ZeroFluxRxns);
0113 
0114 % find MLE reactions
0115 FBAsoln = optimizeCbModel(model);
0116 model.lb(model.c==1) = FBAsoln.f;
0117 [minFlux,maxFlux] = fluxVariability(model,100);
0118 for i = 1:length(minFlux)
0119     tmp(i,1) = max([abs(minFlux(i)) abs(maxFlux(i))])<tol;
0120 end
0121 MLE_Rxns = setdiff(model.rxns(tmp),ZeroFluxRxns);
0122 
0123 % minimize the network flux
0124 [ MinimizedFlux modelIrrevFM]= minimizeModelFlux_local(model,GeneOption);
0125 
0126 % separate pFBA optima rxns from ELE rxns
0127 modelIrrevFM = changeRxnBounds(modelIrrevFM,'netFlux',MinimizedFlux.f,'b');
0128 [minFlux,maxFlux] = fluxVariability(modelIrrevFM,100); 
0129 pFBAopt_Rxns = modelIrrevFM.rxns((abs(minFlux)+abs(maxFlux))>=tol);
0130 ELE_Rxns = modelIrrevFM.rxns((abs(minFlux)+abs(maxFlux))<=tol);
0131 pFBAopt_Rxns = unique(regexprep(pFBAopt_Rxns,'_[f|b]$',''));
0132 ELE_Rxns = unique(regexprep(ELE_Rxns,'_[f|b]$',''));
0133 ELE_Rxns = ELE_Rxns(~ismember(ELE_Rxns,pFBAopt_Rxns));
0134 ELE_Rxns = ELE_Rxns(~ismember(ELE_Rxns,MLE_Rxns));
0135 
0136 % determine pFBA optima genes
0137 pFBAopt_Rxns(ismember(pFBAopt_Rxns,'netFlux'))=[];
0138 [geneList]=findGenesFromRxns(model,pFBAopt_Rxns);
0139 geneList2 = {};
0140 for i = 1:length(geneList),
0141     geneList2(end+1:end+length(geneList{i}),1) = columnVector( geneList{i});
0142 end
0143 pFBAOptima = unique(geneList2);
0144 
0145 % determine Zero Flux genes
0146 Ind2Remove = find(~and(sum(full(model.rxnGeneMat),1),1));
0147 ZeroFluxGenes = unique(model.genes(Ind2Remove));
0148 
0149 % determine ELE genes
0150 [geneList]=findGenesFromRxns(model,ELE_Rxns);
0151 geneList2 = {};
0152 for i = 1:length(geneList)
0153     geneList2(end+1:end+length(geneList{i}),1) = columnVector( geneList{i});
0154 end
0155 ELEGenes = unique(geneList2);
0156 ELEGenes = setdiff(ELEGenes,[pFBAOptima;ZeroFluxGenes]);
0157 
0158 % determine Met ineff genes
0159 MLEGenes = setdiff(model.genes, [pFBAOptima;ZeroFluxGenes;ELEGenes]);
0160 
0161 % clean up lists by removing non-genes
0162 pFBAOptima(~cellfun('isempty',regexp(pFBAOptima,'dead_end')))=[];
0163 ELEGenes(~cellfun('isempty',regexp(ELEGenes ,'dead_end')))=[];
0164 MLEGenes(~cellfun('isempty',regexp(MLEGenes,'dead_end')))=[];
0165 ZeroFluxGenes(~cellfun('isempty',regexp(ZeroFluxGenes,'dead_end')))=[];
0166 
0167 pFBAOptima(cellfun('isempty',pFBAOptima))=[];
0168 ELEGenes(cellfun('isempty',ELEGenes ))=[];
0169 MLEGenes(cellfun('isempty',MLEGenes))=[];
0170 ZeroFluxGenes(cellfun('isempty',ZeroFluxGenes))=[];
0171 
0172 % filter out essential genes from pFBA optima
0173 pFBAOptima(ismember(pFBAOptima,pFBAEssential))=[];
0174 
0175 if nargout > 1
0176 % filter out essential Rxns from pFBA optima
0177 pFBAopt_Rxns(ismember(pFBAopt_Rxns,pFBAEssentialRxns))=[];
0178 end
0179 
0180 % prepare output variables
0181 GeneClasses.pFBAEssential =pFBAEssential;
0182 GeneClasses.pFBAoptima = pFBAOptima;
0183 GeneClasses.ELEGenes = ELEGenes;
0184 GeneClasses.MLEGenes = MLEGenes;
0185 GeneClasses.ZeroFluxGenes = ZeroFluxGenes;
0186 GeneClasses.Blockedgenes = Blocked_genes;
0187 
0188 RxnClasses.Essential_Rxns = pFBAEssentialRxns;
0189 RxnClasses.pFBAOpt_Rxns = pFBAopt_Rxns;
0190 RxnClasses.ELE_Rxns = ELE_Rxns;
0191 RxnClasses.MLE_Rxns = MLE_Rxns;
0192 RxnClasses.ZeroFlux_Rxns = ZeroFluxRxns;
0193 RxnClasses.Blocked_Rxns = Blocked_Rxns;
0194 
0195 if ~isempty(map)
0196     MapVector = zeros(length(model_sav.rxns),1);
0197     MapVector(ismember(model_sav.rxns,Blocked_Rxns))= 1;
0198     MapVector(ismember(model_sav.rxns,ZeroFluxRxns))= 2;
0199     MapVector(ismember(model_sav.rxns,MLE_Rxns))= 3;
0200     MapVector(ismember(model_sav.rxns,ELE_Rxns))= 4;
0201     MapVector(ismember(model_sav.rxns,pFBAopt_Rxns))= 5;
0202     MapVector(ismember(model_sav.rxns,pFBAEssentialRxns))= 6;
0203     
0204     options.lb = 0;
0205     options.ub = 6;
0206     tmpCmap = hsv(18);
0207     tmpCmap = [tmpCmap([1,3,4,6,11,14],:); 0 0 0;];
0208     options.fileName = mapOutName;
0209     options.colorScale = flipud(round(tmpCmap*255));
0210     
0211     global CB_MAP_OUTPUT
0212     if ~exist('CB_MAP_OUTPUT', 'var') || isempty(CB_MAP_OUTPUT)
0213         changeCbMapOutput('svg');
0214     end
0215     drawFlux(map, model_sav, MapVector, options);
0216 end
0217 end
0218 end
0219 
0220 function [ MinimizedFlux modelIrrev]= minimizeModelFlux_local(model,GeneOption)
0221 % This function finds the minimum flux through the network and returns the
0222 % minimized flux and an irreversible model
0223 
0224 % convert model to irrev
0225 modelIrrev = convertToIrreversible(model);
0226 
0227 % add pseudo-metabolite to measure flux through network
0228 if nargin==1,GeneOption=0;
0229 end
0230 if GeneOption==0, % signal that you want to minimize the sum of all gene and non-gene associated fluxes
0231     modelIrrev.S(end+1,:) = ones(size(modelIrrev.S(1,:)));
0232 elseif GeneOption==1, % signal that you want to minimize the sum of only gene-associated fluxes
0233     %find all reactions which are gene associated
0234     Ind=find(sum(modelIrrev.rxnGeneMat,2)>0);
0235     modelIrrev.S(end+1,:) = zeros(size(modelIrrev.S(1,:)));
0236     modelIrrev.S(end,Ind) = 1;
0237 end
0238 modelIrrev.b(end+1) = 0;
0239 modelIrrev.mets{end+1} = 'fluxMeasure';
0240 
0241 % add a pseudo reaction that measures the flux through the network
0242 modelIrrev = addReaction(modelIrrev,'netFlux',{'fluxMeasure'},[-1],false,0,inf,0,'','');
0243 
0244 % set the flux measuring demand as the objective
0245 modelIrrev.c = zeros(length(modelIrrev.rxns),1);
0246 modelIrrev = changeObjective(modelIrrev, 'netFlux');
0247 
0248 % minimize the flux measuring demand (netFlux)
0249 MinimizedFlux = optimizeCbModel(modelIrrev,'min');
0250 end

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003