0001 function [normScore,nRxnsMet,nRxnsMetUni,rawScore] = reporterMets(model,data,nRand,pValFlag,nLayers,metric,dataRxns,inclExchFlag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025 if (nargin < 3)
0026 nRand = 1000;
0027 end
0028
0029 if (nargin < 4)
0030 pValFlag = false;
0031 end
0032
0033 if (nargin < 5)
0034 nLayers = 1;
0035 end
0036
0037 if (nargin < 6)
0038 metric = 'default';
0039 else
0040 if (isempty(metric))
0041 metric = 'default';
0042 end
0043 end
0044
0045 if (nargin < 7)
0046 dataRxns = model.rxns;
0047 else
0048 if isempty(dataRxns)
0049 dataRxns = model.rxns;
0050 end
0051 end
0052
0053 if (nargin < 8)
0054 inclExchFlag = false;
0055 end
0056
0057 if (pValFlag)
0058
0059 minP = min(min(data(data ~= 0)));
0060 data(data == 0) = minP;
0061 data = -norminv(data,0,1);;
0062 end
0063
0064 [nRxnsTot,nData] = size(data);
0065
0066
0067
0068 if (iscell(dataRxns{1}))
0069 for j = 1:length(dataRxns)
0070 dataRxnTmp = dataRxns{j};
0071 selModelRxn = ismember(model.rxns,dataRxnTmp);
0072 dataModelRxnMap(j,:) = selModelRxn';
0073 end
0074 dataModelRxnMap = sparse(dataModelRxnMap);
0075 end
0076
0077 excInd = find(findExcRxns(model),1);
0078
0079
0080 for i = 1:length(model.mets)
0081 rxnInd = find(model.S(i,:) ~= 0)';
0082 for j = 2:nLayers
0083 metInd = find(any(model.S(:,rxnInd) ~= 0,2));
0084 rxnInd = union(rxnInd,find(any(model.S(metInd,:) ~= 0,1)));
0085 end
0086 rxnInd = setdiff(rxnInd,excInd);
0087 if (nargin < 7)
0088 dataRxnInd = false(length(dataRxns),1);
0089 dataRxnInd(rxnInd) = true;
0090 nRxnsMetUni(i) = sum(rxnInd);
0091 else
0092 if (~iscell(dataRxns{1}))
0093 dataRxnInd = ismember(dataRxns,model.rxns(rxnInd));
0094 nRxnsMetUni(i) = length(unique(dataRxns(dataRxnInd)));
0095 else
0096 dataRxnInd = any(full(dataModelRxnMap(:,rxnInd)),2);
0097 nRxnsMetUni(i) = sum(dataRxnInd);
0098 end
0099 end
0100 nRxnsMet(i) = sum(dataRxnInd);
0101
0102 if (nRxnsMet(i) == 0)
0103 rawScore(i,:) = zeros(1,nData);
0104 else
0105 if (sum(dataRxnInd) == 1)
0106 switch metric
0107 case {'default','mean','median'}
0108 rawScore(i,:) = data(dataRxnInd,:);
0109 case 'std'
0110 rawScore(i,:) = zeros(1,nData);
0111 case 'count'
0112 rawScore(i,:) = data(dataRxnInd,:) ~= 0;
0113 end
0114 else
0115 switch metric
0116 case 'default'
0117 rawScore(i,:) = nansum(data(dataRxnInd,:))/sqrt(nRxnsMet(i));
0118 case 'mean'
0119 rawScore(i,:) = nanmean(data(dataRxnInd,:));
0120 case 'median'
0121 rawScore(i,:) = nanmedian(data(dataRxnInd,:));
0122 case 'std'
0123 rawScore(i,:) = nanstd(data(dataRxnInd,:));
0124 case 'count'
0125 rawScore(i,:) = nansum(data(dataRxnInd,:) ~= 0);
0126 end
0127 end
0128 end
0129 end
0130 nRxnsMet = nRxnsMet';
0131 nRxnsMetUni = nRxnsMetUni';
0132
0133 nRxnsUni = unique(nRxnsMet);
0134
0135
0136 normScore = zeros(length(rawScore),nData);
0137 if (~isempty(nRand))
0138 for i = 1:length(nRxnsUni)
0139 nRxns = nRxnsUni(i);
0140 if (nRxns > 0)
0141 randScore = [];
0142 for j = 1:nRand
0143
0144 randInd = randint(nRxns,1,[1 nRxnsTot]);
0145 if (length(randInd) == 1)
0146 switch metric
0147 case {'default','std','mean','median'}
0148 randScore(j,:) = data(randInd,:);
0149 case 'count'
0150 randScore(j,:) = data(randInd,:) ~= 0;
0151 end
0152 else
0153 switch metric
0154 case 'default'
0155 randScore(j,:) = nansum(data(randInd,:))/sqrt(nRxns);
0156 case 'mean'
0157 randScore(j,:) = nanmean(data(randInd,:));
0158 case 'median'
0159 randScore(j,:) = nanmedian(data(randInd,:));
0160 case 'std'
0161 randScore(j,:) = nanstd(data(randInd,:));
0162 case 'count'
0163 randScore(j,:) = nansum(data(randInd,:) ~= 0);
0164 end
0165 end
0166 end
0167 randMean = mean(randScore);
0168 randStd = std(randScore);
0169 metInd = find(nRxnsMet == nRxns);
0170 normScore(metInd,:) = (rawScore(metInd,:)-repmat(randMean,length(metInd),1))./repmat(randStd,length(metInd),1);
0171 end
0172 end
0173 else
0174 normScore = rawScore;
0175 end