findNeighborRxns Identifies the reactions and the corresponding genes that are adjacent (having a common metabolite) to a reaction of interest. Useful for characterizing the network around an orphan reaction. [neighborRxns,neighborGenes,mets] = findNeighborRxns(model,rxn) INPUTS model COBRA model structure rxn the target reaction (only 1 for now) OUTPUTS neighborRxns the neighboring rxns in the network, (having common metabolites) neighborGenes the genes associated with the neighbor rxns mets the metabolites in the target reaction Jeff Orth 10/11/09
0001 function [neighborRxns,neighborGenes,mets] = findNeighborRxns(model,rxn) 0002 %findNeighborRxns Identifies the reactions and the corresponding genes 0003 %that are adjacent (having a common metabolite) to a reaction of interest. 0004 %Useful for characterizing the network around an orphan reaction. 0005 % 0006 % [neighborRxns,neighborGenes,mets] = findNeighborRxns(model,rxn) 0007 % 0008 %INPUTS 0009 % model COBRA model structure 0010 % rxn the target reaction (only 1 for now) 0011 % 0012 %OUTPUTS 0013 % neighborRxns the neighboring rxns in the network, (having common 0014 % metabolites) 0015 % neighborGenes the genes associated with the neighbor rxns 0016 % mets the metabolites in the target reaction 0017 % 0018 % Jeff Orth 0019 % 10/11/09 0020 0021 %get the metabolites in the rxn 0022 metIndex = find(model.S(:,findRxnIDs(model,rxn))); 0023 0024 % exclude common mets (atp, adp, h, h2o, pi) ** make this an input option 0025 iAtpC = findMetIDs(model,'atp[c]'); 0026 iAtpP = findMetIDs(model,'atp[p]'); 0027 iAdpC = findMetIDs(model,'adp[c]'); 0028 iAdpP = findMetIDs(model,'adp[p]'); 0029 iHC = findMetIDs(model,'h[c]'); 0030 iHP = findMetIDs(model,'h[p]'); 0031 iH2oC = findMetIDs(model,'h2o[c]'); 0032 iH2oP = findMetIDs(model,'h2o[p]'); 0033 iPiC = findMetIDs(model,'pi[c]'); 0034 iPiP = findMetIDs(model,'pi[p]'); 0035 metIndex = setdiff(metIndex,[iAtpC,iAtpP,iAdpC,iAdpP,iHC,iHP,iH2oC,iH2oP,iPiC,iPiP]); 0036 0037 %get the rxns for each met 0038 nRxnIndexs = {}; 0039 for i = 1:length(metIndex) 0040 nRxnIndexs{i} = find(model.S(metIndex(i),:)); 0041 end 0042 0043 % remove target rxn from list 0044 for i = 1:length(metIndex); 0045 nRxnIndexs{i} = setdiff(nRxnIndexs{i},findRxnIDs(model,rxn)); 0046 end 0047 0048 neighborRxns = {}; 0049 for i = 1:length(metIndex) 0050 neighborRxns{i} = model.rxns(nRxnIndexs{i}); 0051 end 0052 0053 %get genes for each rxn 0054 neighborGenes = {}; 0055 for i = 1:length(metIndex) 0056 neighborGenes{i} = model.grRules(nRxnIndexs{i}); 0057 end 0058 0059 mets = model.mets(metIndex); 0060 0061