From 2009.igem.org
function [MinimalGenomes, UnEssentialGeneLists]=FindMinimalGenomes(TheModel, NumberOfMinGenomes, FinalGrowthRate, InitialFilterK, NumberOfSteps)
%NumberOfMinGenomes: Number of minimal genomes you want the code to find (since there are many possible minimal genomes)
%InitialFilterK: This constant (<1) times the initial growth rate gives you the minimum survival rate acceptible for the first gene deletion iteration.
%FinalGrowthRate is what the final growth rate will be. Lower values result in a worse final growth rate, but more genes can be deleted
%(resulting in a smaller allowable minimal genome).
disp('May take 5 mins per NumberOfMinGenomes (default is 1 list)')
solutionI=optimizeCbModel(TheModel);
fValueInitial = solutionI.f
switch nargin
case 1
f1 = fValueInitial*0.87
f2 = 0.3
n = 2 % n=2 results in 3 "steps (3 iterations of gene filtering)"
f3 = (f1-(f2)-0.001)/n
NumberOfMinGenomes=1;
case 2
f1 = fValueInitial*0.87
f2 = 0.3
n = 2
f3 = (f1-(f2)-0.001)/n
case 3
f1 = fValueInitial*0.87
f2 = FinalGrowthRate
n = 2
f3 = (f1-(f2)-0.001)/n
case 4
f1 = fValueInitial*InitialFilterK
f2 = FinalGrowthRate
n = 2
f3 = (f1-(f2)-0.001)/n
case 5
f1 = fValueInitial*InitialFilterK
f2 = FinalGrowthRate
n = NumberOfSteps-1
f3 = (f1-(f2)-0.001)/n
otherwise
error('Too many input arguments!')
end
MinimalGenomes=cell(1,NumberOfMinGenomes);
UnEssentialGeneLists=cell(1,NumberOfMinGenomes);
TheModelOriginal = TheModel;
w=1;
while w<=NumberOfMinGenomes
w
f=f1;
geneList=cell(1);
geneListTEST=cell(1);
d=1;
RandomGeneList=cell(1,length(TheModelOriginal.genes));
%this part shuffles the TheModel list for randomization
random=randperm(length(TheModel.genes))';
for x=(1:length(TheModel.genes))
RandomGeneList(x)=TheModelOriginal.genes(random(x));
end
while f>=f2
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
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:'), disp(length(geneList))
f=f-f3;
end
UnEssentialGeneLists{w} = geneList';
%"inverts" the unessential list to obtain the essential list
InvertThis = UnEssentialGeneLists{w};
List2=TheModel.genes;
InvertedList=cell(1);
Match=zeros(length(List2));
for a=(1:length(InvertThis))
for b=(1:length(List2))
if strcmp(InvertThis(a), List2(b)) == 1
Match(b)=1;
end
end
end
d=1;
for c=(1:length(List2))
if Match(c)==0
InvertedList(d)=List2(c);
d=d+1;
end
end
MinimalGenomes{w}=InvertedList'
w=w+1;
end
%University of Alberta 2009 iGEM team
%Eric B
"
