-
Notifications
You must be signed in to change notification settings - Fork 1
/
SteadyComPOACplex.m
338 lines (329 loc) · 13.6 KB
/
SteadyComPOACplex.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
function [POAtable, fluxRange, Stat, GRvector] = SteadyComPOACplex(modelCom,options,solverParam)
%Pairwise POA for community model at community steady-state for a range of growth rates
%[POAtable, fluxRange, Stat, GRvector] = SteadyComPOACplex(modelCom,options,solverParam)
%
%INPUT
% modelCom A community COBRA model structure with the following extra fields:
% (the following fields are required - others can be supplied)
% S Stoichiometric matrix
% b Right hand side
% c Objective coefficients
% lb Lower bounds
% ub Upper bounds
% (at least one of the below two is needed)
% infoCom structure containing community reaction info
% (returned along with the community model created with createCommModel)
% indCom the index structure corresponding to infoCom
%
% options (optional) option structure with the following fields:
% GRmax Maximum growth rate of the model (default to be found
% SteadyComCplex.m)
% optGRpercent A vector of percentages. Perform POA at these percents
% of max. growth rate respectively (Default = 99.99)
% optBMpercent Only consider solutions that yield at least a certain
% percentage of the optimal biomass (Default = 95)
% rxnNameList List of reactions (IDs or .rxns) to be analyzed.
% Use a (N_rxns + N_organism) x K matrix for POA of K
% linear combinations of fluxes and/or abundances
% (Default = biomass reaction of each species)
% pairList Pairs in rxnNameList to be analyzed. N_pair by 2 array of:
% - indices referring to the rxns in rxnNameList, e.g.,
% [1 2] to analyze rxnNameList{1} vs rxnNameList{2}
% - rxn names which are members of rxnNameList, e.g.,
% {'EX_glc-D(e)','EX_ac(e)'}
% If not supplied, analyze all K(K-1) pairs from the K
% targets in rxnNameList.
% symmetric Used only when pairList is not supplied. Avoid running
% symmetric pairs (e.g. j vs k and k vs j)
% Nstep Number of steps for fixing one flux at a value between
% the min. and the max. possible fluxes. Default 10.
% Can also be a vector indicating the fraction of intermediate value to be analyzed
% e.g. [0 0.5 1] means computing at minFlux, 0.5(minFlux + maxFlux) and maxFlux
% NstepScale Used only when Nstep is a single number.
% -'lin' for a linear (uniform) scale of step size
% -'log' for a log scaling of the step sizes
% fluxRange Flux range for each entry in rxnNameList. K x 2 matrix.
% Defaulted to be found by SteadyComFVACplex.m
% (other parameters)
% savePOA Must be non-empty. The filename to save the POA results
% (default 'POAtmp/POA')
% threads > 1 for explicitly stating the no. of threads used,
% 0 or -1 for using all available threads. Default 1.
% verbFlag Verbose output. 0 or 1.
% loadModel String of filename to be loaded. If non-empty, load the
% cplex model ('loadModel.mps'), basis ('loadModel.bas')
% and parameters ('loadModel.prm').
% May add also other parameters in SteadyComCplex for calculating the maximum growth rate.
%
%OUTPUT
% (N_gr = numel(optGRpercent))
% POAtable K x K cells.
% Each (i,i)-cell contains a Nstep x 1 x N_gr matrix of
% the fluxes at which rxnNameList{i} is fixed
% (i,j)-cell contains a Nstep x 2 x N_gr matrix,
% (p,:,q)-entry being the range of rxnNameList{j}
% when rxnNameList{i} is fixed at POAtable{i,i}(p,1,q)
% at growth rate = GRvector(q)
% fluxRange K x 2 x N_gr matrix of flux range for each entry in rxnNameList
% Stat K x K x N_gr structure array with fields:
% -'cor': the slope from linear regression between the
% fluxes of a pair
% -'r2': the corresponding coefficient of determination (R-square)
% GRvector Vector of growth rates being analyzed
%% Initialization
%check required fields for community model
if ~isfield(modelCom,'indCom')
if ~isfield(modelCom,'infoCom') || ~isstruct(modelCom.infoCom) || ...
~all(isfield(modelCom.infoCom,{'spBm','EXcom','EXsp','spAbbr','rxnSps','metSps'}))
error('infoCom must be provided for calculating the max. community growth rate.\n');
end
%get useful reaction indices
modelCom.indCom = infoCom2indCom(modelCom);
end
%get paramters
if ~exist('options', 'var')
options = struct();
end
if ~exist('solverParam', 'var') || isempty(solverParam)
%default Cplex parameters
solverParam = getCobraComParams('CplexParam');
end
param2get = {'GRmax', 'optGRpercent', 'Nstep', 'GRfx','BMmaxLB','BMmaxUB',...
'rxnNameList', 'verbFlag', 'savePOA','loadModel'};
eval(sprintf('[%s] = getCobraComParams(param2get, options, modelCom);', ...
strjoin(param2get, ',')...
)...
);
if isempty(savePOA)
error('A non-empty file name must be provided to save the POA results.');
end
[feasTol, ~] = getCobraSolverParams('LP',{'feasTol'; 'optTol'}, solverParam);
if isfield(solverParam,'simplex') && isfield(solverParam.simplex, 'tolerances')...
&& isfield(solverParam.simplex.tolerances,'feasibility')
%override the feasTol in CobraSolverParam if given in solverParam
feasTol = solverParam.simplex.tolerances.feasibility;
else
%otherwise use the feasTol in COBRA toolbox
solverParam.simplex.tolerances.feasibility = feasTol;
end
init = true;
if exist([savePOA '_MasterModel.mat'], 'file')
load([savePOA '_MasterModel.mat'],'LPstart','LPmodel','GRvector','kDisp','idRow');
LP = Cplex('POA');
LP.Model = LPmodel;
LP.Start = LPstart;
LP = setCplexParam(LP, solverParam);
if ~isfield(options, 'BMmaxLB')
options.BMmaxLB = LP.Model.lhs(idRow);
end
if ~isfield(options, 'BMmaxUB')
options.BMmaxUB = LP.Model.rhs(idRow);
end
init = false;
end
if numel(Nstep) > 1
Nstep = numel(Nstep);
end
if ischar(rxnNameList)
rxnNameList = {rxnNameList};
end
if iscell(rxnNameList)
Ncheck = numel(rxnNameList);
else
Ncheck = size(rxnNameList,2);
end
if init
addRow = false;
%get maximum growth rate
if isempty(GRmax)
if exist('Cplex.p','file') == 6
options.minNorm = false;
[~, result,LP] = SteadyComCplex(modelCom, options,solverParam);
else
%need further achitecture for using COBRA solver
warning('Currently support Cplex only.');
return
end
if strcmp(result.stat,'infeasible')
%infeasible model
warning('Model is infeasible.');
POAtable = cell(Ncheck);
fluxRange = NaN(Ncheck,2);
Stat = repmat(struct('cor',[],'r2',[]),Ncheck,Ncheck);
GRvector = NaN(numel(optGRpercent), 1);
return
end
GRmax = result.GRmax;
idRow = size(LP.Model.A,1);
else
%If GRmax is given, BMmaxLB and BMmaxUB should be included in options in this case to ensure feasibility
if ~isempty(loadModel)
[m, n] = size(modelCom.S);
nRxnSp = sum(modelCom.indCom.rxnSps > 0); %number of species-specific rxns
nSp = numel(modelCom.indCom.spBm); %number of species
% load solution if given and growth rate is known
LP = Cplex('fluxSampling');
LP.readModel([loadSol '.mps']);
LP.readBasis([loadSol '.bas']);
LP.readParam([loadSol '.prm']);
fprintf('Load model ''%s'' successfully.\n', loadModel);
addRow = true;
if size(LP.Model.A,1) > m + 2*nRxnSp + nSp
[ynRow,idRow] = ismember(sparse(ones(nSp,1),n+1:n+nSp,ones(nSp,1),1,n+nSp),...
LP.Model.A(m+2*nRxnSp+nSp+1:end,1:n+nSp),'rows');
if ynRow
idRow = m + 2*nRxnSp + nSp + idRow;
end
addRow = ~ynRow;
end
else
%get LP using SteadyComCplex if only growth rate is given
options2 = options;
options2.LPonly = true;
[~, ~, LP] = SteadyComCplex(modelCom, options2, solverParam);
%no constraint on total biomass using LPonly option
addRow = true;
end
end
if addRow
%add a row for constraining the sum of biomass if not exist
LP.addRows(BMmaxLB, ...
sparse(ones(1, nSp), n + 1: n + nSp, ones(1, nSp), 1, size(LP.Model.A,2)),...
BMmaxUB, 'UnityBiomass');
idRow = size(LP.Model.A,1);
else
%using BMmaxLB and BMmaxUB stored in the LP if not given in options
if ~isfield(options,'BMmaxLB') %take from LP if not supplied
BMmaxLB = LP.Model.lhs(idRow);
end
if ~isfield(options,'BMmaxUB') %take from LP if not supplied
BMmaxUB = LP.Model.rhs(idRow);
end
LP.Model.lhs(idRow) = BMmaxLB;
%not allow the max. biomass to exceed the one at max growth rate,
%can happen if optBMpercent < 100. May dismiss this constraint or
%manually supply BMmaxUB in the options if sum of biomass should be variable
LP.Model.rhs(idRow) = BMmaxUB;
end
%set Cplex parameters
LP = setCplexParam(LP, solverParam);
LP.Model.A = updateLPcom(modelCom, GRmax, GRfx, [], LP.Model.A, []);
LP.Model.sense = 'minimize';
LP.Model.obj(:) = 0;
LP.solve();
%check and adjust for feasibility
dev = checkSolFeas(LP);
kBMadjust = 0;
BMmaxLB = LP.Model.lhs(idRow);
while (~isfield(LP.Solution, 'x') || dev > feasTol) && kBMadjust < 10
kBMadjust = kBMadjust + 1;
%row of biomass constraint should at the end
LP.Model.lhs(idRow) = BMmaxLB * (1 - feasTol/(11 - kBMadjust));
LP.solve();
dev = checkSolFeas(LP);
if verbFlag
fprintf('BMmax adjusment: %d\n',kBMadjust);
end
end
if (~isfield(LP.Solution, 'x') || dev > feasTol)
error('Model not feasible.')
end
if ~isfield(options, 'BMmaxLB')
options.BMmaxLB = LP.Model.lhs(idRow);
end
if ~isfield(options, 'BMmaxUB')
options.BMmaxUB = LP.Model.rhs(idRow);
end
GRvector = GRmax * optGRpercent/100;
if numel(optGRpercent) == 1
kDisp = 2;
else
d = max(GRvector(2:end) - GRvector(1:end-1));
if d < 1
kDisp = abs(floor(log10(abs(d))));
else
kDisp = 0;
end
end
end
nVar = size(LP.Model.A,2);
%handle objective matrix
if isnumeric(rxnNameList)
if size(rxnNameList,1) >= size(modelCom.S,2) && size(rxnNameList,1) <= nVar
%it is a matrix of objective vectors
objList = [sparse(rxnNameList); sparse(nVar - size(rxnNameList,1), size(rxnNameList,2))];
elseif size(rxnNameList,1) == 1 || size(rxnNameList,2) == 1
%reaction index
objList = sparse(rxnNameList, 1:numel(rxnNameList), ones(numel(rxnNameList),1),...
nVar, max(size(rxnNameList)));
else
error('Invalid numerical input of rxnNameList.');
end
elseif iscell(rxnNameList)
objList = sparse(nVar, numel(rxnNameList));
for jRxnName = 1:numel(rxnNameList)
rJ = findRxnIDs(modelCom,rxnNameList{jRxnName});
if ~all(rJ)
error('Invalid names in rxnNameList');
end
objList(rJ,jRxnName) = 1;
end
else
error('Invalid input of rxnNameList');
end
options.rxnNameList = objList;
clear objList
if isempty(savePOA)
%always use save option to reduce memory need
savePOA = ['POAtmp' filesep 'POA'];
end
directory = strsplit(savePOA,filesep);
if numel(directory) > 1
%not saving in the current directory. Check existence
directory = strjoin(directory(1:end-1),filesep);
if ~exist(directory,'dir')
mkdir(directory);
end
end
if init
LPstart = LP.Start;
LPmodel = LP.Model;
save([savePOA '_MasterModel.mat'],'LPstart','LPmodel','GRvector', 'kDisp','idRow','GRmax');
end
for j = 1:numel(GRvector)
optionsJ = options;
optionsJ.GR = GRvector(j);
optionsJ.savePOA = sprintf(['%s_GR%.' num2str(kDisp) 'f'], savePOA, GRvector(j));
%better reset the model to ensure feasibility
if j > 1
load([savePOA '_MasterModel.mat'],'LPstart','LPmodel');
LP.Model = LPmodel;
LP.Start = LPstart;
LP = setCplexParam(LP, solverParam);
end
clear LPmodel LPstart
SteadyComPOAgrCplex(modelCom,optionsJ,solverParam,LP);
end
%collect output from save file
POAtable = cell(Ncheck, Ncheck);
Stat = repmat(struct('cor',0,'r2',0), [Ncheck, Ncheck, numel(GRvector)]);
fluxRange = zeros(Ncheck, 2, numel(GRvector));
for j = 1:numel(GRvector)
data = load(sprintf(['%s_GR%.' num2str(kDisp) 'f.mat'],savePOA,GRvector(j)), 'POAtable', 'fluxRange', 'Stat');
for p = 1:Ncheck
for q = 1:Ncheck
if isempty(data.POAtable{p,q})
POAtable{p,q} = [];
elseif size(data.POAtable{p,q}, 1) == 1
%single point, no range
POAtable{p,q}(:, :, j) = repmat(data.POAtable{p,q}, Nstep, 1);
else
POAtable{p,q}(:, :, j) = data.POAtable{p,q};
end
Stat(p,q,j).cor = data.Stat(p,q).cor;
Stat(p,q,j).r2 = data.Stat(p,q).r2;
fluxRange(:,:, j) = data.fluxRange;
end
end
end
end