growExpMatch run the growthExpMatch algorythm [solution]=growthExpMatch(model, KEGGFilename, compartment, iterations, dictionary, logFile, threshold) INPUTS model COBRA model structure KEGGFilename File name containing Kegg database (.lst file with list of reactions: each listing is reaction name followed by colon followed by the reaction formula) compartment [c] --> transport from cytoplasm [c] to extracellulat space [e] (default), [p] creates transport from [c] to [p] and from [p] to [c] iterations Number of iterations to run dictionary n x 2 cell array of metabolites names for conversion from Kegg ID's to the compound abbreviations from BiGG database (1st column is compound abb. (non-compartmenalized) and 2nd column is Kegg ID) Both columns are the same length. OPTINAL INPUTS logFile solution is printed in this file (name of reaction added and flux of that particular reaction) (Default = GEMLog.txt) threshold threshold number for biomass reaction; model is considered to be growing when the flux of the biomass reaction is above threshold. (Default = 0.05) OUTPUT solution MILP solution that consists of the continuous solution, integer solution, objective value, stat, full solution, and imported reactions %Procedure to run SMILEY: (1) obtain all input files (ie. model, CompAbr, and KeggID are from BiGG, KeggList is from Kegg website) (2) remove desired reaction from model with removeRxns, or set the model on a particular Carbon or Nitrogen source (3) create an SUX Matrix by using the function MatricesSUX = generateSUXMatrix(model,dictionary, KEGGFilename,compartment) (4) run it through SMILEY using [solution,b,solInt]=Smiley(MatricesSUX) (5) solution.importedRxns contains the solutions to all iterations MILPproblem A LHS matrix b RHS vector c Objective coeff vector lb Lower bound vector ub Upper bound vector osense Objective sense (-1 max, +1 min) csense Constraint senses, a string containting the constraint sense for each row in A ('E', equality, 'G' greater than, 'L' less than). vartype Variable types x0 Initial solution
0001 function [solution]=growthExpMatch(model, KEGGFilename, compartment, iterations, dictionary, logFile, threshold) 0002 %growExpMatch run the growthExpMatch algorythm 0003 % 0004 % [solution]=growthExpMatch(model, KEGGFilename, compartment, iterations, dictionary, logFile, threshold) 0005 % 0006 %INPUTS 0007 % model COBRA model structure 0008 % KEGGFilename File name containing Kegg database (.lst file with list of 0009 % reactions: each listing is reaction name followed by colon 0010 % followed by the reaction formula) 0011 % compartment [c] --> transport from cytoplasm [c] to extracellulat space 0012 % [e] (default), [p] creates transport from [c] to [p] 0013 % and from [p] to [c] 0014 % iterations Number of iterations to run 0015 % dictionary n x 2 cell array of metabolites names for conversion from 0016 % Kegg ID's to the compound abbreviations from BiGG database 0017 % (1st column is compound abb. (non-compartmenalized) and 0018 % 2nd column is Kegg ID) Both columns are the same length. 0019 % 0020 %OPTINAL INPUTS 0021 % logFile solution is printed in this file (name of reaction added and 0022 % flux of that particular reaction) (Default = GEMLog.txt) 0023 % 0024 % threshold threshold number for biomass reaction; model is considered 0025 % to be growing when the flux of the biomass reaction is 0026 % above threshold. (Default = 0.05) 0027 % 0028 %OUTPUT 0029 % solution MILP solution that consists of the continuous solution, integer 0030 % solution, objective value, stat, full solution, and 0031 % imported reactions 0032 % 0033 % 0034 %%Procedure to run SMILEY: 0035 %(1) obtain all input files (ie. model, CompAbr, and KeggID are from BiGG, KeggList is from Kegg website) 0036 %(2) remove desired reaction from model with removeRxns, or set the model 0037 % on a particular Carbon or Nitrogen source 0038 %(3) create an SUX Matrix by using the function MatricesSUX = 0039 % generateSUXMatrix(model,dictionary, KEGGFilename,compartment) 0040 %(4) run it through SMILEY using [solution,b,solInt]=Smiley(MatricesSUX) 0041 %(5) solution.importedRxns contains the solutions to all iterations 0042 % 0043 % MILPproblem 0044 % A LHS matrix 0045 % b RHS vector 0046 % c Objective coeff vector 0047 % lb Lower bound vector 0048 % ub Upper bound vector 0049 % osense Objective sense (-1 max, +1 min) 0050 % csense Constraint senses, a string containting the constraint sense for 0051 % each row in A ('E', equality, 'G' greater than, 'L' less than). 0052 % vartype Variable types 0053 % x0 Initial solution 0054 0055 % Based on IT 11/2008 0056 % Edited by xx 0057 0058 if nargin<2 0059 MatricesSUX= model; 0060 iterations=2; 0061 else 0062 MatricesSUX = generateSUXMatrix(model,dictionary, KEGGFilename,compartment); 0063 end 0064 if nargin <6 0065 logFile='GEMLog'; 0066 end 0067 if nargin <7 0068 threshold= 0.05; 0069 end 0070 0071 MatricesSUXori = MatricesSUX; 0072 0073 startKegg = find(MatricesSUX.MatrixPart ==2, 1, 'first'); 0074 stopKegg = find(MatricesSUX.MatrixPart ==2, 1, 'last'); 0075 lengthKegg = stopKegg -startKegg +1; 0076 startEx = find(MatricesSUX.MatrixPart ==3, 1, 'first'); 0077 stopEx = find(MatricesSUX.MatrixPart ==3, 1, 'last'); 0078 lengthEx = stopEx -startEx +1; 0079 if isempty(lengthEx) 0080 lengthEx = 0; 0081 end 0082 0083 vmax = 1000; 0084 0085 [a,b]=size(MatricesSUX.S); 0086 0087 cols = b + lengthKegg + lengthEx; 0088 rows = a + lengthKegg + lengthEx; 0089 0090 [i1,j1,s1] = find(MatricesSUX.S); 0091 0092 % Kegg 0093 %add rows 0094 i2a=(a+1:a+lengthKegg)'; 0095 j2a=(startKegg:stopKegg)'; 0096 s2a=ones(lengthKegg,1); 0097 %add cols 0098 i2b=(a+1:a+lengthKegg)'; 0099 j2b=(b+1:b+lengthKegg)'; 0100 s2b=-vmax*ones(lengthKegg,1); 0101 0102 %Exchange 0103 %add rows 0104 i3a=(a+lengthKegg+1:rows)'; 0105 j3a=(startEx:stopEx)'; 0106 s3a=ones(lengthEx,1); 0107 %add cols 0108 i3b=(a+lengthKegg+1:rows)'; 0109 j3b=(b+lengthKegg+1:cols)'; 0110 s3b=-vmax*ones(lengthEx,1); 0111 0112 if isempty(MatricesSUX.c); 0113 disp('No biomass reaction found'); 0114 else 0115 biomass_rxn_loc = find(MatricesSUX.c); 0116 end 0117 0118 i = [i1;i2a;i2b;i3a;i3b;rows+1]; 0119 j=[j1;j2a;j2b;j3a;j3b;biomass_rxn_loc]; 0120 s=[s1;s2a;s2b;s3a;s3b;1]; 0121 0122 MatricesSUX.A = sparse(i,j,s); 0123 MatricesSUX.b = zeros(size(MatricesSUX.A,1),1); 0124 %%% Set the threshold to 0.05 %%% 0125 MatricesSUX.b(size(MatricesSUX.A,1),1) = threshold; 0126 0127 MatricesSUX.cOuter = zeros(size(MatricesSUX.A,2),1); 0128 MatricesSUX.cOuter(b+1:b+lengthKegg)=1; 0129 MatricesSUX.cOuter(b+lengthKegg+1:cols)=2.1; 0130 0131 MatricesSUX.csense(1:a)='E'; 0132 MatricesSUX.csense(a+1:rows)='L'; 0133 MatricesSUX.csense(rows+1) = 'G'; 0134 0135 MatricesSUX.lb(b+1:cols)=0; 0136 MatricesSUX.ub(b+1:cols)=1; 0137 MatricesSUX.rxns(b+1:cols)=strcat(MatricesSUX.rxns(startKegg:stopEx),'dummy'); 0138 0139 MatricesSUX.MatrixPart(b+1:cols)=4; 0140 Int= find(MatricesSUX.MatrixPart>=4); 0141 0142 spy(MatricesSUX.A) 0143 0144 x0=zeros(size(MatricesSUX.A,2),1); 0145 0146 solInt=zeros(length(Int),1); 0147 0148 %%% setting the MILPproblem %%% 0149 vartype(b+1:cols) = 'B'; 0150 vartype(1:b) = 'C'; 0151 0152 MILPproblem.A = MatricesSUX.A; 0153 MILPproblem.b = MatricesSUX.b; 0154 MILPproblem.c = MatricesSUX.cOuter; 0155 MILPproblem.lb = MatricesSUX.lb; 0156 MILPproblem.ub = MatricesSUX.ub; 0157 MILPproblem.csense = MatricesSUX.csense; 0158 MILPproblem.osense = 1; 0159 MILPproblem.vartype = vartype; 0160 MILPproblem.x0 = x0; 0161 0162 for i = 1: iterations 0163 solution = solveCobraMILP(MILPproblem); 0164 0165 if(solution.obj~=0) 0166 solInt(:,i+1)=solution.int; 0167 printSolutionGEM(MatricesSUX, solution,logFile,i); 0168 MILPproblem.A(rows+i+1,j2b(1):cols) = solInt(:,i+1)'; 0169 MILPproblem.b(rows+i+1) = sum(solInt(:,i+1))-.1; 0170 MILPproblem.csense(rows+i+1) = 'L'; 0171 % save([logFile '_solution_' num2str(i)]); 0172 end 0173 0174 0175 solution.importedRxns = findImportedReactions(solInt, MatricesSUX); 0176 tmp=find(solution.cont); 0177 for j=1:length(tmp) 0178 if(tmp(j)>=Int(1)) 0179 MILPproblem.ub(tmp(j))=solution.cont(tmp(j));%%% 0180 end 0181 end 0182 if (solution.stat~=1) 0183 break 0184 end 0185 end 0186 printSolutionGEM(MatricesSUX, solution); 0187 0188 function importedRxns = findImportedReactions(solInt, MatricesSUX) 0189 importedRxns= {}; 0190 stopModel = find((MatricesSUX.MatrixPart ==1), 1, 'last'); 0191 for i = 1: size(solInt,2)-1 0192 [x,y] = find(solInt(:, i+1)); 0193 for j = 1: size(x) 0194 importedRxns{i, j} = MatricesSUX.rxns(stopModel + x(j)); 0195 MatricesSUX.rxns(stopModel + x(j)); 0196 end 0197 end 0198 0199