enz_cent_RCM

 Builds Binary Stoichiometric Matrix and Enzyme-Enzyme Networks Removing Currency Metabolites
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 The function reads a Metabolic Network SBML file and builds 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(fileName)

INPUTS
 fileName                        The metabolic Network in the SBML format
 
OUTPUTS
 *_Removed_Mets_RCM.dat          file contains removed metabolits from the original model
 *_Enzyme_Cent_RCM.dat           Undirected-Enzyme-Enzyme Network (comma separated Format)
 *_Enzyme_Cent_RCM_Cyt.dat       Undirected-Enzyme-Enzyme Network - Cytoscape Compatible
 
 Yazdan Asgari 12/07/2012        http://lbb.ut.ac.ir
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [Output] = enz_cent_RCM(fileName)
0002 % Builds Binary Stoichiometric Matrix and Enzyme-Enzyme Networks Removing Currency Metabolites
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 % The function reads a Metabolic Network SBML file and builds 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(fileName)
0009 %
0010 %INPUTS
0011 % fileName                        The metabolic Network in the SBML format
0012 %
0013 %OUTPUTS
0014 % *_Removed_Mets_RCM.dat          file contains removed metabolits from the original model
0015 % *_Enzyme_Cent_RCM.dat           Undirected-Enzyme-Enzyme Network (comma separated Format)
0016 % *_Enzyme_Cent_RCM_Cyt.dat       Undirected-Enzyme-Enzyme Network - Cytoscape Compatible
0017 %
0018 % Yazdan Asgari 12/07/2012        http://lbb.ut.ac.ir
0019 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0020 
0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0022 % check validity of input file format
0023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0024 check=regexp(fileName,'.xml');
0025 assert(~isempty(check),'The SBML fileName must contain .xml at its end')
0026 
0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028 % start time evaluation of program
0029 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0030 tic;
0031 
0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0033 % reading the SBML file using COBRA Toolbox Command, and sets size of the S matrix
0034 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0035 model=readCbModel(fileName);
0036 [m,n]=size(model.S);
0037 
0038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0039 % calculate summation of each columns (i.e. How many metabolites each enzyme correlates)
0040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0041 S_bin=zeros(size(model.S));
0042 S_bin(find(model.S))=1;
0043 CB=sum(S_bin,1);
0044 A=zeros(m,n);
0045 B=zeros(m,1);
0046 N3=zeros(m,1);
0047 
0048 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0049 % for each binary S-matrix element, subtracts its value from the column summation and put the result in the A matrix.
0050 % A(q) means the metabolite q connects to how many other metabolites through the enzyme i.
0051 % for each row, sums the binary S-matrix over all columns.
0052 % B(q) means how many enzymes the metabolite q correlates.
0053 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0054 for q=1:m
0055     for i=1:n
0056         if S_bin(q,i)~=0
0057             A(q,i)=CB(1,i)-S_bin(q,i);
0058         end
0059         B(q,1)=B(q,1)+S_bin(q,i);
0060     end
0061 end
0062 
0063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0064 % Assumption: Generally, every metabolite is connected to the other one through a specific enzyme.
0065 % If a metabolite connects to more than one metabolite through an enzyme, this will be considered as a suspicious case.
0066 % Therefore, every N3(q) value equal to 3 will be marked for further analysis.
0067 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0068 for q=1:m
0069     for i=1:n
0070         if A(q,i)==3
0071             N3(q,1)=N3(q,1)+1;
0072         end
0073     end
0074 end
0075 
0076 s=0;
0077 for q=1:m
0078     if N3(q,1)~=0
0079         s=1;
0080     end
0081 end
0082 
0083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0084 % building the output file name for writing removed metabolites
0085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0086 outname1=strrep(fileName,'.xml','_Removed_Mets_RCM.dat')
0087 fout1 = fopen(outname1, 'w+');
0088 fprintf(fout1, 'Metabolite\t\tMetabolite Name\t\tMax1\t\tMax2\n');
0089 fprintf(fout1, '----------------------------------------------\n');
0090 
0091 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0092 % If there is any value for N3 array, the RCM algorithm will be done.
0093 % This algorithm will be deleted the most probable metabolite among all (i.e. the one with the maximum value of N3 and C)
0094 % The selected metabolite will be deleted in the binary S-Matrix, and the "WHILE LOOP" repeated.
0095 % The algorithm is ended if there is not any suspicious case.
0096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0097 while s==1
0098     C=zeros(m,1);
0099     max1=max(N3,[],1);
0100     for q=1:m
0101         if N3(q,1)==max1
0102             C(q,1)=B(q,1);
0103         else
0104             C(q,1)=0;
0105         end
0106     end
0107     max2=max(C,[],1);
0108     for q=1:m
0109         if ( (N3(q,1)==max1) && (C(q,1)==max2) )
0110             for i=1:n
0111                 S_bin(q,i)=0;
0112             end
0113             fprintf(fout1,'%s\t\t%s\t\t%d\t\t%d\n',model.mets{q},model.metNames{q},max1,max2);
0114         end
0115     end
0116     
0117     CB=sum(S_bin,1);
0118     A=zeros(m,n);
0119     B=zeros(m,1);
0120     N3=zeros(m,1);
0121     for q=1:m
0122         for i=1:n
0123             if S_bin(q,i)~=0
0124                 A(q,i)=CB(1,i)-S_bin(q,i);
0125             end
0126             B(q,1)=B(q,1)+S_bin(q,i);
0127         end
0128     end
0129     for q=1:m
0130         for i=1:n
0131             if A(q,i)==3
0132                 N3(q,1)=N3(q,1)+1;
0133             end
0134         end
0135     end
0136     s=0;
0137     for q=1:m
0138         if N3(q,1)~=0
0139             s=1;
0140         end
0141     end
0142 end
0143 
0144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0145 % construction of Undirected-Enzyme-Enzyme Network based on the new binary S-matrix(comma separated Format)
0146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0147 Aenz=S_bin'*S_bin;
0148 outname2=strrep(fileName,'.xml','_Enzyme_Cent_RCM.dat')
0149 dlmwrite(outname2,full(Aenz));
0150 
0151 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0152 % re-format of Undirected-Enzyme-Enzyme Network it to a Cytoscape-compatible file.
0153 % One could import the file using "File/Import/Network from Table(Text/MS Excel)..."
0154 % Select "first column" as "Source Interaction" and "second column" as "Target Interaction"
0155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0156 [m,n]=size(Aenz);
0157 outname3=strrep(fileName,'.xml','_Enzyme_Cent_RCM_Cyt.dat')
0158 fout2 = fopen(outname3, 'w+');
0159 for row=1:m
0160     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0161     % because cell(i,j)=cell(j,i) we must delete duplicate entries by putting
0162     % col=row:n in the second if command. since we must ignor diagonal elements,
0163     % the counter will be col=row+1:n
0164     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0165     for col=row+1:n
0166         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0167         % edge are those which includes number not equal to zero
0168         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0169          if Aenz(row,col)~=0
0170             fprintf(fout2, '%s\t%s\t%d\n',model.rxns{row},model.rxns{col},Aenz(row,col));
0171         end
0172     end
0173 end
0174 fclose(fout1);
0175 fclose(fout2);
0176 
0177 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0178 % End of time evaluation of program
0179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0180 toc;
0181