sampleCbModel Sample the solution-space of a constraint-based model [modelSampling,samples] = sampleCbModel(model,sampleFile,samplerName,options) INPUTS model COBRA model structure sampleFile File names for sampling output files OPTIONAL INPUTS samplerName Name of the sampler to be used to sample the solution (default 'ACHR') options Options for sampling and pre/postprocessing (default values in parenthesis) nWarmupPoints Number of warmup points (5000) nFiles Number of output files (10) nPointsPerFile Number of points per file (1000) nStepsPerPoint Number of sampler steps per point saved (200) nPointsReturned Number of points loaded for analysis (2000) nFilesSkipped Number of output files skipped when loading points to avoid potentially biased initial samples (2) removeLoopsFlag Attempt to remove loops from the model (false) removeLoopSamplesFlag Remove sampling data for reactions involved in loops (true) OUTPUTS modelSampling Cleaned up model used in sampling samples Uniform random samples of the solution space Examples of use: 1) Sample a model called 'superModel' using default settings and save the results in files with the common beginning 'superModelSamples' [modelSampling,samples] = sampleCbModel(superModel,'superModelSamples'); 2) Sample a model called 'hyperModel' using default settings except with a total of 50 sample files saved and with 5000 sample points returned. Save the results in files with the common beginning 'hyperModelSamples' options.nFiles = 50; options.nPointsReturned = 5000; [modelSampling,samples] = sampleCbModel(hyperModel,'hyperModelSamples'); Markus Herrgard 8/14/06
0001 function [modelSampling,samples] = sampleCbModel(model,sampleFile,samplerName,options) 0002 %sampleCbModel Sample the solution-space of a constraint-based model 0003 % 0004 % [modelSampling,samples] = sampleCbModel(model,sampleFile,samplerName,options) 0005 % 0006 %INPUTS 0007 % model COBRA model structure 0008 % sampleFile File names for sampling output files 0009 % 0010 %OPTIONAL INPUTS 0011 % samplerName Name of the sampler to be used to sample the solution 0012 % (default 'ACHR') 0013 % options Options for sampling and pre/postprocessing (default values 0014 % in parenthesis) 0015 % nWarmupPoints Number of warmup points (5000) 0016 % nFiles Number of output files (10) 0017 % nPointsPerFile Number of points per file (1000) 0018 % nStepsPerPoint Number of sampler steps per point saved (200) 0019 % nPointsReturned Number of points loaded for analysis (2000) 0020 % nFilesSkipped Number of output files skipped when loading 0021 % points to avoid potentially biased initial samples (2) 0022 % removeLoopsFlag Attempt to remove loops from the model (false) 0023 % removeLoopSamplesFlag Remove sampling data for reactions involved in 0024 % loops (true) 0025 % 0026 %OUTPUTS 0027 % modelSampling Cleaned up model used in sampling 0028 % samples Uniform random samples of the solution space 0029 % 0030 % Examples of use: 0031 % 0032 % 1) Sample a model called 'superModel' using default settings and save the 0033 % results in files with the common beginning 'superModelSamples' 0034 % 0035 % [modelSampling,samples] = sampleCbModel(superModel,'superModelSamples'); 0036 % 0037 % 2) Sample a model called 'hyperModel' using default settings except with a total of 50 sample files 0038 % saved and with 5000 sample points returned. 0039 % Save the results in files with the common beginning 'hyperModelSamples' 0040 % 0041 % options.nFiles = 50; 0042 % options.nPointsReturned = 5000; 0043 % [modelSampling,samples] = sampleCbModel(hyperModel,'hyperModelSamples'); 0044 % 0045 % Markus Herrgard 8/14/06 0046 0047 % Default options 0048 nWarmupPoints = 5000; 0049 nFiles = 10; 0050 nPointsPerFile = 1000; 0051 nStepsPerPoint = 200; 0052 nPointsReturned = 2000; 0053 nFilesSkipped = 2; 0054 removeLoopsFlag = false; 0055 removeLoopSamplesFlag = true; 0056 0057 if (nargin < 3) 0058 samplerName = 'ACHR'; 0059 end 0060 0061 % Handle options 0062 if (nargin > 3) 0063 if (isfield(options,'nWarmupPoints')) 0064 nWarmupPoints = options.nWarmupPoints; 0065 end 0066 if (isfield(options,'nFiles')) 0067 nFiles = options.nFiles; 0068 end 0069 if (isfield(options,'nPointsPerFile')) 0070 nPointsPerFile = options.nPointsPerFile; 0071 end 0072 if (isfield(options,'nStepsPerPoint')) 0073 nStepsPerPoint = options.nStepsPerPoint; 0074 end 0075 if (isfield(options,'nPointsReturned')) 0076 nPointsReturned = options.nPointsReturned; 0077 end 0078 if (isfield(options,'nFilesSkipped')) 0079 nFilesSkipped = options.nFilesSkipped; 0080 end 0081 if (isfield(options,'removeLoopsFlag')) 0082 removeLoopsFlag = options.removeLoopsFlag; 0083 end 0084 if (isfield(options,'removeLoopSamplesFlag')) 0085 removeLoopSamplesFlag = options.removeLoopSamplesFlag; 0086 end 0087 end 0088 0089 0090 fprintf('Prepare model for sampling\n'); 0091 % Prepare model for sampling by reducing bounds 0092 [nMet,nRxn] = size(model.S); 0093 fprintf('Original model: %d rxns %d metabolites\n',nRxn,nMet); 0094 0095 % Reduce model 0096 fprintf('Reduce model\n'); 0097 modelRed = reduceModel(model, 1e-6, false,false); 0098 [nMet,nRxn] = size(modelRed.S); 0099 fprintf('Reduced model: %d rxns %d metabolites\n',nRxn,nMet); 0100 save modelRedTmp modelRed; 0101 0102 modelSampling = modelRed; 0103 0104 switch samplerName 0105 case 'ACHR' 0106 % Use Artificial Centering Hit-and-run 0107 0108 fprintf('Create warmup points\n'); 0109 % Create warmup points for sampler 0110 warmupPts= createHRWarmup(modelSampling,nWarmupPoints); 0111 0112 save sampleCbModelTmp modelSampling warmupPts 0113 0114 fprintf('Run sampler for a total of %d steps\n',nFiles*nPointsPerFile*nStepsPerPoint); 0115 % Sample model 0116 ACHRSampler(modelSampling,warmupPts,sampleFile,nFiles,nPointsPerFile,nStepsPerPoint); 0117 0118 otherwise 0119 error(['Unknown sampler: ' samplerName]); 0120 end 0121 0122 fprintf('Load samples\n'); 0123 % Load samples 0124 nPointsPerFileLoaded = ceil(nPointsReturned/(nFiles-nFilesSkipped)); 0125 if (nPointsPerFileLoaded > nPointsPerFile) 0126 error('Attempted to return more points than were saved'); 0127 end 0128 samples = loadSamples(sampleFile,nFiles,nPointsPerFileLoaded,nFilesSkipped); 0129 samples = samples(:,round(linspace(1,size(samples,2),nPointsReturned))); 0130 % Fix reaction directions 0131 [modelSampling,samples] = convRevSamples(modelSampling,samples);