Team:Alberta/Project/Modeling/FindMinimalGenome
From 2009.igem.org
(Difference between revisions)
Line 1: | Line 1: | ||
- | function [MinimalGenomes, UnEssentialGeneLists, Models] = FindMinimalGenomes(TheModel, NumberOfMinGenomes, FinalGrowthRate, NumberOfSteps, FirstThreshold) | + | 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) | %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 | %FinalGrowthRate: is what the final growth rate will apporximately be. Lower values result in a lower final | ||
Line 36: | Line 36: | ||
case 1 | case 1 | ||
FinalGrowthRate = 0.2 | FinalGrowthRate = 0.2 | ||
- | + | FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate | |
- | + | ||
n = 2; iterationsPerList = n+1 % n = 2 results in 3 iterations of gene filtering. | n = 2; iterationsPerList = n+1 % n = 2 results in 3 iterations of gene filtering. | ||
NumberOfMinGenomes = 1; | NumberOfMinGenomes = 1; | ||
case 2 | case 2 | ||
FinalGrowthRate = 0.2 | FinalGrowthRate = 0.2 | ||
- | + | FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate | |
- | + | ||
n = 2; iterationsPerList = n+1 | n = 2; iterationsPerList = n+1 | ||
case 3 | case 3 | ||
FinalGrowthRate | FinalGrowthRate | ||
- | + | FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate | |
- | + | ||
n = 2; iterationsPerList = n+1 | n = 2; iterationsPerList = n+1 | ||
case 4 | case 4 | ||
Line 54: | Line 51: | ||
error('NumberOfSteps must be an integer') | error('NumberOfSteps must be an integer') | ||
end | end | ||
- | + | FinalGrowthRate | |
- | + | FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate | |
if NumberOfSteps == 1 | if NumberOfSteps == 1 | ||
n = 0.000000001; %to cause only 1 iteration & to prevent division by 0 (which only Chuck Norris can do). | n = 0.000000001; %to cause only 1 iteration & to prevent division by 0 (which only Chuck Norris can do). | ||
Line 67: | Line 64: | ||
error('NumberOfSteps must be an integer') | error('NumberOfSteps must be an integer') | ||
end | end | ||
- | + | FinalGrowthRate | |
- | + | FirstThreshold | |
if NumberOfSteps == 1 | if NumberOfSteps == 1 | ||
n = 0.000000001; %to cause only 1 iteration & to prevent division by 0 | n = 0.000000001; %to cause only 1 iteration & to prevent division by 0 | ||
Line 75: | Line 72: | ||
n = NumberOfSteps-1; | n = NumberOfSteps-1; | ||
iterationsPerList = NumberOfSteps | iterationsPerList = NumberOfSteps | ||
+ | if FirstThreshold < FinalGrowthRate | ||
+ | error('FirstThreshold must be greater than FinalGrowthRate') | ||
+ | end | ||
end | end | ||
otherwise | otherwise | ||
Line 80: | Line 80: | ||
end | end | ||
- | fStep = ( | + | fStep = (FirstThreshold-(FinalGrowthRate)-0.0001)/n |
MinimalGenomes = cell(1,NumberOfMinGenomes); | MinimalGenomes = cell(1,NumberOfMinGenomes); | ||
UnEssentialGeneLists = cell(1,NumberOfMinGenomes); | UnEssentialGeneLists = cell(1,NumberOfMinGenomes); | ||
Line 91: | Line 91: | ||
while w<=NumberOfMinGenomes | while w<=NumberOfMinGenomes | ||
- | disp(['Finding Minimal Genome List: ' num2str(w)]); | + | disp(['Finding Minimal Genome List: ' num2str(w) ' of ' num2str(NumberOfMinGenomes)]); |
geneList = cell(1); | geneList = cell(1); | ||
geneListTEST = cell(1); | geneListTEST = cell(1); | ||
d = 1; | d = 1; | ||
- | f = | + | k = 1; |
+ | f = FirstThreshold; | ||
%this part shuffles the TheModel list for randomization | %this part shuffles the TheModel list for randomization | ||
Line 107: | Line 108: | ||
end | end | ||
- | while f >= | + | while f >= FinalGrowthRate |
- | disp(['Performing | + | disp(['Performing iteration #' num2str(k) '. Threshold growth rate for this iteration of deletions is ' num2str(f) ' ...']); |
- | + | ||
for c = (1:length(RandomGeneList)) | for c = (1:length(RandomGeneList)) | ||
match = zeros(length(RandomGeneList)); | match = zeros(length(RandomGeneList)); | ||
Line 137: | Line 137: | ||
disp(['Genes Deleted So Far for this List is: ' num2str(length(geneList))]) | disp(['Genes Deleted So Far for this List is: ' num2str(length(geneList))]) | ||
f = f-fStep; | f = f-fStep; | ||
+ | k = k+1; | ||
end | end | ||
Line 155: | Line 156: | ||
end | end | ||
end | end | ||
- | d = 1; | + | |
- | + | d = 1; | |
for c=(1:length(WithThis)) | for c=(1:length(WithThis)) | ||
if Match(c) == 0 | if Match(c) == 0 |
Revision as of 02:42, 9 October 2009
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 nargin >=3 if FinalGrowthRate > InitialGrowthRate disp(['FinalGrowthRate (set by you): ' num2str(FinalGrowthRate)]); error('The FinalGrowthRate must be less than the InitialGrowthRate') end end if nargin >=4 if NumberOfSteps < 1 error('NumberOfSteps must be greater than or equal to 1!') end end if nargin == 5 if FirstThreshold >= InitialGrowthRate error('FirstThreshold must be less than the InitialGrowthRate') end end %Initialization of variables switch nargin case 1 FinalGrowthRate = 0.2 FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate n = 2; iterationsPerList = n+1 % n = 2 results in 3 iterations of gene filtering. NumberOfMinGenomes = 1; case 2 FinalGrowthRate = 0.2 FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate n = 2; iterationsPerList = n+1 case 3 FinalGrowthRate FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate n = 2; iterationsPerList = n+1 case 4 if rem(NumberOfSteps,1) ~= 0 error('NumberOfSteps must be an integer') end FinalGrowthRate FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+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 FinalGrowthRate FirstThreshold if NumberOfSteps == 1 n = 0.000000001; %to cause only 1 iteration & to prevent division by 0 iterationsPerList = 1 else n = NumberOfSteps-1; iterationsPerList = NumberOfSteps if FirstThreshold < FinalGrowthRate error('FirstThreshold must be greater than FinalGrowthRate') end end otherwise error('Too many input arguments! Sorry, try again!') end fStep = (FirstThreshold-(FinalGrowthRate)-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) ' of ' num2str(NumberOfMinGenomes)]); geneList = cell(1); geneListTEST = cell(1); d = 1; k = 1; f = FirstThreshold; %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 >= FinalGrowthRate disp(['Performing iteration #' num2str(k) '. Threshold growth rate for this iteration of deletions is ' 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; k = k+1; end UnEssentialGeneLists{w} = geneList'; %This part "inverts" the unessential list of genes to obtain the essential list (the minimal genome). 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