plotSampleHist

PURPOSE ^

plotSampleHist Compare flux histograms for one or more samples

SYNOPSIS ^

function plotSampleHist(rxnNames,samples,models,nBins,perScreen)

DESCRIPTION ^

plotSampleHist Compare flux histograms for one or more samples
for one or more reactions

 plotSampleHist(rxnNames,samples,models,nBins,perScreen)

INPUTS
 rxnNames  Cell array of reaction abbreviations
 samples   Cell array containing samples
 models    Cell array containing model structures or common model
 structure

OPTIONAL INPUTS
 nBins     Number of bins to be used
 perScreen Number of reactions to show per screen.  Either a number or [nY, nX] vector.
     (press Enter to advance screens)

CONTROLS
 To advance to next screen hit enter/return or type f and hit enter/return
 To rewind to previous screen type r or b and hit enter/return
 To quit script type q and hit enter/return

EXAMPLE
 sampleStructOut1 = gpSampler(model1, 2150);
 sampleStructOut2 = gpSampler(model2, 2150);
 Plot for model 1
 plotSampleHist(model1.rxns,{samplePoints1},{model1})

 Plot reactions reactions in model 1 that also exist in model 2 using 10
 bins and plotting 12 reactions per screen.
 plotSampleHist(model1.rxns,{samplePoints1,samplePoints2},{model1,model2},10,12)

 Markus Herrgard 7/17/06

 Added ability to move in reverse direction. Richard Que (2/05/10)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function plotSampleHist(rxnNames,samples,models,nBins,perScreen)
