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