-
Notifications
You must be signed in to change notification settings - Fork 0
/
Rawls_TIMBR_prediction.m
158 lines (133 loc) · 6.42 KB
/
Rawls_TIMBR_prediction.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
%%Running TIMBR Predictions
% Code Adapted from blais_timbr_prediction.m which can be found at
% https://github.com/csbl/ratcon1
clear all
close all
clc
%% TIMBR software requirements: COBRA toolbox and linear programming solver
% addpath('ncomm_blais_xls2model.m')
% addpath('ncomm_blais_model2irrev.m')
% addpath('timbr.m')
% Specify output directory for TIMBR predictions
path_timbr_directory = './timbr_directory/'; % Specify a folder for TIMBR predictions to be saved
%initCobraToolbox % this only needs to be run once after installation
% Choose a linear programming solver to use
changeCobraSolver('glpk'); % gurobi6 or glpk
%% Load rat-specific metabolic network
load('rno_cobra_load.mat')
% Reactions conisdered "off" have lb and ub equal to 0
off_in_rno = rno_cobra_load.rxns(...
rno_cobra_load.lb == 0 & rno_cobra_load.ub == 0);
% Remove reactions disabled in the model and set the objective to biomass
rno_cobra = changeObjective(removeRxns(...
rno_cobra_load,off_in_rno),'RCR99999');
%% Test biomass production with COBRA model
% Because exchange fluxes were formulated as fmol / cell / hour and
% biomass was formulated in units of fmol / cell,
% the growth rate is measured as growth per hour (1 / hour):
% Note: These values similar to but not equal to values
% reported in the manuscript. By default, exchange reaction boundary
% conditions are set to relaxed physiological constraints.
rno_growth_rate_cobra = optimizeCbModel(rno_cobra);
rno_growth_rate = rno_growth_rate_cobra.f;
rno_doubling_time = log(2) / (rno_growth_rate);
disp([num2str(rno_growth_rate), ' = growth rate of 0.0510 per hour'])
disp([num2str(rno_doubling_time), ' = doubling time of 13.6026 hours'])
%% Convert COBRA models into irreversible format
% Specify default production requirements for all TIMBR simulations
% Originally, we required that the model produced biomass at a
% growth rate of 1 / week with an ATP maintenance requirement.
% However, we found that it was too difficult to resolve
% how much TIMBR production scores reflected the global demand for
% biomass production vs biomarker production.
% obj_value_required = 1 / 7 / 24;
obj_value_required = 0;
% atp_value_required = 5000;
atp_value_required = 0;
% Set unconstrained reaction bounds to +/- 10^6 because the largest
% consumption rate is fairly close to the default 10^3 reaction bound.
default_bound = 1000000; % 10^6
rno_irrev_base = ncomm_blais_model2irrev(rno_cobra);
default_rxns_rno = rno_irrev_base.rxns(rno_irrev_base.ub == 1000);
rno_irrev_default = changeRxnBounds(rno_irrev_base, ...
default_rxns_rno, default_bound, 'u');
rno_irrev_atp = changeRxnBounds(rno_irrev_default,...
'RCR11017_f',atp_value_required,'l');
rno_irrev_biomass = changeRxnBounds(rno_irrev_atp,...
'RCR99999_f',obj_value_required,'l');
% Note: as mentioned above, TIMBR predictions were not performed while
% requiring biomass production.
rno_irrev = rno_irrev_biomass;
%% Determine maxmimum flux through each exchange reaction for iRno
[rno_exchange, rno_demand] = findExcRxns(rno_irrev);
rno_production = rno_irrev.rxns((cellfun(@length, ...
regexp(rno_irrev.rxns,'_r')) == 0) & rno_exchange);
rno_fva = zeros(length(rno_production),1);
for exchange_index = 1:length(rno_production)
rno_consumption = setdiff(intersect(...
regexprep(rno_production(exchange_index),'_f','_r'),...
rno_irrev.rxns),...
rno_production(exchange_index));
if (~isempty(rno_consumption))
rno_production_fba = optimizeCbModel(changeObjective(...
changeRxnBounds(changeRxnBounds(...
rno_irrev,rno_production,default_bound,'u'),...
rno_consumption,0,'b'),...
rno_production(exchange_index)));
rno_fva(exchange_index,1) = rno_production_fba.f;
else
rno_production_fba = optimizeCbModel(changeObjective(...
changeRxnBounds(...
rno_irrev,rno_production,default_bound,'u'),...
rno_production(exchange_index)));
rno_fva(exchange_index,1) = rno_production_fba.f;
end
disp(rno_production_fba)
end
%% Specify minimimum required flux for each exchange reaction
fba_threshold = 1e-4;
rno_production_ok = rno_fva > fba_threshold;
rno_production_id = rno_production(rno_production_ok,1);
rno_production_requirement = min(0.90 * rno_fva(rno_production_ok,1), 100);
rno_production_count = length(rno_production_id);
%% Load TIMBR reaction weights
% These files are generated by the R script:
% Rawls_TIMBRWeights.R
rno_timbr_weights_load = readtable([...
'Rawls_supplement_timbrweights.txt'],'Delimiter','\t');
[rno_model_has_weight, rno_model2weight_index] = ismember(...
rno_irrev.rxns, rno_timbr_weights_load.rxn_irreversible);
[rno_weight_in_model, rno_weight2model_index] = ismember(...
rno_timbr_weights_load.rxn_irreversible, rno_irrev.rxns);
if (all(rno_model_has_weight))
% first 3 columns are organism_id, rxn_id, and rxn_irrerversible
rno_timbr_weights = rno_timbr_weights_load(rno_model2weight_index,4:end);
rno_timbr_id = rno_timbr_weights.Properties.VariableNames';
rno_timbr_count = length(rno_timbr_id);
else
warning('reaction weights not specified for all reactions in rno')
end
%% Run TIMBR algorithm and save results as .txt files
% Estimate the global network demand of producing each exchange metabolite
% under treatment and control conditions for various compounds.
% Saved outputs will be analyzed in the R/Bioconductor script:
% Rawls_TIMBRAnalysis.R
for rno_timbr_index = 1:rno_timbr_count
rno_timbr_network_demand = zeros(rno_timbr_count, 1);
disp(rno_timbr_id{rno_timbr_index});
for rno_production_index = 1:rno_production_count
disp(rno_production_id(rno_production_index,1));
rno_timbr_network_demand(rno_production_index,1) = ...
timbr(rno_irrev, ...
rno_production_id(rno_production_index,1), ...
rno_production_requirement(rno_production_index,1), ...
rno_timbr_weights(:,rno_timbr_index));
end
rno_timbr_file_name = [path_timbr_directory ...
'ncomm_blais_timbr_' rno_timbr_id{rno_timbr_index} '.txt'];
rno_timbr_table = table(...
rno_production_id,...
rno_timbr_network_demand,...
'VariableNames',{'timbr_id' 'timbr_value'});
writetable(rno_timbr_table,rno_timbr_file_name,'Delimiter','\t');
end