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' single gene over expression
'dgo' double gene over expression
(optional, default 'sgd')
analysisType determines whether to use FBA ('fba') or MOMA ('moma')
in the optimization. (optional, default 'fba')
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 (optional, 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/overexpressed
2: Proved lethal in sgd (single gene deletion)
3: - redundant, no longer used -
4: Involved in dead-end reaction
grRatioMuts growth rate ratio between mutated strain and wild type,
matches the originalGenes(genes) mutants. Note that
this does not directly map to model.genes, as is the case
for COBRA getEssentialGenes. However, this can be
obtained by afterwards running:
grRatio=zeros(1,numel(model.genes));
grRatio(genes)=grRatioMuts;
Usage: [genes, fluxes, originalGenes, details, grRatioMuts]=findGeneDeletions(model,testType,analysisType,...
refModel,oeFactor)0001 function [genes, fluxes, originalGenes, details, grRatioMuts]=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' single gene over expression 0012 % 'dgo' double gene over expression 0013 % (optional, default 'sgd') 0014 % analysisType determines whether to use FBA ('fba') or MOMA ('moma') 0015 % in the optimization. (optional, default 'fba') 0016 % refModel MOMA works by fitting the flux distributions of two 0017 % models to be as similar as possible. The most common 0018 % application is where you have a reference model where 0019 % some of the fluxes are constrained from experimental 0020 % data. This model is required when using MOMA 0021 % oeFactor a factor by which the fluxes should be increased if a 0022 % gene is overexpressed (optional, default 10) 0023 % 0024 % genes a matrix with the genes that were deleted in each 0025 % optimization (the gene indexes in originalGenes). Each 0026 % row corresponds to a column in fluxes 0027 % fluxes a matrix with the resulting fluxes. Double deletions 0028 % that result in an unsolvable problem have all zero 0029 % flux. Single deletions that result in an unsolvable 0030 % problem are indicated in details instead 0031 % originalGenes simply the genes in the input model. Included for 0032 % simple presentation of the output 0033 % details not all genes will be deleted in all analyses. It is 0034 % for example not necessary to delete genes for dead end 0035 % reactions. This is a vector with details about 0036 % each gene in originalGenes and why or why not it was 0037 % deleted 0038 % 1: Was deleted/overexpressed 0039 % 2: Proved lethal in sgd (single gene deletion) 0040 % 3: - redundant, no longer used - 0041 % 4: Involved in dead-end reaction 0042 % grRatioMuts growth rate ratio between mutated strain and wild type, 0043 % matches the originalGenes(genes) mutants. Note that 0044 % this does not directly map to model.genes, as is the case 0045 % for COBRA getEssentialGenes. However, this can be 0046 % obtained by afterwards running: 0047 % grRatio=zeros(1,numel(model.genes)); 0048 % grRatio(genes)=grRatioMuts; 0049 % 0050 % Usage: [genes, fluxes, originalGenes, details, grRatioMuts]=findGeneDeletions(model,testType,analysisType,... 0051 % refModel,oeFactor) 0052 0053 originalModel=model; 0054 if nargin<5 0055 oeFactor=10; 0056 end 0057 if nargin<2 0058 testType='sgd'; 0059 else 0060 testType=char(testType); 0061 end 0062 0063 %Check that the test type is correct 0064 if ~strcmpi(testType,'sgd') && ~strcmpi(testType,'dgd') && ~strcmpi(testType,'sgo') && ~strcmpi(testType,'dgo') 0065 EM='Incorrect test type'; 0066 dispEM(EM); 0067 end 0068 0069 %Check that the analysis type is correct 0070 if nargin<3 0071 analysisType = 'fba'; 0072 else 0073 analysisType=char(analysisType); 0074 if ~any(strcmpi(analysisType,{'fba','moma'})) 0075 EM='Incorrect analysis type'; 0076 dispEM(EM); 0077 end 0078 end 0079 0080 if (strcmpi(testType,'sgo') || strcmpi(testType,'dgo')) && strcmpi(analysisType,'fba') 0081 EM='Over expression is only available when using MOMA'; 0082 dispEM(EM); 0083 end 0084 0085 if strcmpi(analysisType,'moma') && nargin<4 0086 EM='A reference model must be supplied when using MOMA'; 0087 dispEM(EM); 0088 end 0089 0090 originalGenes=model.genes; 0091 details=zeros(numel(model.genes),1); 0092 0093 %First simplify the model to reduce the size 0094 model=simplifyModel(model,true,false,true,true); 0095 model=removeReactions(model,{},true,true); %Removes unused genes 0096 details(~ismember(originalGenes,model.genes))=4; 0097 0098 [~, geneMapping]=ismember(model.genes,originalGenes); 0099 growthWT=solveLP(model); 0100 growthWT=growthWT.f; 0101 0102 %Do single deletion/over expression. This is done here since the double 0103 %deletion depends on which single deletions prove lethal (to reduce the 0104 %size of the system) 0105 if strcmpi(testType,'sgd') || strcmpi(testType,'sgo') || strcmpi(testType,'dgd') 0106 fluxes=zeros(numel(model.rxns),numel(model.genes)); 0107 grRatioMuts=zeros(1,numel(model.genes)); 0108 for i=1:numel(model.genes) 0109 if strcmpi(testType,'sgd') || strcmpi(testType,'dgd') 0110 %Constrain all reactions involving the gene to 0 0111 tempModel=removeGenes(model,i,false,false,false); 0112 else 0113 %To over express a gene, the stoichiometry of the corresponding 0114 %reactions are changed so that the same flux leads to a higher 0115 %production 0116 tempModel=model; 0117 I=find(model.rxnGeneMat(:,i)); 0118 tempModel.S(:,I)=tempModel.S(:,I).*oeFactor; 0119 end 0120 if strcmpi(analysisType,'fba') || strcmpi(testType,'dgd') 0121 sol=solveLP(tempModel); 0122 else 0123 [fluxA, ~, flag]=qMOMA(tempModel,refModel); 0124 sol.x=fluxA; 0125 sol.stat=flag; 0126 end 0127 0128 %If the optimization terminated successfully 0129 if sol.stat==1 0130 fluxes(:,i)=sol.x; 0131 grRatioMuts(i)=sol.f/growthWT; 0132 details(geneMapping(i))=1; 0133 else 0134 fluxes(:,i)=0; 0135 grRatioMuts(i)=0; 0136 details(geneMapping(i))=2; 0137 end 0138 end 0139 genes=geneMapping; 0140 end 0141 0142 %Now do for DGO. This is rather straight forward since it is always 0143 %solvable and it doesn't matter if there are iso-enzymes 0144 if strcmpi(testType,'dgo') 0145 genesToModify=nchoosek(1:numel(model.genes),2); 0146 genes=geneMapping(genesToModify); 0147 %Since I assume that this is never lethal I set the details already 0148 details(geneMapping)=1; 0149 0150 fluxes=sparse(numel(model.rxns),size(genesToModify,1)); 0151 for i=1:size(genesToModify,1) 0152 I=find(model.rxnGeneMat(:,genesToModify(i,:))); 0153 %To over express a gene, the stoichiometry of the corresponding 0154 %reactions are changed so that the same flux leads to a higher 0155 %production 0156 tempModel=model; 0157 tempModel.S(:,I)=tempModel.S(:,I).*oeFactor; 0158 fluxA=qMOMA(tempModel,refModel); 0159 fluxes(:,i)=fluxA; 0160 grRatioMuts(i)=fluxA(logical(model.c))/growthWT; 0161 end 0162 end 0163 0164 %For double deletions FBA or MOMA 0165 if strcmpi(testType,'dgd') 0166 %This is a little lazy but it's fine. Check which genes have already 0167 %been deleted in 'sgd' analysis. 0168 [~, I]=ismember(originalGenes(details==1),model.genes); 0169 genesToModify=nchoosek(I,2); 0170 genes=geneMapping(genesToModify); 0171 grRatioMuts=zeros(1,numel(genes)); 0172 fluxes=sparse(numel(model.rxns),size(genesToModify,1)); 0173 for i=1:size(genesToModify,1) 0174 tempModel=removeGenes(model,genesToModify(i,:),false,false,false); 0175 0176 if strcmpi(analysisType,'fba') 0177 sol=solveLP(tempModel); 0178 else 0179 [fluxA, ~, flag]=qMOMA(tempModel,refModel); 0180 sol.x=fluxA; 0181 sol.stat=flag; 0182 end 0183 0184 if sol.stat==1 0185 fluxes(:,i)=sol.x; 0186 grRatioMuts(i)=sol.f/growthWT; 0187 end 0188 end 0189 end 0190 0191 %Map back to the old model 0192 [~, I]=ismember(model.rxns,originalModel.rxns); 0193 temp=fluxes; 0194 fluxes=sparse(numel(originalModel.rxns),size(temp,2)); 0195 fluxes(I,:)=temp; 0196 end 0197