0002 %plotSampleHist Compare flux histograms for one or more samples
0003 %for one or more reactions
0004 %
0005 % plotSampleHist(rxnNames,samples,models,nBins,perScreen)
0006 %
0007 %INPUTS
0008 % rxnNames  Cell array of reaction abbreviations
0009 % samples   Cell array containing samples
0010 % models    Cell array containing model structures or common model
0011 % structure
0012 %
0013 %OPTIONAL INPUTS
0014 % nBins     Number of bins to be used
0015 % perScreen Number of reactions to show per screen.  Either a number or [nY, nX] vector.
0016 %     (press Enter to advance screens)
0017 %
0018 %CONTROLS
0019 % To advance to next screen hit enter/return or type f and hit enter/return
0020 % To rewind to previous screen type r or b and hit enter/return
0021 % To quit script type q and hit enter/return
0022 %
0023 %EXAMPLE
0024 % sampleStructOut1 = gpSampler(model1, 2150);
0025 % sampleStructOut2 = gpSampler(model2, 2150);
0026 % Plot for model 1
0027 % plotSampleHist(model1.rxns,{samplePoints1},{model1})
0028 %
0029 % Plot reactions reactions in model 1 that also exist in model 2 using 10
0030 % bins and plotting 12 reactions per screen.
0031 % plotSampleHist(model1.rxns,{samplePoints1,samplePoints2},{model1,model2},10,12)
0032 %
0033 % Markus Herrgard 7/17/06
0034 %
0035 % Added ability to move in reverse direction. Richard Que (2/05/10)
0036 
0037 if (nargin < 5)
0038     perScreen = 1e5;
0039 end
0040 
0041 if (nargin < 4 || isempty(nBins))
0042     [tmp,nSamples] = size(samples{1});
0043     nBins = round(nSamples/25);
0044 end
0045 
0046 if (~iscell(samples))
0047     singleDataFlag = true;
0048     samplesTmp = samples;
0049     clear samples;
0050     samples{1} = samplesTmp;
0051 end
0052 
0053 if (~iscell(models))
0054     commonModelFlag = true;
0055 else
0056     commonModelFlag = false;
0057 end
0058 
0059 if (~iscell(rxnNames))
0060     rxnNameList{1} = rxnNames;
0061 else
0062     rxnNameList = rxnNames;
0063 end
0064 
0065 nRxns = length(rxnNameList);
0066 
0067 for i = 1:nRxns
0068     if (commonModelFlag)
0069         keepRxn(i) = ismember(rxnNameList{i},models.rxns);
0070     else
0071         for j = 1:length(models)
0072             isIn(j) = ismember(rxnNameList{i},models{j}.rxns);
0073         end
0074         if (all(isIn))
0075             keepRxn(i) = true;
0076         else
0077             keepRxn(i) = false;
0078         end
0079     end
0080 end
0081 
0082 rxnNameList = rxnNameList(keepRxn');
0083 nRxns = length(rxnNameList);
0084 if length(perScreen) ==2
0085     nX = perScreen(2);
0086     nY = perScreen(1);
0087     perScreen = nX*nY;
0088 else
0089     nX = ceil(sqrt(min(nRxns, perScreen)));
0090     nY = ceil(min(nRxns, perScreen)/nX);
0091 end
0092 
0093 j=1;
0094 flagQuit = false;
0095 fig=figure;
0096 while ~flagQuit
0097     clear counts;
0098     currLB = 1e6;
0099     currUB = -1e6;
0100     for i = 1:length(samples)
0101         if (commonModelFlag)
0102             id = findRxnIDs(models,rxnNameList{j});
0103         else
0104             id = findRxnIDs(models{i},rxnNameList{j});
0105         end
0106         if (isempty(id))
0107             if (commonModelFlag)
0108                 id = findRxnIDs(models,[rxnNameList{j} '_r']);
0109             else
0110                 id = findRxnIDs(models{i},[rxnNameList{j} '_r']);
0111             end
0112             if (isempty(id))
0113                 warning('Reaction does not exist');
0114             end
0115         end
0116         currLB = min(currLB,min(samples{i}(id,:)'));
0117         currUB = max(currUB,max(samples{i}(id,:)'));
0118     end
0119     
0120     if currUB-currLB < 8*eps*abs(mean([currUB, currLB])) %rescale currUB, LB if they are the same # within numerical precision (8*eps)
0121         av = abs(mean([currUB, currLB]));
0122         currUB = currUB+8*eps*av;
0123         currLB = currLB-8*eps*av;
0124     end
0125     bins = linspace(currLB,currUB,nBins);
0126 
0127     for i = 1:length(samples)
0128         sampleSign = 1;
0129         if (commonModelFlag)
0130             id = findRxnIDs(models,rxnNameList{j});
0131         else
0132             id = findRxnIDs(models{i},rxnNameList{j});
0133         end
0134         if (isempty(id))
0135             if (commonModelFlag)
0136                 id = findRxnIDs(models,[rxnNameList{j} '_r']);
0137             else
0138                 id = findRxnIDs(models{i},[rxnNameList{j} '_r']);
0139             end
0140             sampleSign = -1;
0141         end
0142         n = hist(sampleSign*samples{i}(id,:),bins);
0143         counts(:,i) = smooth(bins,n');
0144     end
0145     freq = counts./repmat(sum(counts),size(counts,1),1);
0146     
0147     subplot(nY,nX, mod(j-1, perScreen)+1);
0148     plot(bins,freq,'-');
0149     axis([currLB-.0001 currUB+.0001 0 max(max(freq))]);
0150     h = text((currUB+currLB)/2,max(max(freq))*1.1,strrep(rxnNameList{j},'_','-'));
0151     set(h,'HorizontalAlignment','center');
0152     
0153     user_input = 'f';
0154     if j>=nRxns
0155         if nRxns<=perScreen, return; end
0156         user_input = input('End of sampples; reverse (r) or quit (q): ','s');
0157         while isempty(user_input)|| ~(ismember(user_input(1),{'b','r','q','e'}))
0158             user_input = input('End of sampples; reverse (r) or quit (q)): ','s');
0159         end
0160         clear fig
0161     elseif mod(j, perScreen) ==0
0162         user_input = input('Move forward (f), reverse (r) or quit (q): ','s');
0163         if isempty(user_input), user_input='f'; end
0164         clear fig
0165     end
0166     if user_input(1)=='r' || user_input(1)=='b'
0167         j=j-perScreen*2+1;
0168         if j<0, j=1; end
0169     elseif user_input(1)=='q'|| user_input(1)=='e'
0170         flagQuit=true;
0171     elseif user_input(1)=='f'
0172         j=j+1;
0173     end
0174 end
0175

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