Team:Alberta/Project/Modeling/FindMinimalGenomes

From 2009.igem.org

Latest revision as of 00:41, 22 October 2009

function [MinimalGenomes, UnEssentialGeneLists, Models] = FindMinimalGenomes(TheModel, NumberOfMinGenomes, FinalGrowthRate, 
IterationsPerList, 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 IterationsPerList < 1
       error('IterationsPerList 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(IterationsPerList,1) ~= 0
            error('IterationsPerList must be an integer')
        end
        FinalGrowthRate;
        FirstThreshold = ((InitialGrowthRate-FinalGrowthRate)*0.87)+FinalGrowthRate;
        if IterationsPerList == 1
            n = 0.000000001; %Causes only 1 iteration & prevents division by 0 (which only Chuck Norris can do).
        else
        n = IterationsPerList-1;
        end
    case 5
        if rem(IterationsPerList,1) ~= 0
            error('IterationsPerList must be an integer')
        end
        FinalGrowthRate;
        FirstThreshold;
        if IterationsPerList == 1
            n = 0.000000001; %to cause only 1 iteration & to prevent division by 0
        else
        n = IterationsPerList-1;
        if FirstThreshold < FinalGrowthRate
                error('FirstThreshold must be greater than FinalGrowthRate')
        end
        end
    otherwise
        error('Too many input arguments! Sorry, try again!')
end

ThreshDrop = (FirstThreshold-(FinalGrowthRate)-0.0001)/n;  
MinimalGenomes = cell(1,NumberOfMinGenomes);
UnEssentialGeneLists = cell(1,NumberOfMinGenomes);
TheModelOriginal = TheModel;
Models = cell(1);

disp(['Initial Growth Rate is ' num2str(InitialGrowthRate)])
disp(['Final Growth Rate will be ~ ' num2str(FinalGrowthRate)])
disp(['First of ' num2str(IterationsPerList) ' growth rate thresholds is ' num2str(FirstThreshold)])
disp(['The growth rate threshold will drop by this amount for every iteration: ' num2str(ThreshDrop)])
disp(['Number of iterations for every list is ' num2str(IterationsPerList)])
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-ThreshDrop;
        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'
    
    if nargout > 2
        Models{w} = deleteModelGenes(TheModelOriginal, geneList);
    end
    w = w+1;
end

if nargout > 2
    disp('Reassign a model of interest like this: modelMinimal = Models{3}, except replace "Models" with your 3rd output 
 argument.')
end

%University of Alberta 2009 iGEM Team
%Project BioBytes
%EricB

• Recent changes
• What links here
• Related changes
• Special pages
• My preferences

• Printable version
• Permanent link
• Privacy policy
• Disclaimers