identifyCorrelSets

PURPOSE ^

identifyCorrelSets Identify correlated reaction sets from sampling data

SYNOPSIS ^

function [setsSorted,setNoSorted,setSize] = identifyCorrelSets(model,samples,corrThr,R)

DESCRIPTION ^

identifyCorrelSets Identify correlated reaction sets from sampling data

 [sets,setNumber,setSize] =  identifyCorrelSets(model,samples,corrThr,R)

INPUTS
 model         COBRA model structure
 samples       Samples to be used to identify correlated sets

OPTIONAL INPUTS
 corrThr       Minimum correlation (R^2) threshold (Default = 1-1e-8)
 R             Correlation coefficient 

OUTPUTS
 sets          Sorted cell array of sets (largest first)
 setNumber     List of set numbers for each reaction in model (0 indicates
               that there is no set)
 setSize       List of set sizes

 Markus Herrgard 9/15/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [setsSorted,setNoSorted,setSize] = identifyCorrelSets(model,samples,corrThr,R)
0002 %identifyCorrelSets Identify correlated reaction sets from sampling data
0003 %
0004 % [sets,setNumber,setSize] =  identifyCorrelSets(model,samples,corrThr,R)
0005 %
0006 %INPUTS
0007 % model         COBRA model structure
0008 % samples       Samples to be used to identify correlated sets
0009 %
0010 %OPTIONAL INPUTS
0011 % corrThr       Minimum correlation (R^2) threshold (Default = 1-1e-8)
0012 % R             Correlation coefficient
0013 %
0014 %OUTPUTS
0015 % sets          Sorted cell array of sets (largest first)
0016 % setNumber     List of set numbers for each reaction in model (0 indicates
0017 %               that there is no set)
0018 % setSize       List of set sizes
0019 %
0020 % Markus Herrgard 9/15/06
0021 
0022 if (nargin < 3)
0023     corrThr = 1-1e-8;
0024 end
0025 
0026 nRxns = length(model.rxns);
0027 
0028 % Calculate correlation coefficients
0029 if (nargin < 4)
0030     R = corrcoef(samples');
0031     R = R - eye(nRxns);
0032 end
0033 
0034 % Define adjacency matrix
0035 adjMatrix = (abs(R) >= corrThr);
0036 
0037 % Only work with reactions that are correlated with others
0038 selCorrelRxns = any(adjMatrix)';
0039 rxnList = model.rxns(selCorrelRxns);
0040 adjMatrix = adjMatrix(selCorrelRxns,selCorrelRxns);
0041 
0042 % Construct set number index
0043 hasSet = false(size(rxnList));
0044 currSetNo = 0;
0045 setNoTmp = zeros(size(rxnList));
0046 for i = 1:length(rxnList)
0047     if (~hasSet(i))
0048         currSetNo = currSetNo+1;
0049         setMembers = find(adjMatrix(i,:));
0050         hasSet(setMembers) = true;
0051         hasSet(i) = true;
0052         setNoTmp(setMembers) = currSetNo;
0053         setNoTmp(i) = currSetNo;
0054     end
0055 end
0056 setNo = zeros(size(model.rxns));
0057 [tmp,index1,index2] = intersect(model.rxns,rxnList);
0058 setNo(index1) = setNoTmp(index2);  
0059 
0060 % Construct list of sets
0061 for i = 1:max(setNo)
0062     sets{i}.set = find(setNo == i);
0063     sets{i}.names = model.rxns(sets{i}.set);
0064     setSize(i) = length(sets{i}.set);
0065 end
0066 
0067 % Sort everything
0068 [setSize,sortInd] = sort(setSize');
0069 sortInd = flipud(sortInd);
0070 setsSorted = sets(sortInd);
0071 setNoSorted = zeros(size(setNo));
0072 for i = 1:length(sortInd)
0073     setNoSorted(setNo == sortInd(i)) = i;
0074 end
0075 setSize = flipud(setSize);
0076 setsSorted = setsSorted';

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