findMetabolicJunctions

PURPOSE ^

findMetabolicJunctions Find metabolic branchpoints with different numbers

SYNOPSIS ^

function validJunctionMets = findMetabolicJunctions(model,nRxnsJnc)

DESCRIPTION ^

findMetabolicJunctions Find metabolic branchpoints with different numbers
of branches

 validJunctionMets = findMetabolicJunctions(model,nRxnsJnc)

INPUTS
 model                 COBRA model structure
 nRxnJnc               Number of reactions to be considered a junction

OUTPUT
 validJunctionMets     List of junction metabolites

 Markus Herrgard 12/14/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function validJunctionMets = findMetabolicJunctions(model,nRxnsJnc)
0002 %findMetabolicJunctions Find metabolic branchpoints with different numbers
0003 %of branches
0004 %
0005 % validJunctionMets = findMetabolicJunctions(model,nRxnsJnc)
0006 %
0007 %INPUTS
0008 % model                 COBRA model structure
0009 % nRxnJnc               Number of reactions to be considered a junction
0010 %
0011 %OUTPUT
0012 % validJunctionMets     List of junction metabolites
0013 %
0014 % Markus Herrgard 12/14/06
0015 
0016 if (isfield(model,'c'))
0017     selRxnsC = (model.c == 0);
0018 else
0019     selRxnsC = true(length(model.rxns),1);
0020 end
0021 
0022 [baseMetNames,compSymbols,uniqueMetNames,uniqueCompSymbols] = parseMetNames(model.mets);
0023 uniqueMetNames = uniqueMetNames';
0024 for i = 1:length(uniqueMetNames)
0025     sel = ismember(baseMetNames,uniqueMetNames{i});
0026     nRxnsMetUni(i) = sum(any(model.S(sel,selRxnsC) ~= 0,1));
0027 end
0028 nRxnsMetUni = full(nRxnsMetUni');
0029 junctionMets = uniqueMetNames(nRxnsMetUni >= nRxnsJnc);
0030 
0031 validJunctionMets = {};
0032 for i = 1:length(junctionMets)
0033     sel = ismember(baseMetNames,junctionMets{i});
0034     if (length(unique(compSymbols(sel))) == 1)
0035         selRxns = any(model.S(sel,:) ~= 0,1) & selRxnsC';
0036         thisRxns = model.rxns(selRxns);
0037         geneMap = model.rxnGeneMat(findRxnIDs(model,thisRxns),:);
0038         selNonZero = any(geneMap,2);
0039         if (sum(selNonZero) == nRxnsJnc & ...
0040             size(unique(geneMap(selNonZero,:),'rows'),1) == nRxnsJnc)
0041             validJunctionMets{end+1} = junctionMets{i};
0042             if (verbFlag)
0043                 fprintf('*** %s ***\n',junctionMets{i});
0044                 for j = 1:length(thisRxns)
0045                     %fprintf('%s\t',thisRxns{j});
0046                     geneInd = find(model.rxnGeneMat(findRxnIDs(model,thisRxns{j}),:));
0047                     if (~isempty(geneInd))
0048                         thisGenes = model.genes(geneInd);
0049                         for k = 1:length(thisGenes)
0050                             fprintf('%s ',thisGenes{k});
0051                         end
0052                     end
0053                     fprintf('\n');
0054                 end
0055             end
0056         end
0057     end
0058 end
0059 
0060 validJunctionMets = validJunctionMets';
0061 for i = 1:length(validJunctionMets)
0062     validJunctionMets{i} = [validJunctionMets{i} '(c)'];
0063 end

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003