reporterMets

PURPOSE ^

reporterMets Implements the reporter metabolites algorithm by Patil &

SYNOPSIS ^

function [normScore,nRxnsMet,nRxnsMetUni,rawScore] = reporterMets(model,data,nRand,pValFlag,nLayers,metric,dataRxns,inclExchFlag)

DESCRIPTION ^

reporterMets Implements the reporter metabolites algorithm by Patil &
Nielsen

 [normScore,nRxnsMet,nRxnsMetUni,rawScore] = reporterMets(model,data,nRand,pValFlag,nLayers,metric,dataRxns)

 model         Metabolic network reconstruction structure
 data          Data matrix/vector
 nRand         Number of randomizations
 pValFlag      The data are p-values and should be converted to z-scores
 nLayers       Number of reaction layers around each metabolite considered (default = 1)
 metric        Metric used to evaluate score
               ('default','mean','median','std','count')
 dataRxns      Reaction list for the data file (if different from the
 model reactions)
 inclExchFlag

 normScore     Normalized scores for each metabolite
 nRxnsMet      Number of reactions connected to each metabolite
 nRxnsMetUni
 rawScore      Raw unnormalized scores

 Markus Herrgard 7/20/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [normScore,nRxnsMet,nRxnsMetUni,rawScore] = reporterMets(model,data,nRand,pValFlag,nLayers,metric,dataRxns,inclExchFlag)
0002 %reporterMets Implements the reporter metabolites algorithm by Patil &
0003 %Nielsen
0004 %
0005 % [normScore,nRxnsMet,nRxnsMetUni,rawScore] = reporterMets(model,data,nRand,pValFlag,nLayers,metric,dataRxns)
0006 %
0007 % model         Metabolic network reconstruction structure
0008 % data          Data matrix/vector
0009 % nRand         Number of randomizations
0010 % pValFlag      The data are p-values and should be converted to z-scores
0011 % nLayers       Number of reaction layers around each metabolite considered (default = 1)
0012 % metric        Metric used to evaluate score
0013 %               ('default','mean','median','std','count')
0014 % dataRxns      Reaction list for the data file (if different from the
0015 % model reactions)
0016 % inclExchFlag
0017 %
0018 % normScore     Normalized scores for each metabolite
0019 % nRxnsMet      Number of reactions connected to each metabolite
0020 % nRxnsMetUni
0021 % rawScore      Raw unnormalized scores
0022 %
0023 % Markus Herrgard 7/20/06
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     %error('Not implemented yet')
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 % Handle case where more than one rxn is associated with a
0067 % data value
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 % Compute raw scores
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 % Do randomization
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                 % Sample with replacement
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

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