0001 function [grRatioDble,grRateKO,grRateWT] = doubleGeneDeletion(model,method,geneList1,geneList2,verbFlag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
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
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
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
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
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
0117 if (mod(delCounter,1000) == 0)
0118 save doubleGeneDeletionTmp.mat grRateKO
0119 end
0120 [isInModel,geneID2] = ismember(geneList2{geneNo2},model.genes);
0121
0122 rxnInd = find(any(model.rxnGeneMat(:,[geneID1 geneID2]),2));
0123 if (~isempty(rxnInd))
0124
0125 x = true(nGenes,1);
0126 x(geneID1) = false;
0127 x(geneID2) = false;
0128 constrainRxn = false(length(rxnInd),1);
0129
0130 for rxnNo = 1:length(rxnInd)
0131 if (~eval(model.rules{rxnInd(rxnNo)}))
0132 constrainRxn(rxnNo) = true;
0133 end
0134 end
0135
0136 if (any(constrainRxn))
0137 constrRxnInd = rxnInd(constrainRxn);
0138 modelTmp = model;
0139 modelTmp.lb(constrRxnInd) = 0;
0140 modelTmp.ub(constrRxnInd) = 0;
0141
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
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
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;