Compare the multiple sets of samples xglc is optional, a random sugar distribution is calculated if empty expects samp1 and samp2 to have a field named points containing an array of sampled points expects model.rxns to contain a list of rxn names measuredMetabolites is an optional parameter fed to calcMDVfromSamp.m which only calculates the MDVs for the metabolites listed in this array. e.g. totalz is the sum of all zscores zscore is the calculated difference for each mdv element distributed across all the points mdv1,mdv2 each containing fields: - mdv - the calculated mdv distribution converted from the idv solved from each point contained in their respective samples sampX - names - the names of the metabolites - ave - the average of each mdv element across all of the points - stdev - the standard dev for each mdv element across all points Wing Choi 2/11/08
0001 function [totalz,zscore,mdvs] = compareMultSamp(xglc,model,samps,measuredMetabolites) 0002 0003 % Compare the multiple sets of samples 0004 % xglc is optional, a random sugar distribution is calculated if empty 0005 % expects samp1 and samp2 to have a field named points containing 0006 % an array of sampled points 0007 % expects model.rxns to contain a list of rxn names 0008 % measuredMetabolites is an optional parameter fed to calcMDVfromSamp.m 0009 % which only calculates the MDVs for the metabolites listed in this 0010 % array. e.g. 0011 % 0012 % totalz is the sum of all zscores 0013 % zscore is the calculated difference for each mdv element distributed 0014 % across all the points 0015 % mdv1,mdv2 each containing fields: 0016 % - mdv - the calculated mdv distribution converted from the idv 0017 % solved from each point contained in their respective samples sampX 0018 % - names - the names of the metabolites 0019 % - ave - the average of each mdv element across all of the points 0020 % - stdev - the standard dev for each mdv element across all points 0021 % Wing Choi 2/11/08 0022 0023 %glc = zeros(62,1); 0024 %glc = [.2 ;.8 ;glc]; 0025 0026 if (nargin < 3) 0027 disp '[totalz,zscore,mdvs] = compareMultSamp(xglc,model,samps,measuredMetabolites)'; 0028 return; 0029 end 0030 0031 if (nargin < 4) 0032 measuredMetabolites = []; 0033 end 0034 0035 mdvs.ave = []; 0036 mdvs.stdev = []; 0037 0038 if (isempty(xglc)) 0039 % random glucose 0040 xglc = rand(64,1); 0041 xglc = xglc/sum(xglc); 0042 xglc = idv2cdv(6)*xglc; 0043 end 0044 0045 % generate the translation index array 0046 % can shave time by not regenerating this array on every call. 0047 xltmdv = zeros(1,4096); 0048 for i = 1:4096 0049 xltmdv(i) = length(strrep(dec2base(i-1,2),'0','')); 0050 end 0051 0052 % calculate mdv for samp1 and samp2 0053 for i = 1:length(samps) 0054 samp.points = samps(i).points; 0055 mdv = calcMDVfromSamp(model,xglc,samp,measuredMetabolites); 0056 l = length(mdv.ave); 0057 mdvs.ave(1:l,i) = mdv.ave; 0058 mdvs.stdev(1:l,i) = mdv.stdev; 0059 %[mdv2] = calcMDVfromSamp(model,xglc,samp2,measuredMetabolites); 0060 %[totalz,zscore] = compareTwoMDVs(mdv1,mdv2); 0061 end 0062 0063 totalz = 0; 0064 zscore = 0; 0065 %mdvs = mdv; 0066 0067 return 0068 end 0069 0070 0071 %Here's what the function does. 0072 %Apply slvrXXfast to each point 0073 %for each field in the output, apply iso2mdv to get a much shorter vector. 0074 %store all the mdv's for each point and for each metabolite in both sets. 0075 % 0076 %Compute the mean and standard deviation of each mdv and then get a z-score 0077 % between them (=(mean1-mean2)/(sqrt(sd1^2+sd2^2))). 0078 %Add up all the z-scores (their absolute value) and have this function return 0079 % that value. 0080 % 0081 %Intuitively what we're doing here is comparing the two sets based on 0082 % how different the mdv's appear. 0083 %We're going to see if different glucose mixtures result in different values. 0084 %I'm going to rewrite part of slvrXXfast so it doesn't return every metabolite 0085 % but only those we can actually measure, but for now just make it generic. 0086 0087 function mdv = myidv2mdv (idv,xltmdv) 0088 0089 % generate the mdv 0090 len = length(idv); 0091 %disp(sprintf('idv is %d long',len)); 0092 mdv = zeros(1,xltmdv(len)+1); 0093 %disp(sprintf('mdv is %d long',length(mdv))); 0094 for i = 1:len 0095 idx = xltmdv(i) + 1; 0096 %disp(sprintf('idx is %d, currently on %d',idx,i)); 0097 mdv(idx) = mdv(idx) + idv(i); 0098 end 0099 0100 return 0101 end