Builds Directed Enzyme-Enzyme Networks Removing Currency Metabolites %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The function reads a Metabolic Network SBML file and builds Directed Enzyme-Enzyme Networks. The Remove Currency Metabolites (RCM) algorithm removes currency metabolites in the metabolic network automatically. Note: COBRA Toolbox must be installed in MATLAB before running this function [Output] = enz_cent_RCM_dir(fileName) INPUTS fileName The metabolic Network in the SBML format OUTPUTS *_Enzyme_Cent_RCM_Dir_Cyt.sif Directed-Enzyme-Enzyme Network - Cytoscape Compatible (.sif file) Yazdan Asgari 12/07/2012 http://lbb.ut.ac.ir %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 function [Output] = enz_cent_RCM_dir(fileName) 0002 % Builds Directed Enzyme-Enzyme Networks Removing Currency Metabolites 0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0004 % The function reads a Metabolic Network SBML file and builds Directed Enzyme-Enzyme Networks. 0005 % The Remove Currency Metabolites (RCM) algorithm removes currency metabolites in the metabolic network automatically. 0006 % Note: COBRA Toolbox must be installed in MATLAB before running this function 0007 % 0008 % [Output] = enz_cent_RCM_dir(fileName) 0009 % 0010 %INPUTS 0011 % fileName The metabolic Network in the SBML format 0012 % 0013 %OUTPUTS 0014 % *_Enzyme_Cent_RCM_Dir_Cyt.sif Directed-Enzyme-Enzyme Network - Cytoscape Compatible (.sif file) 0015 % 0016 % Yazdan Asgari 12/07/2012 http://lbb.ut.ac.ir 0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0018 0019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0020 % check validity of input file format 0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0022 check=regexp(fileName,'.xml'); 0023 assert(~isempty(check),'The SBML fileName must contain .xml at its end') 0024 0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0026 % start time evaluation of program 0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 tic; 0029 0030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0031 % reading the SBML file using COBRA Toolbox Command, and sets size of the S matrix 0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0033 model=readCbModel(fileName); 0034 [m,n]=size(model.S); 0035 0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0037 % calculate summation of each columns (i.e. How many metabolites each enzyme correlates) 0038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0039 S_bin=zeros(size(model.S)); 0040 S_bin(find(model.S))=1; 0041 CB=sum(S_bin,1); 0042 A=zeros(m,n); 0043 B=zeros(m,1); 0044 N3=zeros(m,1); 0045 0046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0047 % for each binary S-matrix element, subtracts its value from the column summation and put the result in the A matrix. 0048 % A(q) means the metabolite q connects to how many other metabolites through the enzyme i. 0049 % for each row, sums the binary S-matrix over all columns. 0050 % B(q) means how many enzymes the metabolite q correlates. 0051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0052 for q=1:m 0053 for i=1:n 0054 if S_bin(q,i)~=0 0055 A(q,i)=CB(1,i)-S_bin(q,i); 0056 end 0057 B(q,1)=B(q,1)+S_bin(q,i); 0058 end 0059 end 0060 0061 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0062 % Assumption: Generally, every metabolite is connected to the other one through a specific enzyme. 0063 % If a metabolite connects to more than one metabolite through an enzyme, this will be considered as a suspicious case. 0064 % Therefore, every N3(q) value equal to 3 will be marked for further analysis. 0065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0066 for q=1:m 0067 for i=1:n 0068 if A(q,i)==3 0069 N3(q,1)=N3(q,1)+1; 0070 end 0071 end 0072 end 0073 0074 s=0; 0075 for q=1:m 0076 if N3(q,1)~=0 0077 s=1; 0078 end 0079 end 0080 0081 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0082 % If there is any value for N3 array, the RCM algorithm will be done. 0083 % This algorithm will be deleted the most probable metabolite among all (i.e. the one with the maximum value of N3 and C) 0084 % The selected metabolite will be deleted in the binary S-Matrix, and the "WHILE LOOP" repeated. 0085 % The algorithm is ended if there is not any suspicious case. 0086 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0087 while s==1 0088 C=zeros(m,1); 0089 max1=max(N3,[],1); 0090 for q=1:m 0091 if N3(q,1)==max1 0092 C(q,1)=B(q,1); 0093 else 0094 C(q,1)=0; 0095 end 0096 end 0097 max2=max(C,[],1); 0098 for q=1:m 0099 if ( (N3(q,1)==max1) && (C(q,1)==max2) ) 0100 for i=1:n 0101 S_bin(q,i)=0; 0102 model.S(q,i)=0; 0103 end 0104 end 0105 end 0106 CB=sum(S_bin,1); 0107 A=zeros(m,n); 0108 B=zeros(m,1); 0109 N3=zeros(m,1); 0110 for q=1:m 0111 for i=1:n 0112 if S_bin(q,i)~=0 0113 A(q,i)=CB(1,i)-S_bin(q,i); 0114 end 0115 B(q,1)=B(q,1)+S_bin(q,i); 0116 end 0117 end 0118 for q=1:m 0119 for i=1:n 0120 if A(q,i)==3 0121 N3(q,1)=N3(q,1)+1; 0122 end 0123 end 0124 end 0125 s=0; 0126 for q=1:m 0127 if N3(q,1)~=0 0128 s=1; 0129 end 0130 end 0131 end 0132 0133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0134 % building the output file name (Cytoscape compatble file) 0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0136 outname=strrep(fileName,'.xml','_Enzyme_Cent_RCM_Dir_Cyt.sif') 0137 fout = fopen(outname, 'w+'); 0138 0139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0140 % construction of Directed-Enzyme-Enzyme Network based on the new binary S-matrix 0141 % finds non-zero elements of the new S-matrix (in order to make the algorithm faster), 0142 % parses through each row, and considers an edge for every unlike-signs, 0143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0144 for j=1:m 0145 indices=find(model.S(j,:)); 0146 [a,b]=size(indices); 0147 r=0; 0148 if b~=0 0149 r=1; 0150 end 0151 while r<b 0152 i=1; 0153 while i<(b-r+1) 0154 if model.S(j,indices(1,r))<0 && model.S(j,indices(1,r+i))>0 0155 fprintf(fout,'%s\t%s\t%s\n',model.rxns{indices(1,r)},'reaction-product',model.rxns{indices(1,r+i)}); 0156 elseif model.S(j,indices(1,r))>0 && model.S(j,indices(1,r+i))<0 0157 fprintf(fout,'%s\t%s\t%s\n',model.rxns{indices(1,r+i)},'reaction-product',model.rxns{indices(1,r)}); 0158 end 0159 i=i+1; 0160 end 0161 r=r+1; 0162 end 0163 end 0164 fclose(fout); 0165 0166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0167 % End of time evaluation of program 0168 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0169 toc;