% This script computes the maximal sample size we can reconstruct for a given number of lanes,
% or the minimal number of lanes sufficient to reconstruct a certain number of individuals.
% The struct auxStruct contains all relevant simulation parameters and then runs the program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
auxStruct = struct;
auxStruct.simulate_direction = 'lanes_to_people'; % 'lanes_to_people' finds the maximal number of
% individuals that can be reconstructed with zero
% errors for a given number of pools
% 'people_to_lane' find the minimal number of lanes
% needed to reconstruct a certain number of individuals
auxStruct.maxNumTrials = 500; % maximal number of simulations for each parameter setting
auxStruct.sigma_pooling = 0.05; % sigma of DNA preparation errors.
auxStruct.read_error_prob = 0.01; % prob. base read is different from true base
auxStruct.minFractionOfFailure = 0.05; % maximal error level allowed.
%0.05 means that in at least 95% of the simulations one achieved
%zero error reconstruction
auxStruct.forceBBflag = 0; %if equals 'replace half by 2', then half of the
%non-zero entries in the genotype vector are replaces
%by the value 2, to simulate the "BB" case, .i.e homozygous rare-allele.
%Any other value for forceBBflag simulates the regular "AB" case.
%The default value of forceBBflag corresponds to the "AB" case.
auxStruct.single_region_length = 300; % average length of a region
auxStruct.read_length = 100; % length of each read
auxStruct.carrier_freq_vec = 1 ./ 100; % % try different carrier allele frequencies
auxStruct.num_people_vec = [500 1000 2000 4000]; % try different numbers of individuals
auxStruct.num_pools_vec = 10; % number of available pools/lanes
auxStruct.num_regions_vec = 1; % [1 10 100]; % determine regions size (and coverage)
auxStruct.num_people_in_pool_vec = 25; % [25 100 0.5]; % maximal number of people in one pool (fraction means frequency of people in pool)
auxStruct.num_barcodes = 1; % "1" means no barcodes
auxStruct.reads = 4*10^6 * auxStruct.read_length; % number used in the paper (here reads are measured in base-pairs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform the simulation
optimize_experimental_setting(auxStruct)
% the output appears in a file - see directory.