forked from gulermalaria/iPfal17
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmodel_curation_careyBMCG.m
194 lines (171 loc) · 6.99 KB
/
model_curation_careyBMCG.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
%% Setup environment
initCobraToolbox % this only needs to be run once
changeCobraSolver('gurobi5');
changeCobraSolver('gurobi5','MILP')
changeCobraSolver('gurobi5','LP')
%% load plata model
cd('C:\Users\mac9jc\Documents\MATLAB\curation')
plata = readCbModel('plata_orig.xml');
model = plata; % copy plata model for modifications
%% update apicoplast localization nomenclature (necessary for rewriting file at end)
% replace '_ap' with '[ap]'
model.mets = regexprep(model.mets,'_ap\>' ,'[ap]');
model.mets = regexprep(model.mets,'\[a\]\>' ,'[ap]');
model.rxns = regexprep(model.rxns,'_ap\>' ,'[ap]');
model.rxns = regexprep(model.rxns,'\[a\]\>' ,'[ap]');
%% deleting reactions
cd('C:\Users\mac9jc\Documents\MATLAB\curation\final')
file = table2cell(readtable('reaction_delete_edits_final.xls'));
for j = 1: (length(file(:,1)))
% check that row has reaction ID
if length(file{j,1}) <= 1
disp('no reaction name')
break
end
% remove old version of reaction (if present)
vec = file(j,:);
new1 = removeRxns_old(model, vec{1});
model = new1;
end
clearvars vec j new1 file
%% add/modify reactions
file = table2cell(readtable('edits_additions_or_modifications_final.xls'));
for j = 1:(length(file(:,1)))
% check that row has reaction ID
if length(file{j,1}) <= 1
disp('no reaction name')
break
end
% make edit
vec = file(j,:);
new1 = removeRxns_old(model, vec{1}); % remove old version
% add back correct reaction
if vec{8} == 0 % reaction is not reversible
revFlag = 'false';
else % reaction is reversible
revFlag = 'true';
end
lb = vec{9}; ub = vec{10};
rxnName = vec{1}; metaboliteList = vec{3}; objCoeff = vec{11};
subSystem = vec{7}; grRule = vec{4}; %geneNameList = vec(5);
stoichCoeffList= {};
if isnan(grRule) % if no genes, switch from nan to empty cell
grRule = {};
end
new = addReaction(new1,rxnName,...
metaboliteList,stoichCoeffList,revFlag,lb,ub,objCoeff,...
subSystem,grRule); %,geneNameList,systNameList,checkDuplicate)
model = new;
indx = (strcmp(rxnName,model.rxns));
reaction_ec = vec(13); reaction_notes = vec(14); reaction_ref = vec(15);
model.rxnECNumbers(indx) = reaction_ec;
model.rxnNotes(indx) = reaction_notes;
model.rxnReferences(indx) = reaction_ref;
end
clearvars vec j new1 file indx reaction_ec reaction_notes new
clearvars ub lb revFlag rxnName stoichCoeffList objCoeff metaboliteList
clearvars reaction_ref subSystem grRule
%% biomass modification
file = table2cell(readtable('biomass_edits_final.xls'));
for j = 1:length(file(:,1))
% check that row has reaction ID
if length(file{j,1}) <= 1
disp('no reaction name')
break
end
% make edit
vec = file(j,:);
new1 = removeRxns_old(model, vec{1}); % remove old version
% add back correct reaction
if vec{8} == 0 % reaction is not reversible
revFlag = 'false';
else % reaction is reversible
revFlag = 'true';
end
lb = vec{9}; ub = vec{10};
rxnName = vec{1}; metaboliteList = vec{3}; objCoeff = vec{11};
subSystem = vec{7}; grRule = vec{4}; %geneNameList = vec(5);
stoichCoeffList= {};
if isnan(grRule) % if no genes, switch from nan to empty cell
grRule = {};
end
new = addReaction(new1,rxnName,...
metaboliteList,stoichCoeffList,revFlag,lb,ub,objCoeff,...
subSystem,grRule); %,geneNameList,systNameList,checkDuplicate)
model = new;
indx = (strcmp(rxnName,model.rxns));
reaction_ec = vec(13); reaction_notes = vec(14); reaction_ref = vec(15);
model.rxnECNumbers(indx) = reaction_ec;
model.rxnNotes(indx) = reaction_notes;
model.rxnReferences(indx) = reaction_ref;
end
clearvars vec j new1 file indx reaction_ec reaction_notes new
clearvars ub lb revFlag rxnName stoichCoeffList objCoeff metaboliteList
clearvars reaction_ref subSystem grRule
% check model numbers
[r4, ~] = size(model.rxns);
[g4, ~] = size(model.genes);
j = 0; %ECs
for i = 1:length(model.rxnECNumbers)
if length(model.rxnECNumbers{i}) >1
j = j+1;
end
end
EC4 = j; % number of EC numbers in model
clearvars j vec i
%% a couple additional edits easier here for formatting
model = addReaction(model,'EX_lac_L(e)',{'lac_L[e]'},(-1),false);
model.lb(strcmp('EX_pyr(e)',model.rxns)) = -1000;
model = changeRxnBounds(model,'SUCD2_u6m_mt',0,'l');
model.mets = regexprep(model.mets,'_ap\>' ,'[ap]');
%% check duplicate reactions and growth before saving
[model,removed1] = checkDuplicateRxn(model,1);
if ~isempty(removed1)
warning('duplicate reactions, check removed1')
end
[model,removed2] = checkDuplicateRxn(model,2);
if ~isempty(removed2)
warning('duplicate reactions, check removed2')
end
clearvars removed1 removed2
opt = optimizeCbModel(model,'max',0,0);
if opt.f < .1
warning('fba of cobra <.1')
end
% [missing_mets,present_mets] = biomassPrecursorCheck(model);
% Doesnt work because can't add dm reactions for mets that have exchange
% reactions
%% convert to tiger model and check
addpath(genpath('C:/Users/mac9jc/Documents/MATLAB/tiger/'));
start_tiger
set_solver('gurobi');
set_solver_option('MaxTime',60*60); % max time a TIGER simulation is allowed to run
set_solver_option('IntFeasTol',1e-8); % cutoff number for interpreting a value as 0
pf_tiger = cobra_to_tiger(model); % make cobra model into tiger model, use defaults
opt2 = fba(pf_tiger);
if opt2.val < .1
warning('fba of tiger model <.1')
end
%% SAVE MODEL AS MATLAB
cd('C:\Users\mac9jc\Documents\MATLAB\curation\final')
save model_cobra_edited.mat model;
% Convert each COBRA model to a TIGER model % DO NOT DO THIS- PROMOTES ROUNDING ERRORS
% model_tiger_edited = cobra_to_tiger(model_cobra_edited);
%% SAVE AS XML
writeCbModel(model,'sbml','2017_pf_model.xml',{'c','e','fv','m','ap'},...
{'cytoplasm','extracellular','food_vacuole','mitochondria','apicoplast'} );
% %% save model without localization if needed for metdraw purposes
% model_wo_loc = model
% model_wo_loc.mets = regexprep(model_wo_loc.mets, '\[(m)\]', '\[c\]')
% model_wo_loc.mets = regexprep(model_wo_loc.mets, '\[(ap)\]', '\[c\]')
% model_wo_loc.metNames = regexprep(model_wo_loc.mets, '\[(m)\]', '\[c\]')
% model_wo_loc.mets = regexprep(model_wo_loc.mets, '\[(fv)\]', '\[c\]')
% model_wo_loc.metNames = regexprep(model_wo_loc.mets, '\[(fv)\]', '\[c\]')
% model_wo_loc.metNames = regexprep(model_wo_loc.mets, '\[(ap)\]', '\[c\]')
% model_wo_loc.subSystems = regexprep(model_wo_loc.subSystems,';','')
% model_wo_loc.subSystems = regexprep(model_wo_loc.subSystems,'\s(\w+)','')
% writeCbModel(model_wo_loc,'sbml','2017_pf_model_wo_loc.xml',{'c','e'},...
% {'cytoplasm','extracellular'} );
[~, genes_used_by_model] = findUsedGenes(model, model.genes);
genes_used_by_model = table(genes_used_by_model);
writetable(genes_used_by_model,'genes_used_by_model.csv')