/
SCOUR_Yeast_random.m
225 lines (194 loc) · 8.53 KB
/
SCOUR_Yeast_random.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
function SCOUR_Yeast_random(nT,cov,num_IC,rep)
warning('off','all')
rng('shuffle');
fprintf('Starting machine learning framework...\n');
%% Initial Preparation of Hynne yeast model
% Train on autogenerated data, test on Hynne yeast model
% yeast stoichiometric interactions
load('hynneSTM.mat');
Hynne.S = stm;
% Find Hynne mass action (MA) interactions
[row_HynneMA col_HynneMA] = find(Hynne.S < 0);
HynneMA = [row_HynneMA col_HynneMA];
unique_flux = unique(HynneMA(:,2));
count_1contMet = 1;
count_2contMet = 1;
for i = 1:length(unique_flux)
controller_met_idx = find(HynneMA(:,2) == unique_flux(i));
if length(HynneMA(controller_met_idx,1)) == 1
HynneMA_1contMet_List(count_1contMet,:) = [HynneMA(controller_met_idx,1)' unique_flux(i)];
count_1contMet = count_1contMet + 1;
elseif length(HynneMA(controller_met_idx,1)) == 2
HynneMA_2contMet_List(count_2contMet,:) = [HynneMA(controller_met_idx,1)' unique_flux(i)];
count_2contMet = count_2contMet + 1;
end
end
Hynne_prefix = 'hynne';
%% Identify 1 controller metabolite reactions
% List of Hynne interactions with one controller metabolite
true_1contMet_Hynne = [14,11;
17,14;
19,17;
20,19;
3,23;];
count = 1;
for regIdx = 1:length(HynneMA_1contMet_List)
for trueIdx = 1:size(true_1contMet_Hynne,1)
if isequal(HynneMA_1contMet_List(regIdx,:),true_1contMet_Hynne(trueIdx,:))
Hynne_1contMet_trueInRegIdx(count,1) = regIdx;
count = count + 1;
end
end
end
% Hynne
% Setup testing set and testing label
testingLabel_Hynne_1contMet = logical(zeros(size(HynneMA_1contMet_List,1),1));
testingLabel_Hynne_1contMet(Hynne_1contMet_trueInRegIdx) = 1;
predictedLabel_L2_Hynne_1contMet = logical(round(rand(length(testingLabel_Hynne_1contMet),1)))
% Accuracy, sensitivity, and specificity calculations
predictionAccuracy_Hynne_1contMet = sum(predictedLabel_L2_Hynne_1contMet==testingLabel_Hynne_1contMet)/length(testingLabel_Hynne_1contMet)
fp = 0;
fn = 0;
for k = 1:length(testingLabel_Hynne_1contMet)
if testingLabel_Hynne_1contMet(k) == 1 && predictedLabel_L2_Hynne_1contMet(k) == 0
fn = fn + 1;
elseif testingLabel_Hynne_1contMet(k) == 0 && predictedLabel_L2_Hynne_1contMet(k) == 1
fp = fp + 1;
end
end
tp = 0;
tn = 0;
for k = 1:length(testingLabel_Hynne_1contMet)
if testingLabel_Hynne_1contMet(k) == 1 && predictedLabel_L2_Hynne_1contMet(k) == 1
tp = tp + 1;
elseif testingLabel_Hynne_1contMet(k) == 0 && predictedLabel_L2_Hynne_1contMet(k) == 0
tn = tn + 1;
end
end
sensitivity_Hynne_1contMet = tp / (fn+tp)
specificity_Hynne_1contMet = tn / (tn+fp)
ppv_Hynne_1contMet = tp / (tp+fp);
npv_Hynne_1contMet = tp / (tn+fn);
%% Remove 1 controller metabolite reactions
% Remove from Hynne using predicted 1 controller metabolite
% interactions
predicted_Hynne_1contMet = [HynneMA_1contMet_List(find(predictedLabel_L2_Hynne_1contMet==1),:)];
%Hynne_fluxes_to_remove = unique([1 predicted_Hynne_1contMet(:,2)']);
Hynne_fluxes_to_remove = unique([1 11 14 17 19 23]);
if ~isequal(sort(Hynne_fluxes_to_remove),1:size(Hynne.S,2))
Hynne_regScheme_2contMet = createRegSchemeList_hynne(Hynne.S,Hynne_fluxes_to_remove);
end
if exist('Hynne_regScheme_2contMet','var')
%% Identify 2 controller metabolite reactions
% List of Hynne interactions with two controller metabolites
true_2contMet_Hynne = [2,3,3;
4,6,4;
8,9,7;
5,13,10;
12,15,12;
16,17,13;
18,19,16;
15,20,18;
20,21,20;
3,4,22;];
count = 1;
for regIdx = 1:length(Hynne_regScheme_2contMet)
for trueIdx = 1:size(true_2contMet_Hynne,1)
if isequal(Hynne_regScheme_2contMet(regIdx,:),true_2contMet_Hynne(trueIdx,:))
Hynne_2contMet_trueInRegIdx(count,1) = regIdx;
count = count + 1;
end
end
end
% Hynne
% Setup testing set and testing label
testingLabel_Hynne_2contMet = logical(zeros(size(Hynne_regScheme_2contMet,1),1));
if exist('Hynne_2contMet_trueInRegIdx','var')
testingLabel_Hynne_2contMet(Hynne_2contMet_trueInRegIdx) = 1;
end
predictedLabel_L2_Hynne_2contMet = logical(round(rand(length(testingLabel_Hynne_2contMet),1)))
% Accuracy, sensitivity, and specificity calculations
predictionAccuracy_Hynne_2contMet = sum(predictedLabel_L2_Hynne_2contMet==testingLabel_Hynne_2contMet)/length(testingLabel_Hynne_2contMet)
fp = 0;
fn = 0;
for k = 1:length(testingLabel_Hynne_2contMet)
if testingLabel_Hynne_2contMet(k) == 1 && predictedLabel_L2_Hynne_2contMet(k) == 0
fn = fn + 1;
elseif testingLabel_Hynne_2contMet(k) == 0 && predictedLabel_L2_Hynne_2contMet(k) == 1
fp = fp + 1;
end
end
tp = 0;
tn = 0;
for k = 1:length(testingLabel_Hynne_2contMet)
if testingLabel_Hynne_2contMet(k) == 1 && predictedLabel_L2_Hynne_2contMet(k) == 1
tp = tp + 1;
elseif testingLabel_Hynne_2contMet(k) == 0 && predictedLabel_L2_Hynne_2contMet(k) == 0
tn = tn + 1;
end
end
sensitivity_Hynne_2contMet = tp / (fn+tp)
specificity_Hynne_2contMet = tn / (tn+fp)
ppv_Hynne_2contMet = tp / (tp+fp);
npv_Hynne_2contMet = tp / (tn+fn);
%% Remove 2 controller metabolite reactions
% Remove from Hynne using predicted 2 controller metabolite
% interactions
predicted_Hynne_2contMet = [Hynne_regScheme_2contMet(find(predictedLabel_L2_Hynne_2contMet==1),:)];
%Hynne_fluxes_to_remove = unique([1 predicted_Hynne_1contMet(:,2)' predicted_Hynne_2contMet(:,3)']);
Hynne_fluxes_to_remove = unique([1 11 14 17 19 23 3 4 7 10 12 13 16 18 20 22]);
if ~isequal(sort(Hynne_fluxes_to_remove),1:size(Hynne.S,2))
[~,Hynne_regScheme_3contMet] = createRegSchemeList_hynne(Hynne.S,Hynne_fluxes_to_remove);
end
end
if exist('Hynne_regScheme_3contMet','var')
%% Identify 3 controller metabolite reactions
% List of Hynne interactions with three controller metabolites
true_3contMet_Hynne = [1 2 4,2;
3 6 22,5;
7 8 9,6;
9 10 12,15;
3 5 22,24];
count = 1;
for regIdx = 1:length(Hynne_regScheme_3contMet)
for trueIdx = 1:size(true_3contMet_Hynne,1)
if isequal(Hynne_regScheme_3contMet(regIdx,:),true_3contMet_Hynne(trueIdx,:))
Hynne_3contMet_trueInRegIdx(count,1) = regIdx;
count = count + 1;
end
end
end
% Hynne
% Setup testing set and testing label
testingLabel_Hynne_3contMet = logical(zeros(size(Hynne_regScheme_3contMet,1),1));
if exist('Hynne_3contMet_trueInRegIdx','var')
testingLabel_Hynne_3contMet(Hynne_3contMet_trueInRegIdx) = 1;
end
predictedLabel_L2_Hynne_3contMet = logical(round(rand(length(testingLabel_Hynne_3contMet),1)))
% Accuracy, sensitivity, and specificity calculations
predictionAccuracy_Hynne_3contMet = sum(predictedLabel_L2_Hynne_3contMet==testingLabel_Hynne_3contMet)/length(testingLabel_Hynne_3contMet)
fp = 0;
fn = 0;
for k = 1:length(testingLabel_Hynne_3contMet)
if testingLabel_Hynne_3contMet(k) == 1 && predictedLabel_L2_Hynne_3contMet(k) == 0
fn = fn + 1;
elseif testingLabel_Hynne_3contMet(k) == 0 && predictedLabel_L2_Hynne_3contMet(k) == 1
fp = fp + 1;
end
end
tp = 0;
tn = 0;
for k = 1:length(testingLabel_Hynne_3contMet)
if testingLabel_Hynne_3contMet(k) == 1 && predictedLabel_L2_Hynne_3contMet(k) == 1
tp = tp + 1;
elseif testingLabel_Hynne_3contMet(k) == 0 && predictedLabel_L2_Hynne_3contMet(k) == 0
tn = tn + 1;
end
end
sensitivity_Hynne_3contMet = tp / (fn+tp)
specificity_Hynne_3contMet = tn / (tn+fp)
ppv_Hynne_3contMet = tp / (tp+fp);
npv_Hynne_3contMet = tp / (tn+fn);
predicted_Hynne_3contMet = [Hynne_regScheme_3contMet(find(predictedLabel_L2_Hynne_3contMet==1),:)];
end
save(sprintf('yeast_random_results_IC-%02d_nT-%03d_cov-%02d_rep-%02d.mat',num_IC,nT,cov,rep));