Team:Alberta/Project/Modeling/FindMinimalGenome

From 2009.igem.org

(Difference between revisions)
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