From 2009.igem.org
function [MinimalGenomes, UnEssentialGeneLists, Models] = FindMinimalGenomes(TheModel, NumberOfMinGenomes, FinalGrowthRate, NumberOfSteps, FirstThreshold)
%NumberOfMinGenomes: Number of minimal genomes you want the code to find (since there are many possible minimal genomes)
%FinalGrowthRate: is what the final growth rate will apporximately be. Lower values result in a lower final
%growth rate, but this allows for more genes to be deleted (resulting in a smaller allowable minimal genome).
%FirstThreshold: This is the threshold growth rate for the first iteration of gene deletions.
disp('May take 5 mins per NumberOfMinGenomes (default is 1 list)')
solutionI = optimizeCbModel(TheModel);
InitialGrowthRate = solutionI.f
%Input/Output Error Checking
if nargout == 0
error('You must have at least 1 output argument! For example: MinGenomes = FindMinimalGenomes(model)')
end
if FinalGrowthRate > InitialGrowthRate
disp(['FinalGrowthRate (set by you): ' num2str(FinalGrowthRate)]);
error('The FinalGrowthRate must be less than the InitialGrowthRate')
end
if NumberOfSteps < 1
error('NumberOfSteps must be greater than or equal to 1!')
end
if FirstThreshold >= InitialGrowthRate
error('FirstThreshold must be less than the InitialGrowthRate')
end
%Initialization of variables
switch nargin
case 1
FinalGrowthRate = 0.2
f1 = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate
f2 = FinalGrowthRate;
n = 2; iterationsPerList = n+1 % n = 2 results in 3 iterations of gene filtering.
NumberOfMinGenomes = 1;
case 2
FinalGrowthRate = 0.2
f1 = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate
f2 = 0.3
n = 2; iterationsPerList = n+1
case 3
FinalGrowthRate
f1 = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate
f2 = FinalGrowthRate
n = 2; iterationsPerList = n+1
case 4
if rem(NumberOfSteps,1) ~= 0
error('NumberOfSteps must be an integer')
end
f1 = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate
f2 = FinalGrowthRate
if NumberOfSteps == 1
n = 0.000000001; %to cause only 1 iteration & to prevent division by 0 (which only Chuck Norris can do).
iterationsPerList = 1
else
n = NumberOfSteps-1;
iterationsPerList = NumberOfSteps
end
case 5
if rem(NumberOfSteps,1) ~= 0
error('NumberOfSteps must be an integer')
end
f1 = FirstThreshold
f2 = FinalGrowthRate
if NumberOfSteps == 1
n = 0.000000001; %to cause only 1 iteration & to prevent division by 0
iterationsPerList = 1
else
n = NumberOfSteps-1;
iterationsPerList = NumberOfSteps
end
otherwise
error('Too many input arguments! Sorry, try again!')
end
fStep = (f1-(f2)-0.0001)/n
MinimalGenomes = cell(1,NumberOfMinGenomes);
UnEssentialGeneLists = cell(1,NumberOfMinGenomes);
TheModelOriginal = TheModel;
Models = cell(1);
w = 1;
%Algorithm to find minimal genomes.
while w<=NumberOfMinGenomes
disp(['Finding Minimal Genome List: ' num2str(w)]);
geneList = cell(1);
geneListTEST = cell(1);
d = 1;
f = f1;
%this part shuffles the TheModel list for randomization
RandomGeneList = cell(1,length(TheModelOriginal.genes));
random = randperm(length(TheModel.genes))';
for x = (1:length(TheModel.genes))
RandomGeneList(x) = TheModelOriginal.genes(random(x));
end
while f >= f2
disp(['Performing deletions that result in a growth rate of at least ' num2str(f) ' ...']);
for c = (1:length(RandomGeneList))
match = zeros(length(RandomGeneList));
for a = (1:length(geneList)) %check if gene already on geneList. if so, skip analysis
if strcmp(geneList(a),RandomGeneList(c)) == 1
match(a) = 1;
else
match(a) = 0;
end
end
%This part decides whether or not to permanently knockout a gene
if sum(match) == 0 %if gene does not exist in geneList yet, test its deletion
geneListTEST(d) = RandomGeneList(c);
deltamodel = deleteModelGenes(TheModelOriginal,geneListTEST);
solution = optimizeCbModel(deltamodel);
if solution.f > f %"if the deletion causes satisfactory growth, commit gene to geneList
geneList(d) = geneListTEST(d);
%F(d) = solution.f;
d = d+1;
end
end
end
disp(['Genes Deleted So Far for this List is: ' num2str(length(geneList))])
f = f-fStep;
end
UnEssentialGeneLists{w} = geneList';
%This part "inverts" the unessential list to obtain the essential list
InvertThis = UnEssentialGeneLists{w};
WithThis = TheModel.genes;
InvertedList = cell(1);
Match=zeros(length(WithThis));
for a = (1:length(InvertThis))
for b = (1:length(WithThis))
if strcmp(InvertThis(a), WithThis(b)) == 1
Match(b) = 1;
end
end
end
d = 1;
for c=(1:length(WithThis))
if Match(c) == 0
InvertedList(d) = WithThis(c);
d = d+1;
end
end
MinimalGenomes{w} = InvertedList'
Models{w} = deleteModelGenes(TheModelOriginal, geneList);
w = w+1;
end
if nargout == 3
disp('Access a model of interest like this: model3 = Models{3}, except replace "Models" with your 3rd output argument.')
end
%University of Alberta 2009 iGEM Team
%Project BioBytes
%EricB
"
