runINIT
Generates a model using the INIT algorithm, based on proteomics and/or
transcriptomics and/or metabolomics and/or metabolic tasks. This is the
original implementation, which is now replaced with ftINIT.
model a reference model structure
rxnScores a vector of scores for the reactions in the model.
Positive scores are reactions to keep and negative
scores are reactions to exclude (opt, default all 0.0)
presentMets cell array with unique metabolite names that the model
should produce (opt, default [])
essentialRxns cell array of reactions that are essential and that
have to be in the resulting model. This is normally
used when fitting a model to task (see fitTasks) (opt,
default [])
prodWeight a score that determines the value of having
net-production of metabolites. This is a way of having
a more functional network as it provides a reason for
including bad reactions for connectivity reasons. This
score is for each metabolite, and the sum of these weights
and the scores for the reactions is what is optimized
(opt, default 0.5)
allowExcretion true if excretion of all metabolites should be allowed.
This results in fewer reactions being considered
dead-ends, but all reactions in the resulting model may
not be able to carry flux. If this is "false" then the
equality constraints are taken from model.b. If the
input model lacks exchange reactions then this should
probably be "true", or a large proportion of the model
would be excluded for connectivity reasons
(opt, default false)
noRevLoops true if reversible reactions should be constrained to
only carry flux in one direction. This prevents
reversible reactions from being wrongly assigned as
connected (the forward and backward reactions can form a
loop and therefore appear connected), but it makes the
problem significantly more computationally intensive to
solve (two more integer constraints per reversible reaction)
(opt, default false)
params parameter structure as used by getMILPParams (opt,
default [])
outModel the resulting model structure
deletedRxns reactions which were deleted by the algorithm
metProduction array that indicates which of the
metabolites in presentMets that could be
produced
-2: metabolite name not found in model
-1: metabolite found, but it could not be produced
1: metabolite could be produced
fValue objective value (sum of (the negative of)
reaction scores for the included reactions and
prodWeight*number of produced metabolites)
This function is the actual implementation of the algorithm. See
getINITModel for a higher-level function for model reconstruction. See
PLoS Comput Biol. 2012;8(5):e1002518 for details regarding the
implementation.
Usage: [outModel deletedRxns metProduction fValue]=runINIT(model,...
rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,...
noRevLoops,params)0001 function [outModel, deletedRxns, metProduction, fValue]=runINIT(model,rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,noRevLoops,params) 0002 % runINIT 0003 % Generates a model using the INIT algorithm, based on proteomics and/or 0004 % transcriptomics and/or metabolomics and/or metabolic tasks. This is the 0005 % original implementation, which is now replaced with ftINIT. 0006 % 0007 % model a reference model structure 0008 % rxnScores a vector of scores for the reactions in the model. 0009 % Positive scores are reactions to keep and negative 0010 % scores are reactions to exclude (opt, default all 0.0) 0011 % presentMets cell array with unique metabolite names that the model 0012 % should produce (opt, default []) 0013 % essentialRxns cell array of reactions that are essential and that 0014 % have to be in the resulting model. This is normally 0015 % used when fitting a model to task (see fitTasks) (opt, 0016 % default []) 0017 % prodWeight a score that determines the value of having 0018 % net-production of metabolites. This is a way of having 0019 % a more functional network as it provides a reason for 0020 % including bad reactions for connectivity reasons. This 0021 % score is for each metabolite, and the sum of these weights 0022 % and the scores for the reactions is what is optimized 0023 % (opt, default 0.5) 0024 % allowExcretion true if excretion of all metabolites should be allowed. 0025 % This results in fewer reactions being considered 0026 % dead-ends, but all reactions in the resulting model may 0027 % not be able to carry flux. If this is "false" then the 0028 % equality constraints are taken from model.b. If the 0029 % input model lacks exchange reactions then this should 0030 % probably be "true", or a large proportion of the model 0031 % would be excluded for connectivity reasons 0032 % (opt, default false) 0033 % noRevLoops true if reversible reactions should be constrained to 0034 % only carry flux in one direction. This prevents 0035 % reversible reactions from being wrongly assigned as 0036 % connected (the forward and backward reactions can form a 0037 % loop and therefore appear connected), but it makes the 0038 % problem significantly more computationally intensive to 0039 % solve (two more integer constraints per reversible reaction) 0040 % (opt, default false) 0041 % params parameter structure as used by getMILPParams (opt, 0042 % default []) 0043 % 0044 % outModel the resulting model structure 0045 % deletedRxns reactions which were deleted by the algorithm 0046 % metProduction array that indicates which of the 0047 % metabolites in presentMets that could be 0048 % produced 0049 % -2: metabolite name not found in model 0050 % -1: metabolite found, but it could not be produced 0051 % 1: metabolite could be produced 0052 % fValue objective value (sum of (the negative of) 0053 % reaction scores for the included reactions and 0054 % prodWeight*number of produced metabolites) 0055 % 0056 % This function is the actual implementation of the algorithm. See 0057 % getINITModel for a higher-level function for model reconstruction. See 0058 % PLoS Comput Biol. 2012;8(5):e1002518 for details regarding the 0059 % implementation. 0060 % 0061 % Usage: [outModel deletedRxns metProduction fValue]=runINIT(model,... 0062 % rxnScores,presentMets,essentialRxns,prodWeight,allowExcretion,... 0063 % noRevLoops,params) 0064 0065 if nargin<2 0066 rxnScores=zeros(numel(model.rxns),1); 0067 end 0068 if isempty(rxnScores) 0069 rxnScores=zeros(numel(model.rxns),1); 0070 end 0071 if nargin<3 || isempty(presentMets) 0072 presentMets={}; 0073 else 0074 presentMets=convertCharArray(presentMets); 0075 end 0076 if nargin<4 || isempty(essentialRxns) 0077 essentialRxns={}; 0078 else 0079 essentialRxns=convertCharArray(essentialRxns); 0080 end 0081 if nargin<5 || isempty(prodWeight) 0082 prodWeight=0.5; 0083 end 0084 if nargin<6 0085 allowExcretion=false; 0086 end 0087 if nargin<7 0088 noRevLoops=false; 0089 end 0090 if nargin<8 0091 params=[]; 0092 end 0093 0094 if numel(presentMets)~=numel(unique(presentMets)) 0095 EM='Duplicate metabolite names in presentMets'; 0096 dispEM(EM); 0097 end 0098 0099 %Default is that the metabolites cannot be produced 0100 if ~isempty(presentMets) 0101 metProduction=ones(numel(presentMets),1)*-2; 0102 presentMets=upper(presentMets); 0103 pmIndexes=find(ismember(presentMets,upper(model.metNames))); 0104 metProduction(pmIndexes)=-1; %Then set that they are at least found 0105 else 0106 metProduction=[]; 0107 pmIndexes=[]; 0108 end 0109 0110 %The model should be in the reversible format and all relevant exchange 0111 %reactions should be open 0112 if isfield(model,'unconstrained') 0113 EM='Exchange metabolites are still present in the model. Use simplifyModel if this is not intended'; 0114 dispEM(EM,false); 0115 end 0116 0117 %The irreversible reactions that are essential must have a flux and are 0118 %therefore not optimized for using MILP, which reduces the problem size. 0119 %However, reversible reactions must have a flux in one direction, so they 0120 %have to stay in the problem. The essentiality constraint on reversible 0121 %reactions is implemented in the same manner as for reversible reactions 0122 %when noRevLoops==true, but with the additional constraint that C ub=-1. 0123 %This forces one of the directions to be active. 0124 revRxns=find(model.rev~=0); 0125 essentialReversible=find(ismember(model.rxns(revRxns),essentialRxns)); 0126 essentialRxns=intersect(essentialRxns,model.rxns(model.rev==0)); 0127 0128 %Convert the model to irreversible 0129 irrevModel=convertToIrrev(model); 0130 rxnScores=[rxnScores;rxnScores(model.rev==1)]; 0131 %These are used if noRevLoops is true 0132 if noRevLoops==true 0133 forwardIndexes=find(model.rev~=0); 0134 backwardIndexes=(numel(model.rxns)+1:numel(irrevModel.rxns))'; 0135 else 0136 %Then they should only be used for essential reversible reactions 0137 forwardIndexes=revRxns(essentialReversible); 0138 backwardIndexes=essentialReversible+numel(model.rxns); 0139 end 0140 0141 %Get the indexes of the essential reactions and remove them from the 0142 %scoring vector 0143 essentialIndex=find(ismember(irrevModel.rxns,essentialRxns)); 0144 rxnScores(essentialIndex)=[]; 0145 0146 %Go through each of the presentMets (if they exist) and modify the S matrix 0147 %so that each reaction which produces any of them also produces a 0148 %corresponding fake metabolite and the opposite in the reverse direction. 0149 0150 %This is to deal with the fact that there is no compartment info regarding 0151 %the presentMets. This modifies the irrevModel structure, but that is fine 0152 %since it's the model structure that is returned. 0153 if any(pmIndexes) 0154 irrevModel.metNames=upper(irrevModel.metNames); 0155 metsToAdd.mets=strcat({'FAKEFORPM'},num2str(pmIndexes)); 0156 metsToAdd.metNames=metsToAdd.mets; 0157 metsToAdd.compartments=irrevModel.comps{1}; 0158 0159 %There is no constraints on the metabolites yet, since maybe not all of 0160 %them could be produced 0161 irrevModel=addMets(irrevModel,metsToAdd); 0162 end 0163 0164 %Modify the matrix 0165 for i=1:numel(pmIndexes) 0166 %Get the matching mets 0167 I=ismember(irrevModel.metNames,presentMets(pmIndexes(i))); 0168 0169 %Find the reactions where any of them are used. 0170 [~, K, L]=find(irrevModel.S(I,:)); 0171 0172 %This ugly loop is to avoid problems if a metabolite occurs several 0173 %times in one reaction 0174 KK=unique(K); 0175 LL=zeros(numel(KK),1); 0176 for j=1:numel(KK) 0177 LL(j)=sum(L(K==KK(j))); 0178 end 0179 irrevModel.S(numel(irrevModel.mets)-numel(pmIndexes)+i,KK)=LL; 0180 end 0181 0182 %Some nice to have numbers 0183 nMets=numel(irrevModel.mets); 0184 nRxns=numel(irrevModel.rxns); 0185 nEssential=numel(essentialIndex); 0186 nNonEssential=nRxns-nEssential; 0187 nonEssentialIndex=setdiff(1:nRxns,essentialIndex); 0188 S=irrevModel.S; 0189 0190 %Add so that each non-essential reaction produces one unit of a fake 0191 %metabolite 0192 temp=sparse(1:nRxns,1:nRxns,1); 0193 temp(essentialIndex,:)=[]; 0194 S=[S;temp]; 0195 0196 %Add another set of reactions (will be binary) which also produce these 0197 %fake metabolites, but with a stoichiometry of 1000 0198 temp=sparse(1:nNonEssential,1:nNonEssential,1000); 0199 temp=[sparse(nMets,nNonEssential);temp]; 0200 S=[S temp]; 0201 0202 %Add reactions for net-production of (real) metabolites 0203 if prodWeight~=0 0204 temp=[speye(nMets-numel(pmIndexes))*-1;sparse(nNonEssential+numel(pmIndexes),nMets-numel(pmIndexes))]; 0205 S=[S temp]; 0206 %To keep the number of reactions added like this 0207 nNetProd=nMets-numel(pmIndexes); 0208 else 0209 nNetProd=0; 0210 end 0211 0212 %Add constraints so that reversible reactions can only be used in one 0213 %direction. This is done by adding the fake metabolites A, B, C for each 0214 %reversible reaction in the following manner 0215 % forward: A + .. => ... backwards: B + ... => ... int1: C => 1000 A int2: 0216 % C => 1000 B A ub=999.9 B ub=999.9 C lb=-1 int1 and int2 are binary 0217 if any(forwardIndexes) 0218 nRevBounds=numel(forwardIndexes); 0219 0220 %Add the A metabolites for the forward reactions and the B metabolites 0221 %for the reverse reactions 0222 I=speye(numel(irrevModel.rxns))*-1; 0223 temp=[I(forwardIndexes,:);I(backwardIndexes,:)]; 0224 0225 %Padding 0226 temp=[temp sparse(size(temp,1),size(S,2)-numel(irrevModel.rxns))]; 0227 0228 %Add the int1 & int2 reactions that produce A and B 0229 temp=[temp speye(nRevBounds*2)*1000]; 0230 0231 %And add that they also consume C 0232 temp=[temp;[sparse(nRevBounds,size(S,2)) speye(nRevBounds)*-1 speye(nRevBounds)*-1]]; 0233 0234 %Add the new reactions and metabolites 0235 S=[S sparse(size(S,1),nRevBounds*2)]; 0236 S=[S;temp]; 0237 else 0238 nRevBounds=0; 0239 end 0240 0241 %Add so that the essential reactions must have a small flux and that the 0242 %binary ones (and net-production reactions) may have zero flux. The integer 0243 %reactions for reversible reactions have [0 1] 0244 prob.blx=[irrevModel.lb;zeros(nNonEssential+nNetProd+nRevBounds*2,1)]; 0245 prob.blx(essentialIndex)=max(0.1,prob.blx(essentialIndex)); 0246 0247 %Add so that the binary ones and net-production reactions can have at the 0248 %most flux 1.0 0249 prob.bux=[irrevModel.ub;ones(nNonEssential+nNetProd+nRevBounds*2,1)]; 0250 0251 %Add that the fake metabolites must be produced in a small amount and that 0252 %the A and B metabolites for reversible reactions can be [0 999.9] and C 0253 %metabolites [-1 0] 0254 prob.blc=[irrevModel.b(:,1);ones(nNonEssential,1);zeros(nRevBounds*2,1);ones(nRevBounds,1)*-1]; 0255 0256 %Add that normal metabolites can be freely excreted if 0257 %allowExcretion==true, and that the fake ones can be excreted 1000 units at 0258 %most. C metabolites for essential reversible reactions should have an 0259 %upper bound of -1. If noRevLoops is false, then add this constraint for 0260 %all the reactions instead. 0261 if noRevLoops==true 0262 revUB=zeros(nRevBounds,1); 0263 revUB(essentialReversible)=-1; 0264 else 0265 revUB=ones(nRevBounds,1)*-1; 0266 end 0267 if allowExcretion==true 0268 metUB=inf(nMets,1); 0269 else 0270 metUB=irrevModel.b(:,min(size(irrevModel.b,2),2)); 0271 end 0272 prob.buc=[metUB;ones(nNonEssential,1)*1000;ones(nRevBounds*2,1)*999.9;revUB]; 0273 0274 %Add objective coefficients for the binary reactions. The negative is used 0275 %since we're minimizing. The negative is taken for the prodWeight as well, 0276 %in order to be consistent with the syntax that positive scores are good 0277 prob.c=[zeros(nRxns,1);rxnScores;ones(nNetProd,1)*prodWeight*-1;zeros(nRevBounds*2,1)]; 0278 prob.a=S; 0279 0280 % adapt problem structure for cobra-style solver 0281 prob.c=[prob.c;zeros(size(prob.a,1),1)]; 0282 prob.A=[prob.a -speye(size(prob.a,1))]; 0283 prob.b=zeros(size(prob.a,1), 1); 0284 prob.ub=[prob.bux; prob.buc]; 0285 prob.osense=1; 0286 prob.csense=char(zeros(1,size(prob.a,1))); 0287 prob.csense(:)='E'; 0288 0289 %We still don't know which of the presentMets that can be produced. Go 0290 %through them, force production, and see if the problem can be solved 0291 for i=1:numel(pmIndexes) 0292 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=1; 0293 prob.lb=[prob.blx; prob.blc]; 0294 res=optimizeProb(prob,params); 0295 isFeasible=checkSolution(res); 0296 if ~isFeasible 0297 %Reset the constraint again 0298 prob.blc(numel(irrevModel.mets)-numel(pmIndexes)+i)=0; 0299 else 0300 %Metabolite produced 0301 metProduction(pmIndexes(i))=1; 0302 end 0303 end 0304 prob.lb=[prob.blx; prob.blc]; 0305 0306 %Add that the binary reactions may only take integer values. 0307 prob.vartype = repmat('C', 1, size(prob.A, 2)); 0308 allInt=[(nRxns+1):(nRxns+nNonEssential) size(S,2)-nRevBounds*2+1:size(S,2)]; 0309 prob.vartype(allInt) = 'B'; 0310 0311 % solve problem 0312 res=optimizeProb(prob,params); 0313 0314 %Problem should not be infeasible, but it is possible that the time limit 0315 %was reached before finding any solutions. 0316 if ~checkSolution(res) 0317 if strcmp(res.origStat, 'TIME_LIMIT') 0318 EM='Time limit reached without finding a solution. Try increasing the TimeLimit parameter.'; 0319 else 0320 EM='The problem is infeasible'; 0321 end 0322 dispEM(EM); 0323 end 0324 0325 fValue=res.obj; 0326 0327 %Get all reactions used in the irreversible model 0328 usedRxns=(nonEssentialIndex(res.full(nRxns+1:nRxns+nNonEssential)<0.1))'; 0329 0330 %Map to reversible model IDs 0331 usedRxns=[usedRxns(usedRxns<=numel(model.rxns));revRxns(usedRxns(usedRxns>numel(model.rxns))-numel(model.rxns))]; 0332 0333 %Then get the ones that are not used in either direction or is essential 0334 I=true(numel(model.rxns),1); 0335 I(usedRxns)=false; 0336 I(essentialIndex)=false; 0337 deletedRxns=model.rxns(I); 0338 outModel=removeReactions(model,I,true,true); 0339 end