-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathpathsNum.m
executable file
·174 lines (160 loc) · 6.78 KB
/
pathsNum.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
function selectedPaths = pathsNum(L_DirAdjMatWn, L_DirAdjMatC, ...
L_DirAdjMatMets, pairs, source_mets, target_mets, startFromMin, D, ...
ApplyShortestDistanceOfSubsystems, ThrowErrorOnDViolation)
% This function finds the shortest distance from subsystem to subsystem.
% ATTENTION: THIS IS DIRECTION SPECIFIC !!!
%
% INPUTS
% - L_DirAdjMatWn: (Calculated in allpaths)
% This is like DirAdjMatWn, but it refers to connections of degree L.
% - L_DirAdjMatC: (Calculated in allpaths)
% This is like DirAdjMatC, but it refers to connections of degree L.
% - L_DirAdjMatMets: (Calculated in allpaths)
% This is like DirAdjMatC, but it explicitly lists all the participating
% metabolites for each connection.
% - pairs:
% - source_mets:
% - target_mets:
% - startFromMin:
% - D: (Defined in redGEM_150112)
%
% OUTPUTS
% - selectedPaths: This is a structure containing the information about the path
% from subsystem to subsystem, with fields:
% .sub2subNumP: Number of pathways from source subsystem to target
% subsystem
% .sub2subNumR: Total number of reactions from source subsystem to
% target subsystem
% .joiningReacs: The actual reactions
% .joiningPaths: The actual pathways
% .joiningMets: The actual metabolites
% .ModelReacs: All reactions per level (D)
% .ModelMets: All metabolites per level (D)
% direction = 'bothways'; % 'bothways' or 'eachdirection' here we specify if the shortest distance is for both ways (a to b and b to a) or for each direction
direction = ApplyShortestDistanceOfSubsystems;
numSS=length(unique(pairs));
sub2subNumP=zeros(numSS,numSS,D);
% Adding d_max for the cases where subsystems are not connected even after
% L iterations (it happens), causing the while loop to go look for indices
% larger than L in L_DirAdjMatWn (which is a tensor with 1:L as 3rd
% dimension index)
d_max = size(L_DirAdjMatWn,3);
if strcmp(startFromMin,'yes')
shortestLsub2subNumP=zeros(numSS,numSS);%find shortest L between SS
switch direction
case {'eachdirection'}
for k = 1:size(pairs,1)
d=1;
while d < d_max &&...
all(~any(L_DirAdjMatWn(source_mets{k},target_mets{k},d)))
d=d+1;
end
if d <= d_max
shortestLsub2subNumP(pairs(k,1),pairs(k,2))=d;
end
end
for k = 1:size(pairs,1)
d=1;
while d < d_max &&...
all(~any(L_DirAdjMatWn(target_mets{k},source_mets{k},d)))%switch direction!
d=d+1;
end
if d <= d_max
shortestLsub2subNumP(pairs(k,2),pairs(k,1))=d; %switch direction!
end
end
case {'bothways'}
for k = 1:size(pairs,1)
d=1;
while d <= d_max &&...
(all(~any(L_DirAdjMatWn(target_mets{k},source_mets{k},d))) &&...
all(~any(L_DirAdjMatWn(source_mets{k},target_mets{k},d))))
d=d+1;
end
if d <= d_max
shortestLsub2subNumP(pairs(k,1),pairs(k,2))=d;
shortestLsub2subNumP(pairs(k,2),pairs(k,1))=d; %switch direction!
end
end
end
elseif strcmp(startFromMin,'no')
shortestLsub2subNumP=ones(numSS,numSS);%find shortest L between SS
else
error('Wrong option for startFromMin')
end
% Try to find any subsystems that are not connected by L
[source_ix, target_ix] = find(shortestLsub2subNumP == 0);
if ~isempty(source_ix) || ~isempty(target_ix)
if strcmp(ThrowErrorOnDViolation,'error')
error('L value exceeded in matching subsystems :')
disp(numSS(source_ix));
fprintf('and\n');
disp(numSS(target_ix));
else
warning('L value exceeded in matching subsystems :')
disp(source_ix);
fprintf('and\n');
disp(target_ix);
end
end
clear direction d i
%% get the reactions and pathways joining the subsystems
% this section simply condenses all the met to met pathways(reactions) and represents
% them as subsystem to subsystem
joiningReacs=cell(numSS,numSS,D);
joiningMets=cell(numSS,numSS,D);
joiningPaths=cell(numSS,numSS,D);
ModelPaths=cell(D,1);
ModelReacs=cell(D,1);
ModelMets=cell(D,1);
for el=1:D
for k=1:size(pairs,1)
% Added a if condition, in the case that subsystems are not linked,
% d will stay at 0
if shortestLsub2subNumP(pairs(k,1),pairs(k,2)) > 0
d=el-1+shortestLsub2subNumP(pairs(k,1),pairs(k,2));
a=L_DirAdjMatC(source_mets{k}, target_mets{k},d);
emptyCells = cellfun('isempty',a);
a(emptyCells) = [];
uniqueReacs=unique([a{:}],'stable');
joiningReacs{pairs(k,1),pairs(k,2),el}=uniqueReacs;
joiningPaths{pairs(k,1),pairs(k,2),el}=[a{:}];
sub2subNumP(pairs(k,1),pairs(k,2),el)=size([a{:}],2);
sub2subNumR(pairs(k,1),pairs(k,2),el)=length(uniqueReacs);
ModelReacs{el}=unique( [ModelReacs{el}; uniqueReacs(:) ],'stable');
%collect metabolites
m=L_DirAdjMatMets(source_mets{k}, target_mets{k},d);
m(emptyCells) = [];
uniqueMets=unique([m{:}],'stable');
joiningMets{pairs(k,1),pairs(k,2),el}=uniqueMets;
ModelMets{el}=unique( [ModelMets{el}; uniqueMets(:) ],'stable');
end
if shortestLsub2subNumP(pairs(k,2),pairs(k,1)) > 0
%other direction
d=el-1+shortestLsub2subNumP(pairs(k,2),pairs(k,1)); %switch!
a=L_DirAdjMatC(target_mets{k},source_mets{k},d); %switch!
emptyCells = cellfun('isempty',a);
a(emptyCells) = [];
uniqueReacs=unique([a{:}],'stable');
joiningReacs{pairs(k,2),pairs(k,1),el}=uniqueReacs; %switch!
joiningPaths{pairs(k,2),pairs(k,1),el}=[a{:}]; %switch!
sub2subNumP(pairs(k,2),pairs(k,1),el)=size([a{:}],2);
sub2subNumR(pairs(k,2),pairs(k,1),el)=length(uniqueReacs);
ModelReacs{el}=unique( [ModelReacs{el}; uniqueReacs(:) ],'stable');
%collect metabolites
m=L_DirAdjMatMets(target_mets{k}, source_mets{k},d);
m(emptyCells) = [];
uniqueMets=unique([m{:}],'stable');
joiningMets{pairs(k,2),pairs(k,1),el}=uniqueMets;
ModelMets{el}=unique( [ModelMets{el}; uniqueMets(:) ],'stable');
end
end
end
selectedPaths.sub2subNumP = sub2subNumP;
selectedPaths.sub2subNumR = sub2subNumR;
selectedPaths.joiningReacs = joiningReacs;
selectedPaths.joiningPaths = joiningPaths;
selectedPaths.joiningMets = joiningMets;
selectedPaths.ModelReacs = ModelReacs;
selectedPaths.ModelMets = ModelMets;
clear a emptyCells d i