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