0001 function [GeneClasses RxnClasses modelIrrevFM] = pFBA(model, varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043 if nargin < 2
0044 tol = 1e-6;
0045 GeneOption = 1;
0046 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
0073
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
0083 model_sav = model;
0084
0085
0086 [selExc,selUpt] = findExcRxns(model,0,0);
0087 tempmodel = changeRxnBounds(model,model.rxns(selExc),-1000,'l');
0088 tempmodel = changeRxnBounds(tempmodel,model.rxns(selExc),1000,'u');
0089 tempmodel = reduceModel(tempmodel,tol);
0090 Blocked_Rxns = setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''));
0091 model = removeRxns(model,setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$','')));
0092 Ind2Remove = find(~and(sum(full(model.rxnGeneMat),1),1));
0093 Blocked_genes = model.genes(Ind2Remove);
0094 model.genes(Ind2Remove)={'dead_end'};
0095
0096
0097 grRatio = singleGeneDeletion(model);
0098 grRatio(isnan(grRatio))=0;
0099 pFBAEssential = model.genes(grRatio<tol);
0100
0101 if nargout > 1
0102
0103 RxnRatio = singleRxnDeletion(model);
0104 RxnRatio(isnan(RxnRatio))=0;
0105 pFBAEssentialRxns = model.rxns(RxnRatio<tol);
0106 end
0107
0108
0109
0110 tempmodel = reduceModel(model,tol);
0111 ZeroFluxRxns = setdiff(model.rxns,regexprep(tempmodel.rxns,'_r$',''));
0112 model = removeRxns(model,ZeroFluxRxns);
0113
0114
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
0124 [ MinimizedFlux modelIrrevFM]= minimizeModelFlux_local(model,GeneOption);
0125
0126
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
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
0146 Ind2Remove = find(~and(sum(full(model.rxnGeneMat),1),1));
0147 ZeroFluxGenes = unique(model.genes(Ind2Remove));
0148
0149
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
0159 MLEGenes = setdiff(model.genes, [pFBAOptima;ZeroFluxGenes;ELEGenes]);
0160
0161
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
0173 pFBAOptima(ismember(pFBAOptima,pFBAEssential))=[];
0174
0175 if nargout > 1
0176
0177 pFBAopt_Rxns(ismember(pFBAopt_Rxns,pFBAEssentialRxns))=[];
0178 end
0179
0180
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
0222
0223
0224
0225 modelIrrev = convertToIrreversible(model);
0226
0227
0228 if nargin==1,GeneOption=0;
0229 end
0230 if GeneOption==0,
0231 modelIrrev.S(end+1,:) = ones(size(modelIrrev.S(1,:)));
0232 elseif GeneOption==1,
0233
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
0242 modelIrrev = addReaction(modelIrrev,'netFlux',{'fluxMeasure'},[-1],false,0,inf,0,'','');
0243
0244
0245 modelIrrev.c = zeros(length(modelIrrev.rxns),1);
0246 modelIrrev = changeObjective(modelIrrev, 'netFlux');
0247
0248
0249 MinimizedFlux = optimizeCbModel(modelIrrev,'min');
0250 end