gapFind

PURPOSE ^

gapFind Identifies all blocked metabolites (anything downstream of a gap)

SYNOPSIS ^

function [allGaps,rootGaps,downstreamGaps] = gapFind(model,findNCgaps,verbFlag)

DESCRIPTION ^

gapFind Identifies all blocked metabolites (anything downstream of a gap) 
in a model.  MILP algorithm that finds gaps that may be missed by simple 
inspection of the S matrix. To find every gap in a model, change the rxn
bounds on all exchange reactions to allow uptake of every metabolite.

 [allGaps,rootGaps,downstreamGaps] = gapFind(model,findNCgaps,verbFlag)

INPUT
 model             a COBRA model

OPTIONAL INPUTS
 findNCgaps        find no consupmption gaps as well as no production gaps
                   (default false)   
 verbFlag          verbose flag (default false)

OUTPUTS
 allGaps           all gaps found by GapFind
 rootGaps          all root no production (and consumption) gaps
 downstreamGaps    all downstream gaps

 based on:
 Kumar, V. et al. BMC Bioinformatics. 2007 Jun 20;8:212.

 solve problem:
   max ||xnp||
       s.t. S(i,j)*v(j) >= e*w(i,j)        S(i,j) > 0, j in IR
            S(i,j)*v(j) <= M*w(i,j)        S(i,j) > 0, j in IR
            S(i,j)*v(j) >= e - M(1-w(i,j)) S(i,j) ~= 0, j in R
            S(i,j)*v(j) <= M*w(i,j)        S(i,j) ~= 0, j in R
            ||w(i,j)|| >= xnp(i)
            lb <= v <= ub
            S*v >= 0
            xnp(i) = {0,1}
            w(i,j) = {0,1}

 reformulated for COBRA MILP as:
   max sum(xnp(:))
       s.t. S*v >= 0   (or = 0 if findNCgaps = true)               (1)
            S(i,j)*v(j) - e*w(i,j) >= 0     S(i,j) > 0, j in IR    (2)
            S(i,j)*v(j) - M*w(i,j) <= 0     S(i,j) > 0, j in IR    (3)
            S(i,j)*v(j) - M*w(i,j) >= e-M   S(i,j) ~= 0, j in R    (4)
            S(i,j)*v(j) - M*w(i,j) <= 0     S(i,j) ~= 0, j in R    (5)
            sum(w(i,:)) - xnp(i) >= 0                              (6)
            lb <= v <= ub
            xnp and w are binary variables, v are continuous


 Jeff Orth 7/6/09

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [allGaps,rootGaps,downstreamGaps] = gapFind(model,findNCgaps,verbFlag)
0002 %gapFind Identifies all blocked metabolites (anything downstream of a gap)
0003 %in a model.  MILP algorithm that finds gaps that may be missed by simple
0004 %inspection of the S matrix. To find every gap in a model, change the rxn
0005 %bounds on all exchange reactions to allow uptake of every metabolite.
0006 %
0007 % [allGaps,rootGaps,downstreamGaps] = gapFind(model,findNCgaps,verbFlag)
0008 %
0009 %INPUT
0010 % model             a COBRA model
0011 %
0012 %OPTIONAL INPUTS
0013 % findNCgaps        find no consupmption gaps as well as no production gaps
0014 %                   (default false)
0015 % verbFlag          verbose flag (default false)
0016 %
0017 %OUTPUTS
0018 % allGaps           all gaps found by GapFind
0019 % rootGaps          all root no production (and consumption) gaps
0020 % downstreamGaps    all downstream gaps
0021 %
0022 % based on:
0023 % Kumar, V. et al. BMC Bioinformatics. 2007 Jun 20;8:212.
0024 %
0025 % solve problem:
0026 %   max ||xnp||
0027 %       s.t. S(i,j)*v(j) >= e*w(i,j)        S(i,j) > 0, j in IR
0028 %            S(i,j)*v(j) <= M*w(i,j)        S(i,j) > 0, j in IR
0029 %            S(i,j)*v(j) >= e - M(1-w(i,j)) S(i,j) ~= 0, j in R
0030 %            S(i,j)*v(j) <= M*w(i,j)        S(i,j) ~= 0, j in R
0031 %            ||w(i,j)|| >= xnp(i)
0032 %            lb <= v <= ub
0033 %            S*v >= 0
0034 %            xnp(i) = {0,1}
0035 %            w(i,j) = {0,1}
0036 %
0037 % reformulated for COBRA MILP as:
0038 %   max sum(xnp(:))
0039 %       s.t. S*v >= 0   (or = 0 if findNCgaps = true)               (1)
0040 %            S(i,j)*v(j) - e*w(i,j) >= 0     S(i,j) > 0, j in IR    (2)
0041 %            S(i,j)*v(j) - M*w(i,j) <= 0     S(i,j) > 0, j in IR    (3)
0042 %            S(i,j)*v(j) - M*w(i,j) >= e-M   S(i,j) ~= 0, j in R    (4)
0043 %            S(i,j)*v(j) - M*w(i,j) <= 0     S(i,j) ~= 0, j in R    (5)
0044 %            sum(w(i,:)) - xnp(i) >= 0                              (6)
0045 %            lb <= v <= ub
0046 %            xnp and w are binary variables, v are continuous
0047 %
0048 %
0049 % Jeff Orth 7/6/09
0050 
0051 if nargin < 2
0052     findNCgaps = false;
0053 end
0054 if nargin < 3
0055     verbFlag = false;
0056 end
0057 
0058 M = length(model.rxns); %this was set to 100 in GAMS GapFind implementation
0059 N = length(model.mets);
0060 R = model.rev ~= 0; %reversible reactions
0061 R_index = find(R);
0062 IR = model.rev == 0; %irreversible reactions
0063 IR_index = find(IR);
0064 e = 0.0001;
0065 S = model.S;
0066 lb = model.lb;
0067 ub = model.ub;
0068 
0069 % MILPproblem
0070 %  A      LHS matrix
0071 %  b      RHS vector
0072 %  c      Objective coeff vector
0073 %  lb     Lower bound vector
0074 %  ub     Upper bound vector
0075 %  osense Objective sense (-1 max, +1 min)
0076 %  csense Constraint senses, a string containting the constraint sense for
0077 %         each row in A ('E', equality, 'G' greater than, 'L' less than).
0078 %  vartype Variable types
0079 %  x0      Initial solution
0080 
0081 % initialize MILP fields
0082 
0083 % get number of rows and cols for each constraint
0084 %rows
0085 m_c1 = N; %number of metabolites
0086 m_c2 = length(find(S(:,IR) > 0)); %number of Sij>0 in irreverisible reactions
0087 m_c3 = m_c2;
0088 m_c4 = length(find(S(:,R))); %number of Sij>0 in reversible reactions
0089 m_c5 = m_c4;
0090 m_c6 = N; %number of xnp (metabolites)
0091 %columns
0092 n_v = M; %number of reactions
0093 n_wij_IR = m_c2;
0094 n_wij_R = m_c4;
0095 n_xnp = N;
0096 
0097 % LHS matrix A
0098 
0099 % constraint 1
0100 A = [S sparse(m_c1,(n_wij_IR+n_wij_R+n_xnp))];
0101 
0102 % constraint 2
0103 % create Sij IR matrix and wij IR matrix
0104 Sij_IR = sparse(m_c2,n_v);
0105 wij_IR = sparse(m_c6,n_wij_IR);
0106 row = 1;
0107 for i = 1:length(IR_index)
0108     rxn_index = IR_index(i);
0109     met_index = find(S(:,rxn_index) > 0);
0110     for j = 1:length(met_index)
0111         Sij_IR(row,rxn_index) = S(met_index(j),rxn_index);
0112         wij_IR(met_index(j),row) = 1;
0113         row = row + 1;
0114     end
0115 end  
0116 
0117 A = [A ; Sij_IR -e*speye(m_c2,n_wij_IR) sparse(m_c2,(n_wij_R+n_xnp))];
0118 
0119 % constraint 3
0120 A = [A ; Sij_IR -M*speye(m_c3,n_wij_IR) sparse(m_c3,(n_wij_R+n_xnp))];
0121 
0122 % constraint 4
0123 % create Sij R matrix
0124 Sij_R = sparse(m_c4,n_v);
0125 wij_R = sparse(m_c6,n_wij_R);
0126 row = 1;
0127 for i = 1:length(R_index)
0128     rxn_index = R_index(i);
0129     met_index = find(S(:,rxn_index) ~= 0);
0130     for j = 1:length(met_index)
0131         Sij_R(row,rxn_index) = S(met_index(j),rxn_index);
0132         wij_R(met_index(j),row) = 1;
0133         row = row + 1;
0134     end
0135 end  
0136 
0137 A = [A ; Sij_R sparse(m_c4,n_wij_IR) -M*speye(m_c4,n_wij_R) sparse(m_c4,n_xnp)];
0138 
0139 % constraint 5
0140 A = [A ; Sij_R sparse(m_c5,n_wij_IR) -M*speye(m_c5,n_wij_R) sparse(m_c5,n_xnp)];
0141 
0142 % constraint 6
0143 A = [A ; sparse(m_c6,n_v) wij_IR wij_R -1*speye(m_c6,n_xnp)];
0144 
0145 % RHS vector b
0146 b = [zeros(m_c1+m_c2+m_c3,1);(e-M)*ones(m_c4,1);zeros(m_c5+m_c6,1)];
0147 
0148 % objective coefficient vector c
0149 c = [zeros(n_v+n_wij_IR+n_wij_R,1);ones(n_xnp,1)];
0150 
0151 % upper and lower bounds on variables (v,w,xnp)
0152 lb = [lb;zeros(n_wij_IR+n_wij_R+n_xnp,1)];
0153 ub = [ub;ones(n_wij_IR+n_wij_R+n_xnp,1)];
0154 
0155 % objective sense osense
0156 osense = -1; %want to maximize objective
0157 
0158 % constraint senses csense
0159 if findNCgaps
0160     csense(1:m_c1) = 'E';
0161 else
0162     csense(1:m_c1) = 'G';
0163 end
0164 csense((m_c1+1):(m_c1+m_c2)) = 'G';
0165 csense((m_c1+m_c2+1):(m_c1+m_c2+m_c3)) = 'L';
0166 csense((m_c1+m_c2+m_c3+1):(m_c1+m_c2+m_c3+m_c4)) = 'G';
0167 csense((m_c1+m_c2+m_c3+m_c4+1):(m_c1+m_c2+m_c3+m_c4+m_c5)) = 'L';
0168 csense((m_c1+m_c2+m_c3+m_c4+m_c5+1):(m_c1+m_c2+m_c3+m_c4+m_c5+m_c6)) = 'G';
0169 
0170 % variable types vartype
0171 vartype(1:n_v) = 'C';
0172 vartype((n_v+1):(n_v+n_wij_IR+n_wij_R+n_xnp)) = 'B';
0173 
0174 % inital solution x0
0175 x0 = [];
0176 
0177 
0178 % run COBRA MILP solver
0179 gapFindMILPproblem.A = A;
0180 gapFindMILPproblem.b = b;
0181 gapFindMILPproblem.c = c;
0182 gapFindMILPproblem.lb = lb;
0183 gapFindMILPproblem.ub = ub;
0184 gapFindMILPproblem.osense = osense;
0185 gapFindMILPproblem.csense = csense;
0186 gapFindMILPproblem.vartype = vartype;
0187 gapFindMILPproblem.x0 = x0;
0188 
0189 if verbFlag
0190     parameters.printLevel = 3; 
0191 else
0192     parameters.printLevel = 0;
0193 end
0194 
0195 solution = solveCobraMILP(gapFindMILPproblem,parameters);
0196 
0197 % get the list of gaps from MILP solution
0198 metsProduced = solution.full((n_v+n_wij_IR+n_wij_R+1):(n_v+n_wij_IR+n_wij_R+n_xnp),1);
0199 allGaps = model.mets(~metsProduced);
0200 rootGaps = findRootNPmets(model,findNCgaps); %identify root gaps using findRootNPmets
0201 downstreamGaps = allGaps(~ismember(allGaps,rootGaps));
0202 
0203 
0204 
0205

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