Team:Alberta/Project/Modeling/FindMinimalGenome
From 2009.igem.org
(Difference between revisions)
(New page: function [MinimalGenomes, UnEssentialGeneLists]=FindMinimalGenomes(TheModel, NumberOfMinGenomes, FinalGrowthRate, InitialFilterK, NumberOfSteps) %NumberOfMinGenomes: Number of minimal ge...) |
|||
Line 1: | Line 1: | ||
- | function [MinimalGenomes, UnEssentialGeneLists]=FindMinimalGenomes(TheModel, NumberOfMinGenomes, FinalGrowthRate | + | I will likely replace this with a slightly updated version soon. |
+ | |||
+ | 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 |
- | + | %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)') | disp('May take 5 mins per NumberOfMinGenomes (default is 1 list)') | ||
solutionI=optimizeCbModel(TheModel); | 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') %(unless you expect some gene deletions to cause better growth) | ||
+ | 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 | 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 "steps (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 | |
- | + | fStep = (f1-(f2)-0.0001)/n | |
- | + | 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 | end | ||
+ | fStep = (f1-(f2)-0.0001)/n | ||
MinimalGenomes=cell(1,NumberOfMinGenomes); | MinimalGenomes=cell(1,NumberOfMinGenomes); | ||
UnEssentialGeneLists=cell(1,NumberOfMinGenomes); | UnEssentialGeneLists=cell(1,NumberOfMinGenomes); | ||
TheModelOriginal = TheModel; | TheModelOriginal = TheModel; | ||
+ | Models=cell(1); | ||
+ | |||
w=1; | w=1; | ||
+ | %Algorithm to find minimal genomes. | ||
+ | |||
while w<=NumberOfMinGenomes | 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 | if solution.f>f %"if the deletion causes satisfactory growth, commit gene to geneList | ||
geneList(d)=geneListTEST(d); | geneList(d)=geneListTEST(d); | ||
Line 85: | Line 129: | ||
d=d+1; | d=d+1; | ||
end | 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 | end | ||
- | %University of Alberta 2009 iGEM | + | %University of Alberta 2009 iGEM Team |
- | % | + | %Project BioBytes |
+ | %EricB |
Revision as of 06:01, 8 October 2009
I will likely replace this with a slightly updated version soon.
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') %(unless you expect some gene deletions to cause better growth) 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 "steps (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 fStep = (f1-(f2)-0.0001)/n 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