/
addFlux.m
158 lines (142 loc) · 5.08 KB
/
addFlux.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
function [ref,normalizedFlux,newListRxn] = addFlux(model,FBA_result,parsed,listRxn)
% Changes the thickness attributes of the reaction links in a `CellDesigner`
% model structure and make them propotional to the flux values obtained
% from COBRA functions
%
% USAGE:
%
% [ref, normalizedFlux, newListRxn] = addFlux(model, FBA_result, parsed, listRxn)
%
% INPUTS:
%
% model: A COBRA model structure
% FBAresult: FBA results of a COBRA simulation by the
% `optimizeCbModel` function
% parsed: The CD model structure outputed by the `parseCD`
% function
% listRxn: A list of reaction IDs, after which the flux values
% are modified or added. P.s., the reaction IDs must be
% present in both the parsed CD model and the COBRA
% model structures. the reactions can be examined by
% `cmpR` function.
%
% OUTPUTS:
% ref: An updated parsed `CellDesigner` model that include
% information about the width of reaction links
% normalizedFlux: A list of normalised flux values generated based on
% `FBAresult`.
% newListRxn: New list of reaction IDs
%
% EXAMPLE:
%
% [parsed_1] = addFlux(recon2, fba_results, parsePD, listRxn);
%
% .. Author: - Longfei Mao Oct/2014
ref=parsed;
FBAsolution=FBA_result;
%%%% 1)
% if ~isempty(listRxn)
% for ll=1:length(listRxn);
% ind(ll)=find(strcmp(listRxn(ll),model.rxns(:))); % record every index number for each reaction present in the model.rxn.
% end
% else
% error('no reaction in the reaction list can be found');
% end
%%%% 2)
% for 1:length(ind(:))
%%%%%%%%%%%% obtain the reaction indecies in the model
if (iscell(listRxn))
[tmp,rxnID] = ismember(listRxn,model.rxns);
%%%% check if the reaction names of the listRxn appear in the list of
%%%% the rection name list of the model.
if rxnID==0;
[tmp,rxnID] = ismember(listRxn,model.rxnNames);
end
else
rxnID = find(strcmp(model.rxns,listRxn));
if (isempty(rxnID))
rxnID = 0;
end
if (length(rxnID) > 1)
rxnID = rxnID(1);
end
%%%% check if the reaction names of the listRxn appear in the list of
%%%% the rection name list of the model.
if rxnID==0
rxnID = find(strcmp(model.rxnNames,listRxn));
end
end
if rxnID==0
disp(rxnID)
error('cannot find the reaction IDs in the model');
end
flux(:,1)=FBAsolution.x(rxnID)
%%%%%% normalise the flux values.
absFlux=abs(flux);
rxnWidth=absFlux/max(absFlux);
rxnWidth(rxnWidth>=1)=8; % normalize the values in the descending order
rxnWidth(rxnWidth>0.8 & rxnWidth<1)=6;
rxnWidth(rxnWidth>0.5 & rxnWidth<=0.8)=5;
rxnWidth(rxnWidth>0.2 & rxnWidth<=0.5)=4;
rxnWidth(rxnWidth>1e-3 & rxnWidth<=0.2)=3;
rxnWidth(rxnWidth<1e-3)=0;
% normalized values
normalizedFlux=flux(:,1);
normalizedFlux(:,2)=rxnWidth(:,1);
%results=[];
[ID_row,ID_Column]=size(ref.r_info.ID);
for m=1:ID_row;
for n=1:ID_Column;
r=iscellstr(ref.r_info.ID(m,n));
if ~r;
%results(or,1)=m;
%results(or,2)=n;
ref.r_info.ID{m,n}=' '
end;
end;
end
for r=1: length(flux)
rxnList_width(r)=rxnWidth(r);
id=find(ismember(ref.r_info.ID,listRxn(r))) % a situation where the third column stores the reaction IDs.
if id
[m,n]=size(ref.r_info.ID);
if id>m*(n-1) % the third column of ref.r_info.ID contains reaction ID; for example: {'re5160','re8','DESAT16_2'}
newRxnName=ref.r_info.ID{id-2*m}; % newRxnName is defined to be the ID in the first column of ref.r_info.ID.
end
try
[rw,cw]=size(ref.(newRxnName).width)
catch
return
end
for ddr=1:rw
for ddc=1:cw
ref.(newRxnName).width{ddr,ddc}=rxnList_width(r)
fprintf('set %s ''s width to %d \n',newRxnName,rxnList_width(r));
end
end
newListRxn{r,1}=newRxnName;
else
% newRxnName=listRxn{r};
%%%%%%%% modifying the parsed CD model
% ref.(listRxn{r}) %
if ~isfield(ref,listRxn{r})
newRxnName=strcat('R_',listRxn{r});
if isempty(strfind(newRxnName,'(e)'))
newRxnName=strrep(newRxnName,'(e)','_e');
end
if ~isfield(ref,newRxnName)
disp(listRxn{r});
fprintf('error ! the listRxn{%d}',r);
r=r+1 % procede to the next reaction
else
[rw,cw]=size(ref.(newRxnName).width)
for ddr=1:rw
for ddc=1:cw
ref.(newRxnName).width{ddr,ddc}=rxnList_width(r)
fprintf('set %s ''s width to %d \n',newRxnName,rxnList_width(r));
end
end
end
end
end
end