0001 function validJunctionMets = findMetabolicJunctions(model,nRxnsJnc)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
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
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