-
Notifications
You must be signed in to change notification settings - Fork 0
/
genCytoscapeNetwork.m
132 lines (115 loc) · 3.73 KB
/
genCytoscapeNetwork.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
function param = genCytoscapeNetwork(param, v)
%% generate a .sif file for Cytoscape for pathway name(s)
%
% Author: Jiun Yen
% Date: 2017.4.14
% Version: 2017.4.22
%
% param structure
% .m COBRA-format GEM model
% .pathways structure
% .rids
% .names
% .values
% .counts
% .metFreq frequency of metabolite used in reactions
% .rids indices of reactions with fluxes
% .pName path name
% .fName file name
% .maxMetFreq upper limit for metabolite frequency in S matrix
% .excludeMets metabolites excluded from network
% .includeMets metabolites to include if they are in the pathways
%% find pathway rids
p_rids = param.pathways.rids(ismember(param.pathways.values,find(ismember(param.pathways.names, param.pName))));
%% remove inactive rids
try
rids = intersect(param.rids, p_rids);
catch
error('No rids');
end
%% Update param
S = full(param.m.S);
param.networkRids = rids;
param.networkRxns = cell(length(rids), 1);
param.networkInteractions = cell(2*length(rids), 1);
allmets = logical(sum(abs(S(:,rids)),2));
excludemets = ismember(param.m.metNames, param.excludeMets);
includemets = ismember(param.m.metNames, param.includeMets);
hifreqmets = param.metFreq > param.maxMetFreq;
excludemets = (excludemets | (allmets & hifreqmets)) & (allmets & ~includemets);
param.allMets = param.m.metNames(allmets);
param.excludeMets = param.m.metNames(excludemets);
param.networkMets = param.m.metNames(allmets & ~excludemets);
param.rids = rids;
%% construct network
% remove mets that are exclude from S-matrix
S(excludemets, :) = [];
param.m.metNames(excludemets) = [];
% LHS S-matrix (reactants)
SL = S < 0;
% RHS S-matrix (products)
SR = S > 0;
% write to file
f = fopen([param.fName '_network.sif'], 'w');
for i = 1:length(rids)
rxnName = param.m.rxnNames{rids(i)};
if isempty(rxnName)
rxnName = param.m.rxns{rids(i)};
end
rxnName = editName(rxnName);
metL = param.m.metNames(SL(:,rids(i)));
metR = param.m.metNames(SR(:,rids(i)));
tmp = sprintf('r%u',rids(i));
param.networkRxns(i) = {[tmp ':' rxnName]};
param.networkInteractions(2*i-1) = {[tmp 'L']};
param.networkInteractions(2*i) = {[tmp 'R']};
if ~isempty(metL) && ~isempty(metR)
for a = 1:length(metL)
fprintf(f, [param.networkRxns{i} '\t' param.networkInteractions{2*i-1} '\t' metL{a} '\n']);
end
for b = 1:length(metR)
fprintf(f, [param.networkRxns{i} '\t' param.networkInteractions{2*i} '\t' metR{b} '\n']);
end
end
end
fclose(f);
%% construct node properties
f = fopen([param.fName '_node.txt'], 'w');
fprintf(f, 'name\ttype\n');
for i = 1:length(param.networkMets)
fprintf(f, [param.networkMets{i} '\tmet\n']);
end
for i = 1:length(param.networkRxns)
fprintf(f, [param.networkRxns{i} '\trxn\n']);
end
fclose(f);
%% construct edge properties
% check if necessary
if nargin < 2 || isempty(v)
return;
end
% build directionality vector
dL = (v(rids) < 0) + 0;
dR = (v(rids) > 0) + 0;
v = abs(v(rids));
toDraw = (v ~= 0) + 0;
% build edge properties
f = fopen([param.fName '_edge.txt'], 'w');
fprintf(f, 'interaction\tdirection\tdraw\tweight\n');
for i = 1:length(rids)
fprintf(f, [param.networkInteractions{2*i-1} '\t%u\t%u\t%4.4f\n'], [dL(i) toDraw(i) v(i)]);
fprintf(f, [param.networkInteractions{2*i} '\t%u\t%u\t%4.4f\n'], [dR(i) toDraw(i) v(i)]);
end
fclose(f);
% save param
save([param.fName '_param']);
function name = editName(name)
%% trim rxnName to critical component
% also because Cytoscape is awful with pre-treating special characters
tmp = regexp(name,'[ -}]*','match');
tmp = strsplit(tmp{1}, ' / ');
tmp = strrep(tmp{1}, ', putative', '');
tmp = strsplit(tmp, ';');
tmp = regexprep(tmp{1},'[.*]','');
tmp = strrep(tmp,' ','');
name = strtrim(tmp);