0001 function [outModel geneLocalization transportStruct scores removedRxns]=predictLocalization(model,GSS,defaultCompartment,transportCost,maxTime,plotResults)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054 if nargin<4
0055 transportCost=0.5;
0056 end
0057
0058 if nargin<5
0059 maxTime=15;
0060 end
0061
0062 if nargin<6
0063 plotResults=false;
0064 end
0065
0066 I=ismember(defaultCompartment,GSS.compartments);
0067 if I==false
0068 throw(MException('','defaultCompartment not found in GSS'));
0069 end
0070
0071 if numel(model.comps)>1
0072 fprintf('WARNING: The model has several compartments. All compartments will be merged.\n');
0073 model=mergeCompartments(model,true,true);
0074 end
0075
0076
0077
0078
0079 model=expandModel(model);
0080
0081
0082
0083 removedRxns={};
0084 while 1
0085 irrevModel=convertToIrrev(model);
0086
0087 I=sum(irrevModel.S>0,2);
0088
0089
0090 if isfield(irrevModel,'unconstrained')
0091 I(irrevModel.unconstrained~=0)=2;
0092 end
0093 metsToDelete=false(numel(model.mets),1);
0094
0095
0096
0097 for i=1:numel(irrevModel.mets)
0098
0099
0100
0101 if I(i)<2
0102 if I(i)==1
0103
0104 [crap J]=find(irrevModel.S(i,:)>0);
0105
0106
0107
0108
0109 K=irrevModel.S(:,J)<0;
0110 check=find(K & I<=1);
0111
0112 for j=1:numel(check)
0113
0114 [crap L]=find(irrevModel.S(check(j),:)>0);
0115
0116 if ~isempty(L)
0117 rxn=irrevModel.rxns(J);
0118 rxnRev=irrevModel.rxns(L);
0119 if strcmp(strrep(rxn,'_REV',''),strrep(rxnRev,'_REV',''))
0120 metsToDelete(i)=true;
0121 end
0122 else
0123
0124
0125 continue;
0126 end
0127 end
0128 else
0129
0130 metsToDelete(i)=true;
0131 end
0132 end
0133 end
0134
0135 if any(metsToDelete)
0136
0137 [crap I]=find(model.S(metsToDelete,:));
0138 removedRxns=[removedRxns;model.rxns(I)];
0139 model=removeRxns(model,I,true,true);
0140 else
0141
0142 break;
0143 end
0144 end
0145
0146
0147
0148 I=find(sum(model.rxnGeneMat,2)==0);
0149 for i=1:numel(I)
0150 model.genes=[model.genes;['&&FAKE&&' num2str(i)]];
0151 if isfield(model,'geneShortNames')
0152 model.geneShortNames=[model.geneShortNames;{''}];
0153 end
0154 if isfield(model,'geneMiriams')
0155 model.geneMiriams=[model.geneMiriams;{[]}];
0156 end
0157 if isfield(model,'geneFrom')
0158 model.geneFrom=[model.geneFrom;{{'FAKE'}}];
0159 end
0160 model.rxnGeneMat(I(i),numel(model.genes))=1;
0161 model.grRules{I(i)}='';
0162 end
0163
0164
0165
0166
0167 I=setdiff(model.genes,GSS.genes);
0168 GSS.genes=[GSS.genes;I];
0169 GSS.scores=[GSS.scores;ones(numel(I),numel(GSS.compartments))*0.5];
0170
0171
0172
0173
0174
0175
0176
0177 genes=unique(model.grRules);
0178 nGenes=strrep(genes,'(','');
0179 nGenes=strrep(nGenes,')','');
0180
0181 complexes=setdiff(nGenes,model.genes);
0182 if isempty(complexes{1})
0183 complexes(1)=[];
0184 end
0185 cScores=zeros(numel(complexes),numel(GSS.compartments));
0186 for i=1:numel(complexes)
0187 genesInComplex=regexp(complexes{i},' and ','split');
0188
0189
0190 [I J]=ismember(genesInComplex,GSS.genes);
0191
0192 if any(I)
0193
0194 mScores=mean(GSS.scores(J(I),:));
0195
0196
0197
0198 mScores=(mScores.*sum(I)+(numel(genesInComplex)-sum(I))*0.5)/numel(genesInComplex);
0199 else
0200 fprintf(['WARNING: Could not parse grRule "' complexes{i} '". Assigning score 0.0 in all compartments\n']);
0201 mScores=ones(1,numel(genesInComplex))*0.5;
0202 end
0203 cScores(i,:)=mScores;
0204
0205
0206 model.genes=[model.genes;complexes{i}];
0207 if isfield(model,'geneMiriams')
0208 model.geneMiriams=[model.geneMiriams;{[]}];
0209 end
0210 if isfield(model,'geneShortNames')
0211 model.geneShortNames=[model.geneShortNames;{''}];
0212 end
0213 if isfield(model,'geneFrom')
0214 model.geneFrom=[model.geneFrom;{'COMPLEX'}];
0215 end
0216
0217
0218 I=ismember(model.grRules,['(' complexes{i} ')']);
0219
0220
0221 if ~isempty(I)
0222 model.rxnGeneMat(I,:)=0;
0223 model.rxnGeneMat(I,numel(model.genes))=1;
0224 end
0225 end
0226
0227
0228 GSS.genes=[GSS.genes;complexes];
0229 GSS.scores=[GSS.scores;cScores];
0230
0231
0232
0233 model=removeRxns(model,{},true,true);
0234
0235
0236
0237
0238
0239
0240 [crap I]=getExchangeRxns(model);
0241
0242
0243 J=1:numel(model.rxns);
0244 J(I)=[];
0245 model=permuteModel(model,[I;J'],'rxns');
0246
0247
0248 nER=numel(I);
0249
0250
0251 if isfield(model,'unconstrained')
0252 I=find(model.unconstrained);
0253 J=1:numel(model.mets);
0254 J(I)=[];
0255 model=permuteModel(model,[I;J'],'mets');
0256
0257 nEM=numel(I);
0258 else
0259 nEM=0;
0260 end
0261
0262
0263
0264 model.rxnGeneMat(1:nER,:)=0;
0265 model.grRules(1:nER)={''};
0266
0267
0268 model=removeRxns(model,{},true,true);
0269
0270
0271
0272
0273
0274 [crap J]=ismember(model.genes,GSS.genes);
0275 GSS.genes=model.genes;
0276 GSS.scores=GSS.scores(J,:);
0277
0278
0279
0280 [crap J]=ismember(defaultCompartment,GSS.compartments);
0281 reorder=1:numel(GSS.compartments);
0282 reorder(J)=[];
0283 reorder=[J reorder];
0284 GSS.scores=GSS.scores(:,reorder);
0285 GSS.compartments=GSS.compartments(reorder);
0286
0287
0288
0289
0290 oldS=model.S;
0291 model.S(model.S>0)=1;
0292 model.S(model.S<0)=-1;
0293
0294
0295
0296
0297 model.S(:,model.rev==1)=model.S(:,model.rev==1).*10;
0298
0299
0300
0301
0302 nRxns=numel(model.rxns)-nER;
0303 nMets=numel(model.mets)-nEM;
0304 nGenes=numel(model.genes);
0305 nComps=numel(GSS.compartments);
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317 s=repmat(eye(nMets)*0,1,nComps-1);
0318 s=[zeros(numel(model.mets)-nMets,size(s,2));s];
0319 S=[model.S sparse(numel(model.mets),nRxns*(nComps-1)) s];
0320 s=[sparse(nMets*(nComps-1),numel(model.rxns)+nRxns*(nComps-1)) eye(nMets*(nComps-1))*0];
0321 S=[S;s];
0322
0323
0324
0325 g2c=false(nGenes,nComps);
0326
0327 g2c(:,1)=true;
0328
0329
0330 tic;
0331 bestScore=-inf;
0332 bestS=[];
0333 bestg2c=[];
0334
0335
0336 plotScore=[];
0337 nTrans=[];
0338 totScore=[];
0339 minScore=sum(min(GSS.scores,[],2));
0340 maxScore=sum(max(GSS.scores,[],2));
0341 while toc<maxTime*60
0342
0343
0344
0345
0346
0347 [I J]=find(g2c);
0348 geneToMove=randsample(nGenes,1,true,max(GSS.scores(I,:),[],2)-GSS.scores(sub2ind(size(g2c),I,J))+0.1);
0349
0350
0351
0352
0353 toComp=randsample(nComps,1,true,GSS.scores(geneToMove,:)+0.2);
0354
0355
0356 if toComp==find(g2c(geneToMove,:))
0357 continue;
0358 end
0359
0360
0361 [newS newg2c]=moveGene(S,model,g2c,geneToMove,toComp,nRxns,nMets);
0362
0363
0364
0365
0366 wasConnected=false;
0367 for j=1:5
0368
0369 unconnected=findUnconnected(newS,nEM);
0370
0371
0372 if any(unconnected)
0373
0374
0375
0376 [geneIndex moveTo deltaConnected deltaScore]=selectGenes(newS,nEM,nMets,nER,nRxns,model,unconnected,g2c,GSS);
0377
0378
0379
0380
0381
0382 [score I]=max(1.5*deltaScore+deltaConnected);
0383
0384
0385
0386 hasToAddTransport=true;
0387 if ~isempty(deltaConnected)
0388 if score>0
0389 hasToAddTransport=false;
0390 end
0391 end
0392
0393
0394
0395 if hasToAddTransport==false;
0396 [newS newg2c]=moveGene(newS,model,g2c,geneIndex(I),moveTo(I),nRxns,nMets);
0397 else
0398
0399
0400 transMet=unconnected(randsample(numel(unconnected),1));
0401
0402
0403 comps=ceil((transMet-nEM)/((size(S,1)-nEM)/nComps));
0404
0405
0406
0407 dcIndex=transMet-(comps-1)*nMets;
0408
0409
0410
0411 allIndexes=dcIndex;
0412 for k=1:nComps-1
0413 allIndexes=[allIndexes;dcIndex+nMets*k];
0414 end
0415
0416
0417
0418 I=sum(newS(allIndexes,:)~=0,2)>0;
0419
0420
0421
0422
0423 connectedUsed=setdiff(allIndexes(I),unconnected);
0424
0425
0426
0427 if isempty(connectedUsed)
0428 break;
0429 end
0430
0431
0432
0433 if transMet==dcIndex
0434 newS=addTransport(newS,nRxns,nER,nMets,nEM,nComps,transMet,connectedUsed(randsample(numel(connectedUsed),1)));
0435 else
0436
0437
0438 I=connectedUsed(connectedUsed<(nMets+nEM));
0439 if any(I)
0440 newS=addTransport(newS,nRxns,nER,nMets,nEM,nComps,transMet,I(randsample(numel(I),1)));
0441 else
0442
0443
0444
0445 break;
0446 end
0447 end
0448 end
0449 else
0450 wasConnected=true;
0451 break;
0452 end
0453 end
0454
0455
0456
0457 if wasConnected==true
0458
0459 activeTransport=find(sum(newS(:,nER+nRxns*nComps+1:end),2));
0460
0461
0462 unconnected=findUnconnected(newS(:,1:nER+nRxns*nComps),nEM);
0463
0464
0465 I=setdiff(activeTransport,unconnected);
0466
0467
0468
0469 newS(:,find(sum(newS(I,nER+nRxns*nComps+1:end))==4)+nER+nRxns*nComps)=0;
0470
0471
0472 [score geneScore trCost]=scoreModel(newS,newg2c,GSS,transportCost);
0473
0474
0475 if score>bestScore
0476 bestScore=score;
0477 bestS=newS;
0478 bestg2c=newg2c;
0479 end
0480
0481
0482 if score>=bestScore
0483 plotScore=[plotScore;geneScore];
0484 nTrans=[nTrans;trCost];
0485 totScore=[totScore;score];
0486 S=newS;
0487 g2c=newg2c;
0488
0489 if plotResults==true
0490 subplot(3,2,1);
0491 spy(S);
0492 subplot(3,2,2);
0493 plot(plotScore,'r');
0494 xlabel('Gene score');
0495 subplot(3,2,3);
0496 plot((plotScore-minScore)/(maxScore-minScore),'r');
0497 xlabel('Gene score relative to predictions');
0498 subplot(3,2,4);
0499 plot(nTrans,'g');
0500 xlabel('Transport cost');
0501 subplot(3,2,5);
0502 plot(totScore,'b');
0503 xlabel('Total score');
0504 subplot(3,2,6);
0505 pause(0.2);
0506 end
0507 end
0508 end
0509 end
0510 scores.totScore=score;
0511 scores.geneScore=geneScore;
0512 scores.transCost=trCost;
0513
0514
0515 [I J]=find(bestS(nEM+1:nEM+nMets,end-nMets*(nComps-1)+1:end));
0516 J=ceil(J/nMets+1);
0517 transportStruct.mets=model.metNames(I+nEM);
0518 transportStruct.toComp=GSS.compartments(J);
0519
0520 [I J]=find(bestg2c);
0521 geneLocalization.genes=GSS.genes(I);
0522 geneLocalization.comps=GSS.compartments(J);
0523
0524
0525 [crap I]=sort(geneLocalization.genes);
0526 geneLocalization.genes=geneLocalization.genes(I);
0527 geneLocalization.comps=geneLocalization.comps(I);
0528
0529
0530 I=strmatch('&&FAKE&&',geneLocalization.genes);
0531 geneLocalization.genes(I)=[];
0532 geneLocalization.comps(I)=[];
0533
0534
0535
0536
0537
0538
0539
0540 outModel=model;
0541 outModel.S=oldS;
0542
0543
0544 copyPart=outModel.S(nEM+1:end,nER+1:end);
0545
0546
0547 copyRxnGeneMat=outModel.rxnGeneMat(nER+1:end,:);
0548 outModel.rxnGeneMat=[outModel.rxnGeneMat;repmat(copyRxnGeneMat,nComps-1,1)];
0549
0550
0551
0552
0553 nStartComps=numel(outModel.comps);
0554 if nStartComps==1
0555 outModel.comps={'1'};
0556 outModel.compNames=GSS.compartments(1);
0557 else
0558 if model.metComps(1)==1
0559 outModel.compNames(1)=GSS.compartments(1);
0560 else
0561 outModel.compNames(2)=GSS.compartments(1);
0562 end
0563 end
0564 outModel.compNames=[outModel.compNames;GSS.compartments(2:end)];
0565
0566
0567 for i=1:numel(GSS.compartments)-1
0568 outModel.comps=[outModel.comps;num2str(numel(outModel.comps)+1)];
0569 end
0570
0571 outModel.compOutside=cell(numel(outModel.comps),1);
0572 outModel.compOutside(:)={''};
0573
0574 for i=1:nComps-1
0575 outModel.S=[outModel.S sparse(size(outModel.S,1),nRxns)];
0576 outModel.S=[outModel.S;[sparse(nMets,nRxns*i+nER) copyPart]];
0577 outModel.rxns=[outModel.rxns;strcat(outModel.rxns(nER+1:nER+nRxns),'_',GSS.compartments{i+1})];
0578 outModel.rxnNames=[outModel.rxnNames;strcat(outModel.rxnNames(nER+1:nER+nRxns),' (',GSS.compartments{i+1},')')];
0579 outModel.lb=[outModel.lb;outModel.lb(nER+1:nER+nRxns)];
0580 outModel.ub=[outModel.ub;outModel.ub(nER+1:nER+nRxns)];
0581 outModel.rev=[outModel.rev;outModel.rev(nER+1:nER+nRxns)];
0582 outModel.c=[outModel.c;outModel.c(nER+1:nER+nRxns)];
0583 if isfield(outModel,'grRules')
0584 outModel.grRules=[outModel.grRules;outModel.grRules(nER+1:nER+nRxns)];
0585 end
0586 if isfield(outModel,'subSystems')
0587 outModel.subSystems=[outModel.subSystems;outModel.subSystems(nER+1:nER+nRxns)];
0588 end
0589 if isfield(outModel,'eccodes')
0590 outModel.eccodes=[outModel.eccodes;outModel.eccodes(nER+1:nER+nRxns)];
0591 end
0592 if isfield(outModel,'rxnFrom')
0593 outModel.rxnFrom=[outModel.rxnFrom;outModel.rxnFrom(nER+1:nER+nRxns)];
0594 end
0595 outModel.mets=[outModel.mets;strcat(outModel.mets(nEM+1:nEM+nMets),'_',GSS.compartments{i+1})];
0596 outModel.metNames=[outModel.metNames;outModel.metNames(nEM+1:nEM+nMets)];
0597 outModel.b=[outModel.b;outModel.b(nEM+1:nEM+nMets,:)];
0598 I=ones(nMets,1)*nStartComps+i;
0599 outModel.metComps=[outModel.metComps;I];
0600 if isfield(outModel,'inchis')
0601 outModel.inchis=[outModel.inchis;outModel.inchis(nEM+1:nEM+nMets)];
0602 end
0603 if isfield(outModel,'unconstrained')
0604 outModel.unconstrained=[outModel.unconstrained;outModel.unconstrained(nEM+1:nEM+nMets)];
0605 end
0606 if isfield(outModel,'metMiriams')
0607 outModel.metMiriams=[outModel.metMiriams;outModel.metMiriams(nEM+1:nEM+nMets)];
0608 end
0609 if isfield(outModel,'metFormulas')
0610 outModel.metFormulas=[outModel.metFormulas;outModel.metFormulas(nEM+1:nEM+nMets)];
0611 end
0612 if isfield(outModel,'metFrom')
0613 outModel.metFrom=[outModel.metFrom;outModel.metFrom(nEM+1:nEM+nMets)];
0614 end
0615 end
0616
0617
0618 transS=bestS(:,numel(outModel.rxns)+1:end);
0619 J=sum(transS)>0;
0620
0621
0622
0623 transS(transS~=0)=1;
0624 transS(1:nEM+nMets,:)=transS(1:nEM+nMets,:)*-1;
0625 I=find(sum(transS>0,2));
0626 nTransRxns=numel(I);
0627 outModel.S=[outModel.S transS(:,J)];
0628 filler=ones(nTransRxns,1);
0629 outModel.lb=[outModel.lb;filler*-1000];
0630 outModel.ub=[outModel.ub;filler*1000];
0631 outModel.rev=[outModel.rev;filler];
0632 outModel.c=[outModel.c;filler*0];
0633 outModel.rxnGeneMat=[outModel.rxnGeneMat;sparse(nTransRxns,numel(outModel.genes))];
0634
0635 for i=1:numel(I)
0636 outModel.rxns=[outModel.rxns;strcat('transport',num2str(i))];
0637 outModel.rxnNames=[outModel.rxnNames;['Transport of ',outModel.metNames{I(i)}]];
0638 if isfield(outModel,'grRules')
0639 outModel.grRules=[outModel.grRules;{''}];
0640 end
0641 if isfield(outModel,'subSystems')
0642 outModel.subSystems=[outModel.subSystems;'Inferred transport reactions'];
0643 end
0644 if isfield(outModel,'eccodes')
0645 outModel.eccodes=[outModel.eccodes;{''}];
0646 end
0647 if isfield(outModel,'rxnFrom')
0648 outModel.rxnFrom=[outModel.rxnFrom;{''}];
0649 end
0650 end
0651
0652
0653
0654 [I J]=find(bestS(:,1:nER+nComps*nRxns));
0655 K=true(numel(outModel.rxns),1);
0656 K(J)=false;
0657 K(end-nTransRxns+1:end)=false;
0658 outModel=removeRxns(outModel,K,true);
0659
0660
0661 I=strmatch('&&FAKE&&',outModel.genes);
0662 outModel.genes(I)=[];
0663 if isfield(outModel,'geneMiriams')
0664 outModel.geneMiriams(I)=[];
0665 end
0666 if isfield(outModel,'geneShortNames')
0667 outModel.geneShortNames(I)=[];
0668 end
0669 outModel.rxnGeneMat(:,I)=[];
0670 end
0671
0672
0673 function [S g2c]=moveGene(S,model,g2c,geneToMove,toComp,nRxns,nMets)
0674
0675 currentComp=find(g2c(geneToMove,:));
0676 g2c(geneToMove,:)=false;
0677 g2c(geneToMove,toComp)=true;
0678
0679
0680 [I crap]=find(model.rxnGeneMat(:,geneToMove));
0681
0682
0683 oldRxns=I+(currentComp-1)*nRxns;
0684
0685
0686 newRxns=I+(toComp-1)*nRxns;
0687
0688
0689
0690 metChange=nMets*(toComp-currentComp);
0691
0692
0693 [I J K]=find(S(:,oldRxns));
0694 I=I+metChange;
0695
0696
0697 S(:,oldRxns)=0;
0698 S(sub2ind(size(S),I,newRxns(J)))=K;
0699 end
0700
0701
0702
0703
0704
0705
0706
0707 function unconnected=findUnconnected(S,nEM,metsToCheck)
0708 if nargin>2
0709
0710
0711 I=false(size(S,1),1);
0712 I(1:nEM)=true;
0713 I(metsToCheck)=true;
0714 S=S(I,:);
0715 end
0716
0717 em=false(size(S,1),1);
0718 em(1:nEM)=true;
0719
0720
0721 I=sum(S>2,1) | sum(S>2,1);
0722 revS=S;
0723 revS(:,I)=revS(:,I)*-1;
0724
0725
0726
0727
0728 connected=sum(S>0,2)>1 | em | sum(S~=0,2)==0 | sum(S(:,~I)>0,2)>0 | sum(S(:,I)~=0,2)>1;
0729
0730
0731 unconnected=sum(S>0 | revS>0,2)==0 & connected==false;
0732
0733
0734 maybeUnconnected=~connected & ~unconnected;
0735
0736
0737
0738
0739
0740
0741
0742
0743 deadRxns=any(S(maybeUnconnected,:)>0) & any(S(maybeUnconnected,:)<0);
0744
0745
0746 problematic=any(S(:,deadRxns)~=0,2);
0747
0748
0749
0750 unconnected(problematic & maybeUnconnected)=true;
0751
0752
0753 if nargin>2
0754 unconnected=metsToCheck(unconnected(nEM+1:end));
0755 else
0756 unconnected=find(unconnected);
0757 end
0758 end
0759
0760
0761
0762
0763
0764
0765 function [geneIndex moveTo deltaConnected deltaScore]=selectGenes(S,nEM,nMets,nER,nRxns,model,unconnected,g2c,GSS)
0766
0767 moveTo=zeros(numel(model.genes),1);
0768 deltaConnected=zeros(numel(model.genes),1);
0769
0770
0771 nComps=size(g2c,2);
0772 comps=ceil((unconnected-nEM)/((size(S,1)-nEM)/nComps));
0773
0774
0775
0776 dcIndexes=unconnected-(comps-1)*nMets;
0777
0778
0779 allIndexes=dcIndexes;
0780 for i=1:nComps-1
0781 allIndexes=[allIndexes;dcIndexes+nMets*i];
0782 end
0783
0784
0785 I=sum(S>2,1) | sum(S>2,1);
0786 revS=S;
0787 revS(:,I)=revS(:,I)*-1;
0788
0789
0790
0791 newMets=setdiff(allIndexes,unconnected);
0792 potential=find(sum(S(newMets,:)>0) | sum(revS(newMets,:)>0));
0793 potential(potential<=nER | potential>nER+nRxns*nComps)=[];
0794
0795
0796 rxnComps=ceil((potential-nER)/(nRxns));
0797
0798
0799
0800 dcRxnIndexes=potential-(rxnComps-1)*nRxns;
0801
0802
0803 genes=find(sum(model.rxnGeneMat(dcRxnIndexes,:)>0,1));
0804
0805
0806
0807
0808
0809
0810
0811
0812 [crap J]=find(g2c(genes,:));
0813
0814
0815
0816 geneScores=GSS.scores(sub2ind(size(g2c),genes(:),J));
0817 modGeneScores=1.1-geneScores;
0818 if numel(genes)>25
0819 rGenes=genes(randsample(numel(genes),min(numel(genes),25),true,modGeneScores));
0820
0821
0822 rGenes=unique(rGenes);
0823
0824
0825 [crap I]=ismember(rGenes,genes);
0826 geneScores=geneScores(I);
0827 genes=rGenes;
0828 end
0829 for i=1:numel(genes)
0830
0831
0832
0833
0834 rxns=find(model.rxnGeneMat(:,genes(i)));
0835
0836
0837 mets=find(sum(model.S(:,rxns)~=0,2)>0);
0838
0839
0840 allIndexes=mets;
0841 for j=1:nComps-1
0842 allIndexes=[allIndexes;mets+nMets*j];
0843 end
0844
0845
0846
0847
0848
0849
0850 [I crap]=find(model.S(:,rxns));
0851 moveToComps=unique(comps(ismember(dcIndexes,I)));
0852
0853
0854 bestMove=-inf;
0855 bestComp=[];
0856 for j=1:numel(moveToComps)
0857 newS=moveGene(S,model,g2c,genes(i),moveToComps(j),nRxns,nMets);
0858
0859
0860
0861 dConnected=numel(unconnected)-numel(findUnconnected(newS,nEM,[allIndexes;unconnected]));
0862 if dConnected>bestMove
0863 bestMove=dConnected;
0864 bestComp=moveToComps(j);
0865 end
0866 end
0867
0868
0869
0870 moveTo(genes(i))=bestComp;
0871 deltaConnected(genes(i))=bestMove;
0872 end
0873
0874
0875 geneIndex=genes(:);
0876 moveTo=moveTo(geneIndex);
0877 deltaConnected=deltaConnected(geneIndex);
0878 deltaScore=GSS.scores(sub2ind(size(g2c),geneIndex(:),moveTo))-geneScores;
0879 end
0880
0881
0882
0883
0884
0885 function S=addTransport(S,nRxns,nER,nMets,nEM,nComps,metA,metB)
0886 mets=[metA;metB];
0887
0888 comps=ceil((mets-nEM)/((size(S,1)-nEM)/nComps));
0889
0890 if sum(comps==1)~=1
0891 throw(MException('','Tried to create a transport reaction from a non-default compartment'));
0892 end
0893
0894
0895 rIndex=(nER+nRxns*nComps)+mets(comps~=1)-nEM-nMets;
0896
0897 S(mets,rIndex)=2;
0898 end
0899
0900
0901
0902 function [score geneScore transportCost]=scoreModel(S,g2c,GSS,transportCost)
0903 [I J]=find(g2c);
0904 geneScore=sum(GSS.scores(sub2ind(size(g2c),I,J)));
0905 transportCost=sum(sum(S==2,1)>0)*transportCost;
0906 score=geneScore-transportCost;
0907 end