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
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