checkRxn Checks which reactants in a reaction that can be synthesized and which products that can be consumed. This is primarily for debugging reactions which cannot have flux model a model structure rxn the id of one reaction to check cutoff minimal flux for successful production/consumption (optional, default 10^-7) revDir true if the reaction should be reversed (optional, default false) printReport print a report (optional, default true) report reactants array with reactant indexes canMake boolean array, true if the corresponding reactant can be synthesized by the rest of the metabolic network products array with product indexes canConsume boolean array, true if the corresponding product can be consumed by the rest of the metabolic network Usage: report=checkRxn(model,rxn,cutoff,revDir,printReport)
0001 function report=checkRxn(model,rxn,cutoff,revDir,printReport) 0002 % checkRxn 0003 % Checks which reactants in a reaction that can be synthesized and which 0004 % products that can be consumed. This is primarily for debugging 0005 % reactions which cannot have flux 0006 % 0007 % model a model structure 0008 % rxn the id of one reaction to check 0009 % cutoff minimal flux for successful production/consumption (optional, 0010 % default 10^-7) 0011 % revDir true if the reaction should be reversed (optional, default 0012 % false) 0013 % printReport print a report (optional, default true) 0014 % 0015 % report 0016 % reactants array with reactant indexes 0017 % canMake boolean array, true if the corresponding reactant can 0018 % be synthesized by the rest of the metabolic network 0019 % products array with product indexes 0020 % canConsume boolean array, true if the corresponding product can 0021 % be consumed by the rest of the metabolic network 0022 % 0023 % Usage: report=checkRxn(model,rxn,cutoff,revDir,printReport) 0024 0025 rxn=char(rxn); 0026 if nargin<3 0027 cutoff=10^-7; 0028 end 0029 if nargin<4 0030 revDir=false; 0031 end 0032 if isempty(cutoff) 0033 cutoff=10^-7; 0034 end 0035 if nargin<5 0036 printReport=true; 0037 end 0038 0039 [I, rxnID]=ismember(rxn,model.rxns); 0040 0041 if ~I 0042 EM='Reaction ID not found'; 0043 dispEM(EM); 0044 end 0045 0046 if revDir==false 0047 report.reactants=find(model.S(:,rxnID)<0); 0048 report.products=find(model.S(:,rxnID)>0); 0049 else 0050 report.reactants=find(model.S(:,rxnID)>0); 0051 report.products=find(model.S(:,rxnID)<0); 0052 end 0053 report.canMake=false(numel(report.reactants),1); 0054 report.canConsume=false(numel(report.products),1); 0055 0056 %Remove this field as it would give an annoying note otherwise 0057 if isfield(model,'rxnComps') 0058 model=rmfield(model,'rxnComps'); 0059 end 0060 0061 %There are several ways to do this. Here I choose to add the reactions one 0062 %by one and checking their bounds. This might not be optimal 0063 for i=1:numel(report.reactants) 0064 [tempModel, testRxn]=addExchangeRxns(model,'out',report.reactants(i)); 0065 tempModel=setParam(tempModel,'obj',testRxn,1); 0066 sol=solveLP(tempModel); 0067 if sol.f>cutoff 0068 report.canMake(i)=true; 0069 else 0070 if printReport==true 0071 fprintf(['Failed to make ' model.metNames{report.reactants(i)} '[' model.comps{model.metComps(report.reactants(i))} ']\n']); 0072 end 0073 end 0074 end 0075 0076 for i=1:numel(report.products) 0077 [tempModel, testRxn]=addExchangeRxns(model,'in',report.products(i)); 0078 tempModel=setParam(tempModel,'obj',testRxn,1); 0079 sol=solveLP(tempModel); 0080 if sol.f>cutoff 0081 report.canConsume(i)=true; 0082 else 0083 if printReport==true 0084 fprintf(['Failed to consume ' model.metNames{report.products(i)} '[' model.comps{model.metComps(report.products(i))} ']\n']); 0085 end 0086 end 0087 end 0088 end