doubleGeneDeletion

PURPOSE ^

doubleGeneDeletion Performs double gene deletion analysis using FBA, MOMA,

SYNOPSIS ^

function [grRatioDble,grRateKO,grRateWT] = doubleGeneDeletion(model,method,geneList1,geneList2,verbFlag)

DESCRIPTION ^

doubleGeneDeletion Performs double gene deletion analysis using FBA, MOMA,
or linear MOMA

 [grRatioDble,grRateKO,grRateWT] =
     doubleGeneDeletion(model,method,geneList1,geneList2,verbFlag)

INPUT
 model         COBRA model structure

OPTIONAL INPUTS
 method        Either 'FBA' (default) 'MOMA', or 'lMOMA'
 geneList1     List of query genes to be deleted (default = all genes)
 geneList2     List of target genes to be deleted (default = geneList1)
 verbFlag      Verbose output (default = false)

OUTPUTS
 grRatioDble   Computed growth rate ratio between double deletion strain 
               and wild type
 grRateKO      Double deletion strain growth rates (1/h)
 grRateWT      Wild type growth rate (1/h)

 Markus Herrgard 8/8/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [grRatioDble,grRateKO,grRateWT] = doubleGeneDeletion(model,method,geneList1,geneList2,verbFlag)
0002 %doubleGeneDeletion Performs double gene deletion analysis using FBA, MOMA,
0003 %or linear MOMA
0004 %
0005 % [grRatioDble,grRateKO,grRateWT] =
0006 %     doubleGeneDeletion(model,method,geneList1,geneList2,verbFlag)
0007 %
0008 %INPUT
0009 % model         COBRA model structure
0010 %
0011 %OPTIONAL INPUTS
0012 % method        Either 'FBA' (default) 'MOMA', or 'lMOMA'
0013 % geneList1     List of query genes to be deleted (default = all genes)
0014 % geneList2     List of target genes to be deleted (default = geneList1)
0015 % verbFlag      Verbose output (default = false)
0016 %
0017 %OUTPUTS
0018 % grRatioDble   Computed growth rate ratio between double deletion strain
0019 %               and wild type
0020 % grRateKO      Double deletion strain growth rates (1/h)
0021 % grRateWT      Wild type growth rate (1/h)
0022 %
0023 % Markus Herrgard 8/8/06
0024 
0025 if (nargin < 2)
0026     method = 'FBA';
0027 end
0028 if (nargin < 3)
0029     geneList1 = model.genes;
0030     differentSetsFlag = false;
0031 else
0032     if (isempty(geneList1))
0033         geneList1 = model.genes;
0034     end
0035 end
0036 
0037 if (nargin < 4)
0038     geneList2 = geneList1;
0039     differentSetsFlag = false;
0040 else
0041     if (isempty(geneList2))
0042         geneList2 = geneList1;
0043         differentSetsFlag = false;
0044     else
0045         differentSetsFlag = true;
0046     end
0047 end
0048 if (nargin < 5)
0049     verbFlag = false;
0050 end
0051 
0052 nGenes = length(model.genes);
0053 
0054 % Run single gene deletions first to figure out lethal genes
0055 fprintf('Single deletion analysis to remove lethal genes\n');
0056 [singleRatio1,singleRate1,grRateWT] = singleGeneDeletion(model,method,geneList1,verbFlag);
0057 singleLethal1 = (singleRatio1 < 1e-9);
0058 geneListOrig1 = geneList1;
0059 geneListOrig2 = geneList1;
0060 geneList1 = geneList1(~singleLethal1);
0061 singleRate = singleRate1(~singleLethal1);
0062 [tmp,listMap1] = ismember(geneListOrig1,geneList1);
0063 fprintf('%d non-lethal genes\n',length(geneList1));
0064 
0065 % Repeat the analysis for the second set of genes
0066 if (differentSetsFlag)
0067     fprintf('Single deletion analysis to remove lethal genes from gene set 2\n');
0068     [singleRatio2,singleRate2,grRateWT] = singleGeneDeletion(model,method,geneList2,verbFlag);
0069     singleLethal2 = (singleRatio2 < 1e-9);
0070     geneListOrig2 = geneList2;
0071     geneList2 = geneList2(~singleLethal2);
0072     [tmp,listMap2] = ismember(geneListOrig2,geneList2);
0073     fprintf('%d non-lethal genes\n',length(geneList2));
0074 else
0075     geneList2 = geneList1;
0076     listMap2 = listMap1;
0077 end
0078 
0079 nDelGenes1 = length(geneList1);
0080 nDelGenes2 = length(geneList2);
0081 
0082 grRateKO = ones(nDelGenes1,nDelGenes2)*grRateWT;
0083 grRatioDble = ones(nDelGenes1,nDelGenes2);
0084 
0085 if (differentSetsFlag)
0086     nTotalPairs = nDelGenes1*nDelGenes2;
0087 else
0088     nTotalPairs = nDelGenes1*(nDelGenes1-1)/2;
0089 end
0090 
0091 % Run double deletion analysis
0092 delCounter = 0;
0093 fprintf('Double gene deletion analysis\n');
0094 fprintf('Total of %d pairs to analyze\n',nTotalPairs);
0095 h = waitbar(0,'Double gene deletion analysis in progress ...');
0096 t = cputime;
0097 fprintf('Perc complete\tCPU time\n');
0098 for geneNo1 = 1:nDelGenes1
0099    
0100     % Find gene index
0101     [isInModel,geneID1] = ismember(geneList1{geneNo1},model.genes);
0102     if (~differentSetsFlag)
0103         grRateKO(geneNo1,geneNo1) = singleRate(geneNo1);
0104         initID = geneNo1+1;
0105     else
0106         initID = 1;
0107     end
0108     for geneNo2 = initID:nDelGenes2
0109         delCounter = delCounter + 1;
0110         if (mod(delCounter,10) == 0)
0111             waitbar(delCounter/nTotalPairs,h);
0112         end
0113         if (mod(delCounter,100) == 0)
0114             fprintf('%5.2f\t%8.1f\n',100*delCounter/nTotalPairs,cputime-t);
0115         end
0116         % Save results every 1000 steps
0117         if (mod(delCounter,1000) == 0)
0118             save doubleGeneDeletionTmp.mat grRateKO
0119         end
0120         [isInModel,geneID2] = ismember(geneList2{geneNo2},model.genes);
0121         % Find rxns associated with this gene
0122         rxnInd = find(any(model.rxnGeneMat(:,[geneID1 geneID2]),2));
0123         if (~isempty(rxnInd))
0124             % Initialize the state of all genes
0125             x = true(nGenes,1);
0126             x(geneID1) = false;
0127             x(geneID2) = false;
0128             constrainRxn = false(length(rxnInd),1);
0129             % Figure out if any of the reaction states is changed
0130             for rxnNo = 1:length(rxnInd)
0131                 if (~eval(model.rules{rxnInd(rxnNo)}))
0132                     constrainRxn(rxnNo) = true;
0133                 end
0134             end
0135             % Use FBA/MOMA/lMOMA to calculate deletion strain growth rate
0136             if (any(constrainRxn))
0137                 constrRxnInd = rxnInd(constrainRxn);
0138                 modelTmp = model;
0139                 modelTmp.lb(constrRxnInd) = 0;
0140                 modelTmp.ub(constrRxnInd) = 0;
0141                 % Get double deletion growth rate
0142                 switch method
0143                     case 'lMOMA'
0144                         solKO = linearMOMA(model,modelTmp,'max');
0145                     case 'MOMA'
0146                         solKO = MOMA(model,modelTmp,'max',false,true);
0147                     otherwise
0148                         solKO = optimizeCbModel(modelTmp,'max');
0149                 end
0150                 %solKO = optimizeCbModel(modelTmp,'max');
0151                 if (solKO.stat > 0)
0152                     grRateKO(geneNo1,geneNo2) = solKO.f;
0153                     grRateKO(geneNo2,geneNo1) = solKO.f; 
0154                 else
0155                     grRateKO(geneNo1,geneNo2) = 0;
0156                     grRateKO(geneNo2,geneNo1) = 0;
0157                 end
0158             end
0159         end
0160         if (verbFlag)
0161             fprintf('%4d\t%4.1f\t%10s\t%10s\t%9.3f\t%9.3f\n',delCounter,100*delCounter/nTotalPairs,geneList1{geneNo1},...
0162                 geneList2{geneNo2},grRateKO(geneNo1,geneNo2),grRateKO(geneNo1,geneNo2)/grRateWT*100);
0163         end
0164         if (differentSetsFlag)
0165             grRateKO(geneNo2,geneNo1) = grRateKO(geneNo1,geneNo2);
0166         end
0167     end
0168 end
0169 if ( regexp( version, 'R20') )
0170         close(h);
0171 end
0172 
0173 % Reconstruct the entire matrix
0174 for i = 1:length(geneListOrig1)
0175     for j = 1:length(geneListOrig2)
0176         if (listMap1(i) > 0 & listMap2(j) > 0)
0177             allGrRateKO(i,j) = grRateKO(listMap1(i),listMap2(j));
0178         else
0179             allGrRateKO(i,j) = 0;
0180         end
0181     end
0182 end
0183 
0184 grRatioDble = allGrRateKO/grRateWT;
0185 
0186 grRateKO = allGrRateKO;

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003