growthExpMatch

PURPOSE ^

growExpMatch run the growthExpMatch algorythm

SYNOPSIS ^

function [solution]=growthExpMatch(model, KEGGFilename, compartment, iterations, dictionary, logFile, threshold)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003