|
|
(6 intermediate revisions not shown) |
Line 1: |
Line 1: |
- | 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
| + | |
- | 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 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
| + | |