Home > RAVEN > findGeneDeletions.m

findGeneDeletions

PURPOSE ^

findGeneDeletions

SYNOPSIS ^

function [genes fluxes originalGenes details]=findGeneDeletions(model,testType,analysisType,refModel,oeFactor)

DESCRIPTION ^

 findGeneDeletions
   Deletes genes, optimizes the model, and keeps track of the resulting
   fluxes. This is used for identifying gene deletion targets.
   
   model           a model structure
   testType        single/double gene deletions/over expressions. Over
                   expression only available if using MOMA
                   'sgd'   single gene deletion
                   'dgd'   double gene deletion
                   'sgo'   singel gene over expression
                   'dgo'   double gene over expression
   analysisType    determines whether to use FBA ('fba') or MOMA ('moma') 
                   in the optimization
   refModel        MOMA works by fitting the flux distributions of two
                   models to be as similar as possible. The most common
                   application is where you have a reference model where
                   some of the fluxes are constrained from experimental
                   data. This model is required when using MOMA
   oeFactor        a factor by which the fluxes should be increased if a
                   gene is overexpressed (opt, default 10)

   genes           a matrix with the genes that were deleted in each
                   optimization (the gene indexes in originalGenes). Each
                   row corresponds to a column in fluxes
   fluxes          a matrix with the resulting fluxes. Double deletions
                   that result in an unsolvable problem have all zero
                   flux. Single deletions that result in an unsolvable
                   problem are indicated in details instead
   originalGenes   simply the genes in the input model. Included for
                   simple presentation of the output
   details         not all genes will be deleted in all analyses. It is
                   for example not necessary to delete genes for dead end
                   reactions. This is a vector with details about 
                   each gene in originalGenes and why or why not it was 
                   deleted
                   1: Was deleted
                   2: Proved lethal in SGD
                   3: Only involved in reactions with too many iso-enzymes
                   4: Involved in dead-end reaction

   NOTE: This function disregards complexes. Any one gene can encode a
         reaction even if parts of the complex is deleted.

   Usage: [genes fluxes]=findGeneDeletions(model,testType,analysisType,...
           refModel,oeFactor)

   Rasmus Agren, 2010-12-16

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [genes fluxes originalGenes details]=findGeneDeletions(model,testType,analysisType,refModel,oeFactor)
0002 % findGeneDeletions
0003 %   Deletes genes, optimizes the model, and keeps track of the resulting
0004 %   fluxes. This is used for identifying gene deletion targets.
0005 %
0006 %   model           a model structure
0007 %   testType        single/double gene deletions/over expressions. Over
0008 %                   expression only available if using MOMA
0009 %                   'sgd'   single gene deletion
0010 %                   'dgd'   double gene deletion
0011 %                   'sgo'   singel gene over expression
0012 %                   'dgo'   double gene over expression
0013 %   analysisType    determines whether to use FBA ('fba') or MOMA ('moma')
0014 %                   in the optimization
0015 %   refModel        MOMA works by fitting the flux distributions of two
0016 %                   models to be as similar as possible. The most common
0017 %                   application is where you have a reference model where
0018 %                   some of the fluxes are constrained from experimental
0019 %                   data. This model is required when using MOMA
0020 %   oeFactor        a factor by which the fluxes should be increased if a
0021 %                   gene is overexpressed (opt, default 10)
0022 %
0023 %   genes           a matrix with the genes that were deleted in each
0024 %                   optimization (the gene indexes in originalGenes). Each
0025 %                   row corresponds to a column in fluxes
0026 %   fluxes          a matrix with the resulting fluxes. Double deletions
0027 %                   that result in an unsolvable problem have all zero
0028 %                   flux. Single deletions that result in an unsolvable
0029 %                   problem are indicated in details instead
0030 %   originalGenes   simply the genes in the input model. Included for
0031 %                   simple presentation of the output
0032 %   details         not all genes will be deleted in all analyses. It is
0033 %                   for example not necessary to delete genes for dead end
0034 %                   reactions. This is a vector with details about
0035 %                   each gene in originalGenes and why or why not it was
0036 %                   deleted
0037 %                   1: Was deleted
0038 %                   2: Proved lethal in SGD
0039 %                   3: Only involved in reactions with too many iso-enzymes
0040 %                   4: Involved in dead-end reaction
0041 %
0042 %   NOTE: This function disregards complexes. Any one gene can encode a
0043 %         reaction even if parts of the complex is deleted.
0044 %
0045 %   Usage: [genes fluxes]=findGeneDeletions(model,testType,analysisType,...
0046 %           refModel,oeFactor)
0047 %
0048 %   Rasmus Agren, 2010-12-16
0049 %
0050 originalModel=model;
0051 
0052 if nargin<5
0053     oeFactor=10;
0054 end
0055 
0056 %Check that the test type is correct
0057 if ~strcmpi(testType,'sgd') && ~strcmpi(testType,'dgd') && ~strcmpi(testType,'sgo') && ~strcmpi(testType,'dgo')
0058    throw(MException('','Incorrect test type')); 
0059 end
0060 
0061 %Check that the analysis type is correct
0062 if ~strcmpi(analysisType,'fba') && ~strcmpi(analysisType,'moma')
0063    throw(MException('','Incorrect analysis type')); 
0064 end
0065 
0066 if (strcmpi(testType,'sgo') || strcmpi(testType,'dgo')) && strcmpi(analysisType,'fba')
0067     throw(MException('','Over expression is only available when using MOMA')); 
0068 end
0069 
0070 if strcmpi(analysisType,'moma') && nargin<4
0071     throw(MException('','A reference model must be supplied when using MOMA')); 
0072 end
0073 
0074 originalGenes=model.genes;
0075 details=zeros(numel(model.genes),1);
0076 
0077 %First simplify the model to reduce the size
0078 model=simplifyModel(model,true,false,true,true);
0079 model=removeRxns(model,{},true,true); %Removes unused genes
0080 details(~ismember(originalGenes,model.genes))=4;
0081 
0082 [crap geneMapping]=ismember(model.genes,originalGenes);
0083 
0084 %Get the genes that should be deleted. Not all genes are deleted in all
0085 %optimizations. For example, if all reactions involving a gene have
0086 %iso-enzymes a single deletion of that gene will never have an effect.
0087 %FBA and SGD:   All genes that encode reactions on their own
0088 %FBA and DGD:   All genes that have at most one other gene encoding for at
0089 %               least one of its reactions and doesn't participate
0090 %               exclusively in reactions where a single deletion results in
0091 %               an unsolvable problem
0092 %MOMA and SGD:  All genes that encode reactions on their own
0093 %MOMA and DGD:  All combinations of genes that have at most one other gene
0094 %               encoding for at east one of its reactions and doesn't
0095 %               participate exclusively in reactions where a single
0096 %               deletion results in an unsolvable problem
0097 %MOMA and SGO:  All genes
0098 %MOMA and DGO:  All combinations of genes
0099 genesForSGD=[];
0100 if strcmpi(testType,'sgd') || strcmpi(testType,'dgd')
0101     if strcmpi(testType,'sgd')
0102         I=sum(model.rxnGeneMat,2)>1; %Reactions with iso-enzymes
0103     else
0104         I=sum(model.rxnGeneMat,2)>2; %Reactions with more than one iso-enzyme
0105     end
0106     
0107     %Find genes involved only in these reactions
0108     [crap inIsoRxns]=find(model.rxnGeneMat(I,:));
0109     [crap inNonIsoRxns]=find(model.rxnGeneMat(~I,:));
0110     genesForSGD=unique(inNonIsoRxns);
0111     details(geneMapping(setdiff(inIsoRxns,inNonIsoRxns)))=3;
0112 end
0113 
0114 %Get the genes that should be deleted in a SGD
0115 if strcmpi(testType,'sgd') || strcmpi(testType,'dgd')
0116    genesToModify=genesForSGD; 
0117 end
0118 
0119 %Get the genes that should be deleted in SGO
0120 if strcmpi(testType,'sgo')
0121    genesToModify=1:numel(model.genes);
0122    genesToModify=genesToModify(:);
0123 end
0124 
0125 %Do single deletion/over expression. This is done here since the double
0126 %deletion depends on which single deletions prove lethal
0127 %(to reduce the size of the system)
0128 if strcmpi(testType,'sgd') || strcmpi(testType,'sgo') || strcmpi(testType,'dgd')
0129     fluxes=zeros(numel(model.rxns),numel(genesToModify));
0130     solvable=true(numel(genesToModify),1);
0131     for i=1:numel(genesToModify)
0132        I=find(model.rxnGeneMat(:,genesToModify(i)));
0133        if strcmpi(testType,'sgd') || strcmpi(testType,'dgd')
0134            %Constrain all reactions involving the gene to 0
0135            tempModel=setParam(model,'eq',model.rxns(I),0);
0136        else
0137            %To over express a gene, the stoichiometry of the corresponding
0138            %reactions are changed so that the same flux leads to a higher
0139            %production
0140            tempModel=model;
0141            tempModel.S(:,I)=tempModel.S(:,I).*oeFactor;
0142        end
0143        if strcmpi(analysisType,'fba') || strcmpi(testType,'dgd')
0144             sol=solveLP(tempModel);
0145        else
0146             [fluxA crap flag]=qMOMA(tempModel,refModel);
0147             sol.x=fluxA;
0148             sol.stat=flag;
0149        end
0150        
0151        %If the optimization terminated successfully
0152        if sol.stat==1
0153            fluxes(:,i)=sol.x;
0154            details(geneMapping(genesToModify(i)))=1;
0155        else
0156            solvable(i)=false;
0157            details(geneMapping(genesToModify(i)))=2;
0158        end
0159     end
0160     
0161     fluxes=fluxes(:,solvable);
0162     genes=geneMapping(genesToModify(solvable));
0163 end
0164 
0165 %Now do for DGO. This is rather straight forward since it is always
0166 %solvable and it doesn't matter if there are iso-enzymes
0167 if strcmpi(testType,'dgo')
0168     genesToModify=nchoosek(1:numel(model.genes),2);
0169     genes=geneMapping(genesToModify);
0170     %Since I assume that this is never lethal I set the details already
0171     details(geneMapping)=1;
0172     
0173     fluxes=sparse(numel(model.rxns),size(genesToModify,1));
0174     for i=1:size(genesToModify,1)
0175        [I crap]=find(model.rxnGeneMat(:,genesToModify(i,:)));
0176        %To over express a gene, the stoichiometry of the corresponding
0177        %reactions are changed so that the same flux leads to a higher
0178        %production
0179        tempModel=model;
0180        tempModel.S(:,I)=tempModel.S(:,I).*oeFactor;
0181        [fluxA crap crap]=qMOMA(tempModel,refModel);
0182        fluxes(:,i)=fluxA;
0183     end
0184 end
0185 
0186 %For double deletions FBA or MOMA
0187 if strcmpi(testType,'dgd')
0188     %This is a little lazy but it's fine. Check which genes that have
0189     %already been labled as either ony with too many iso-enzymes or
0190     %non-solveable as a single deletion.
0191     [crap I]=ismember(originalGenes(details==1),model.genes);
0192     genesToModify=nchoosek(I,2);
0193     genes=geneMapping(genesToModify);
0194     
0195     fluxes=sparse(numel(model.rxns),size(genesToModify,1));
0196     for i=1:size(genesToModify,1)
0197        [I crap]=find(model.rxnGeneMat(:,genesToModify(i,:)));
0198 
0199        %Constrain all reactions involving the gene to 0
0200        tempModel=setParam(model,'eq',model.rxns(I),0);
0201        
0202        if strcmpi(analysisType,'fba')
0203             sol=solveLP(tempModel);
0204        else
0205             [fluxA crap flag]=qMOMA(tempModel,refModel);
0206             sol.x=fluxA;
0207             sol.stat=flag;
0208        end
0209        
0210        if sol.stat==1
0211            fluxes(:,i)=sol.x;
0212        end
0213     end
0214 end
0215 
0216 %Map back to the old model
0217 [crap I]=ismember(model.rxns,originalModel.rxns);
0218 temp=fluxes;
0219 fluxes=sparse(numel(originalModel.rxns),size(temp,2));
0220 fluxes(I,:)=temp;
0221 end

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005