diff --git a/core/replaceMets.m b/core/replaceMets.m index 9301d8628..9cfa0804c 100755 --- a/core/replaceMets.m +++ b/core/replaceMets.m @@ -20,7 +20,7 @@ % Output: % model model structure with selected metabolites replaced % removedRxns identifiers of duplicate reactions that were removed -% idxDuplRxns index of removedRxns in original model +% idxDuplRxns indices of removedRxns in the input model % % Note: This function is useful when the model contains both 'oxygen' and % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, @@ -59,10 +59,11 @@ error('The to-be-replaced metabolite cannot be found in the model.'); end -rxnsWithMet = find(model.S(metIdx,:)); +[~, rxnsWithMet] = find(model.S(metIdx,:)); +rxnsWithMet = unique(rxnsWithMet); if verbose==true - fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') - fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) + fprintf('\n\nThe following reactions involve the to-be-replaced metabolite:\n') + fprintf('%s\n',model.rxns{rxnsWithMet}) fprintf('\n') end @@ -88,8 +89,11 @@ idxDelete=[]; if identifiers + if any(ismember(metIdx,repIdx)) + error('The to-be-replaced metabolite and the replacement metabolite are the same.'); + end originalStoch = model.S(metIdx,rxnsWithMet); - model.S(repIdx,rxnsWithMet) = originalStoch; + model.S(repIdx,rxnsWithMet) = model.S(repIdx,rxnsWithMet) + originalStoch; model.S(metIdx,rxnsWithMet) = 0; idxDelete = metIdx; else @@ -99,54 +103,38 @@ % the to-be-replace metabolite ID deleted. % Build list of metaboliteName[compartment] - metCompsN =cellstr(num2str(model.metComps)); - map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); - metCompsN = map.values(metCompsN); - metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); + comps = model.comps(:); + metCompsN = strcat(lower(model.metNames(:)),'[',comps(model.metComps(:)),']'); - for i = 1:length(repIdx) - metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); - if length(metCompsNidx)>1 - for j = 2:length(metCompsNidx) - model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); - idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete - end + % Within each metaboliteName[compartment] group, the first replacement + % metabolite is kept and the other members are merged into it. + groups = unique(metCompsN(repIdx),'stable'); + for i = 1:numel(groups) + metCompsNidx = find(strcmp(groups{i}, metCompsN)); + keepIdx = metCompsNidx(find(ismember(metCompsNidx,repIdx),1)); + others = metCompsNidx(metCompsNidx ~= keepIdx); + if ~isempty(others) + model.S(keepIdx,:) = model.S(keepIdx,:) + sum(model.S(others,:),1); + idxDelete = [idxDelete; others]; % Make list of metabolite IDs to delete end end end - if ~isempty(idxDelete) - model.S(idxDelete,:) =[]; - model.mets(idxDelete) = []; - model.metNames(idxDelete) = []; - model.metComps(idxDelete) = []; - model.b(idxDelete) = []; - if isfield(model,'metFormulas') - model.metFormulas(idxDelete) = []; - end - if isfield(model,'unconstrained') - model.unconstrained(idxDelete) = []; - end - if isfield(model,'metMiriams') - model.metMiriams(idxDelete) = []; - end - if isfield(model,'metCharges') - model.metCharges(idxDelete) = []; - end - if isfield(model,'metDeltaG') - model.metDeltaG(idxDelete) = []; - end - if isfield(model,'inchis') - model.inchis(idxDelete) = []; - end - if isfield(model,'metSmiles') - model.metSmiles(idxDelete) = []; - end - if isfield(model,'metFrom') - model.metFrom(idxDelete) = []; - end + repMets = model.mets(repIdx); + model = removeMets(model,idxDelete); + % Indices of the replacement metabolites shift when metabolites with a + % lower index are deleted. + [~, repIdx] = ismember(repMets,model.mets); + repIdx(repIdx==0) = []; end % This could now have created duplicate reactions. Contract model. -model=contractModel(model,[],repIdx); +if isempty(idxDelete) + removedRxns = cell(0,1); + idxDuplRxns = zeros(0,1); +else + originalRxns = model.rxns; + [model, removedRxns] = contractModel(model,[],repIdx); + [~, idxDuplRxns] = ismember(removedRxns,originalRxns); +end end diff --git a/doc/core/removeMets.html b/doc/core/removeMets.html index 4a3d2d269..2df79cdf2 100644 --- a/doc/core/removeMets.html +++ b/doc/core/removeMets.html @@ -58,7 +58,7 @@

CROSS-REFERENCE INFORMATION ^ <li><a href=convertCharArray convertCharArray
  • getIndexes getIndexes
  • removeReactions removeReactions
    • This function is called by: +
  • mergeCompartments mergeCompartments
  • removeReactions removeReactions
  • replaceMets replaceMets
  • simplifyModel simplifyModel
    • diff --git a/doc/core/replaceMets.html b/doc/core/replaceMets.html index 6f4319513..db6d7f9b3 100644 --- a/doc/core/replaceMets.html +++ b/doc/core/replaceMets.html @@ -48,7 +48,7 @@

    DESCRIPTION ^DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

    This function calls: +
  • contractModel contractModel
  • removeMets removeMets
    • This function is called by: @@ -91,7 +91,7 @@

    SOURCE CODE ^% Output: 0021 % model model structure with selected metabolites replaced 0022 % removedRxns identifiers of duplicate reactions that were removed -0023 % idxDuplRxns index of removedRxns in original model +0023 % idxDuplRxns indices of removedRxns in the input model 0024 % 0025 % Note: This function is useful when the model contains both 'oxygen' and 0026 % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, @@ -130,97 +130,85 @@

    SOURCE CODE ^'The to-be-replaced metabolite cannot be found in the model.'); 0060 end 0061 -0062 rxnsWithMet = find(model.S(metIdx,:)); -0063 if verbose==true -0064 fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') -0065 fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) -0066 fprintf('\n') -0067 end -0068 -0069 model.metNames(metIdx) = model.metNames(repIdx(1)); -0070 if isfield(model,'metFormulas') -0071 model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); -0072 end -0073 if isfield(model,'metMiriams') -0074 model.metMiriams(metIdx) = model.metMiriams(repIdx(1)); -0075 end -0076 if isfield(model,'metCharges') -0077 model.metCharges(metIdx) = model.metCharges(repIdx(1)); -0078 end -0079 if isfield(model,'metDeltaG') -0080 model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1)); -0081 end -0082 if isfield(model,'inchis') -0083 model.inchis(metIdx) = model.inchis(repIdx(1)); -0084 end -0085 if isfield(model,'metSmiles') -0086 model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); -0087 end -0088 -0089 idxDelete=[]; -0090 if identifiers -0091 originalStoch = model.S(metIdx,rxnsWithMet); -0092 model.S(repIdx,rxnsWithMet) = originalStoch; -0093 model.S(metIdx,rxnsWithMet) = 0; -0094 idxDelete = metIdx; -0095 else -0096 % Run through replacement metabolites and their compartments. If any of the -0097 % to-be-replaced metabolites is already present (checked by -0098 % metaboliteName[compartment], then the replacement metabolite is kept and -0099 % the to-be-replace metabolite ID deleted. -0100 -0101 % Build list of metaboliteName[compartment] -0102 metCompsN =cellstr(num2str(model.metComps)); -0103 map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); -0104 metCompsN = map.values(metCompsN); -0105 metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); -0106 -0107 for i = 1:length(repIdx) -0108 metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); -0109 if length(metCompsNidx)>1 -0110 for j = 2:length(metCompsNidx) -0111 model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); -0112 idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete -0113 end -0114 end -0115 end -0116 end -0117 -0118 if ~isempty(idxDelete) -0119 model.S(idxDelete,:) =[]; -0120 model.mets(idxDelete) = []; -0121 model.metNames(idxDelete) = []; -0122 model.metComps(idxDelete) = []; -0123 model.b(idxDelete) = []; -0124 if isfield(model,'metFormulas') -0125 model.metFormulas(idxDelete) = []; -0126 end -0127 if isfield(model,'unconstrained') -0128 model.unconstrained(idxDelete) = []; -0129 end -0130 if isfield(model,'metMiriams') -0131 model.metMiriams(idxDelete) = []; -0132 end -0133 if isfield(model,'metCharges') -0134 model.metCharges(idxDelete) = []; -0135 end -0136 if isfield(model,'metDeltaG') -0137 model.metDeltaG(idxDelete) = []; -0138 end -0139 if isfield(model,'inchis') -0140 model.inchis(idxDelete) = []; -0141 end -0142 if isfield(model,'metSmiles') -0143 model.metSmiles(idxDelete) = []; -0144 end -0145 if isfield(model,'metFrom') -0146 model.metFrom(idxDelete) = []; -0147 end -0148 end -0149 -0150 % This could now have created duplicate reactions. Contract model. -0151 model=contractModel(model,[],repIdx); -0152 end +0062 [~, rxnsWithMet] = find(model.S(metIdx,:)); +0063 rxnsWithMet = unique(rxnsWithMet); +0064 if verbose==true +0065 fprintf('\n\nThe following reactions involve the to-be-replaced metabolite:\n') +0066 fprintf('%s\n',model.rxns{rxnsWithMet}) +0067 fprintf('\n') +0068 end +0069 +0070 model.metNames(metIdx) = model.metNames(repIdx(1)); +0071 if isfield(model,'metFormulas') +0072 model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); +0073 end +0074 if isfield(model,'metMiriams') +0075 model.metMiriams(metIdx) = model.metMiriams(repIdx(1)); +0076 end +0077 if isfield(model,'metCharges') +0078 model.metCharges(metIdx) = model.metCharges(repIdx(1)); +0079 end +0080 if isfield(model,'metDeltaG') +0081 model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1)); +0082 end +0083 if isfield(model,'inchis') +0084 model.inchis(metIdx) = model.inchis(repIdx(1)); +0085 end +0086 if isfield(model,'metSmiles') +0087 model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); +0088 end +0089 +0090 idxDelete=[]; +0091 if identifiers +0092 if any(ismember(metIdx,repIdx)) +0093 error('The to-be-replaced metabolite and the replacement metabolite are the same.'); +0094 end +0095 originalStoch = model.S(metIdx,rxnsWithMet); +0096 model.S(repIdx,rxnsWithMet) = model.S(repIdx,rxnsWithMet) + originalStoch; +0097 model.S(metIdx,rxnsWithMet) = 0; +0098 idxDelete = metIdx; +0099 else +0100 % Run through replacement metabolites and their compartments. If any of the +0101 % to-be-replaced metabolites is already present (checked by +0102 % metaboliteName[compartment], then the replacement metabolite is kept and +0103 % the to-be-replace metabolite ID deleted. +0104 +0105 % Build list of metaboliteName[compartment] +0106 comps = model.comps(:); +0107 metCompsN = strcat(lower(model.metNames(:)),'[',comps(model.metComps(:)),']'); +0108 +0109 % Within each metaboliteName[compartment] group, the first replacement +0110 % metabolite is kept and the other members are merged into it. +0111 groups = unique(metCompsN(repIdx),'stable'); +0112 for i = 1:numel(groups) +0113 metCompsNidx = find(strcmp(groups{i}, metCompsN)); +0114 keepIdx = metCompsNidx(find(ismember(metCompsNidx,repIdx),1)); +0115 others = metCompsNidx(metCompsNidx ~= keepIdx); +0116 if ~isempty(others) +0117 model.S(keepIdx,:) = model.S(keepIdx,:) + sum(model.S(others,:),1); +0118 idxDelete = [idxDelete; others]; % Make list of metabolite IDs to delete +0119 end +0120 end +0121 end +0122 if ~isempty(idxDelete) +0123 repMets = model.mets(repIdx); +0124 model = removeMets(model,idxDelete); +0125 % Indices of the replacement metabolites shift when metabolites with a +0126 % lower index are deleted. +0127 [~, repIdx] = ismember(repMets,model.mets); +0128 repIdx(repIdx==0) = []; +0129 end +0130 +0131 % This could now have created duplicate reactions. Contract model. +0132 if isempty(idxDelete) +0133 removedRxns = cell(0,1); +0134 idxDuplRxns = zeros(0,1); +0135 else +0136 originalRxns = model.rxns; +0137 [model, removedRxns] = contractModel(model,[],repIdx); +0138 [~, idxDuplRxns] = ismember(removedRxns,originalRxns); +0139 end +0140 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/testing/unit_tests/modelCurationTests.html b/doc/testing/unit_tests/modelCurationTests.html index dde7fe80c..0846789f0 100644 --- a/doc/testing/unit_tests/modelCurationTests.html +++ b/doc/testing/unit_tests/modelCurationTests.html @@ -42,7 +42,7 @@

    CROSS-REFERENCE INFORMATION ^ <h2><a name=SUBFUNCTIONS ^

    +
  • function addExchangeRxnsTest(testCase)
  • function addTransportTest(testCase)
  • function addGenesRavenTest(testCase)
  • function removeGenesTest(testCase)
  • function addMetsTest(testCase)
  • function addMets_oneCompTest(testCase)
  • function removeMetsTest(testCase)
  • function addRxnsTest(testCase)
  • function removeReactionsTest(testCase)
  • function getMetsInCompTest(testCase)
  • function getRxnsInCompTest(testCase)
  • function getAllRxnsFromGenesTest(testCase)
  • function changeRxnsTest(testCase)
  • function removeBadRxnsTest(testCase)
  • function changeGrRulesTest(testCase)
  • function copyToCompsTest(testCase)
  • function setExchangeBoundsTest(testCase)
  • function addRxnsGenesMetsTest(testCase)
  • function replaceMetsOutputsTest(testCase)
  • function replaceMetsSharedReactionTest(testCase)
  • function replaceMetsNamesTest(testCase)
  • function replaceMetsNamesRepeatedReplacementTest(testCase)

    • SOURCE CODE ^

    0001 %run this test case with the command
@@ -537,7 +537,95 @@ 
    SOURCE CODE ^end
    +0493 end +0494 +0495 function replaceMetsOutputsTest(testCase) +0496 sourceDir = fileparts(which(mfilename)); +0497 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0498 +0499 %Add a duplicate of atp_c and a copy of ATPM that uses it +0500 metsToAdd.mets={'atpDup_c'}; +0501 metsToAdd.metNames={'atp duplicate'}; +0502 metsToAdd.compartments={'c'}; +0503 evalc('model=addMets(model,metsToAdd);'); +0504 rxnsToAdd.rxns={'ATPM_dup'}; +0505 rxnsToAdd.equations={'atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +0506 evalc('model=addRxns(model,rxnsToAdd,1);'); +0507 +0508 [modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'atpDup_c','atp_c',false,true); +0509 +0510 verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +0511 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +0512 verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +0513 verifyEqual(testCase,idxDuplRxns,numel(model.rxns)) +0514 end +0515 +0516 function replaceMetsSharedReactionTest(testCase) +0517 %The replaced metabolite occurs in a reaction that already contains the +0518 %replacement metabolite; the stoichiometries are summed. +0519 sourceDir = fileparts(which(mfilename)); +0520 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0521 +0522 metsToAdd.mets={'atpDup_c'}; +0523 metsToAdd.metNames={'atp duplicate'}; +0524 metsToAdd.compartments={'c'}; +0525 evalc('model=addMets(model,metsToAdd);'); +0526 rxnsToAdd.rxns={'ATPshared'}; +0527 rxnsToAdd.equations={'atp_c + atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +0528 evalc('model=addRxns(model,rxnsToAdd,1);'); +0529 +0530 modelNew=replaceMets(model,'atpDup_c','atp_c',false,true); +0531 +0532 rxnIdx=strcmp(modelNew.rxns,'ATPshared'); +0533 verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),rxnIdx)),-2) +0534 verifyEqual(testCase,numel(modelNew.mets),numel(model.mets)-1) +0535 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +0536 end +0537 +0538 function replaceMetsNamesTest(testCase) +0539 %The replacement metabolite has a higher index than the to-be-replaced +0540 %metabolite, and the removed duplicate reaction is not the last reaction. +0541 sourceDir = fileparts(which(mfilename)); +0542 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0543 +0544 metsToAdd.mets={'atpNew_c'}; +0545 metsToAdd.metNames={'ATPnew'}; +0546 metsToAdd.compartments={'c'}; +0547 evalc('model=addMets(model,metsToAdd);'); +0548 rxnsToAdd.rxns={'ATPM_dup','ATPM_other'}; +0549 rxnsToAdd.equations={'atpNew_c + h2o_c => adp_c + h_c + pi_c','h2o_c => h_c + pi_c'}; +0550 evalc('model=addRxns(model,rxnsToAdd,1);'); +0551 +0552 [modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'ATP','ATPnew'); +0553 +0554 verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +0555 verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +0556 verifyEqual(testCase,idxDuplRxns,numel(model.rxns)-1) +0557 verifyEqual(testCase,model.rxns(idxDuplRxns),removedRxns) +0558 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp_c'))) +0559 end +0560 +0561 function replaceMetsNamesRepeatedReplacementTest(testCase) +0562 %The replacement name occurs twice in the same compartment; the +0563 %stoichiometry is added once. +0564 sourceDir = fileparts(which(mfilename)); +0565 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0566 +0567 metsToAdd.mets={'atp2_c','atpNew_c'}; +0568 metsToAdd.metNames={'ATP2','ATPnew'}; +0569 metsToAdd.compartments={'c','c'}; +0570 evalc('model=addMets(model,metsToAdd);'); +0571 model.metNames{strcmp(model.mets,'atp2_c')}='ATP'; +0572 rxnsToAdd.rxns={'rNew'}; +0573 rxnsToAdd.equations={'atpNew_c => adp_c'}; +0574 evalc('model=addRxns(model,rxnsToAdd,1);'); +0575 +0576 modelNew=replaceMets(model,'ATPnew','ATP'); +0577 +0578 verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),strcmp(modelNew.rxns,'rNew'))),-1) +0579 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpNew_c'))) +0580 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp2_c'))) +0581 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/testing/unit_tests/modelCurationTests.m b/testing/unit_tests/modelCurationTests.m index 62e3dc962..bd6f80900 100755 --- a/testing/unit_tests/modelCurationTests.m +++ b/testing/unit_tests/modelCurationTests.m @@ -491,3 +491,91 @@ function addRxnsGenesMetsTest(testCase) verifyEqual(testCase,modelNew,modelManual) end + +function replaceMetsOutputsTest(testCase) +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +%Add a duplicate of atp_c and a copy of ATPM that uses it +metsToAdd.mets={'atpDup_c'}; +metsToAdd.metNames={'atp duplicate'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPM_dup'}; +rxnsToAdd.equations={'atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +[modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'atpDup_c','atp_c',false,true); + +verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +verifyEqual(testCase,idxDuplRxns,numel(model.rxns)) +end + +function replaceMetsSharedReactionTest(testCase) +%The replaced metabolite occurs in a reaction that already contains the +%replacement metabolite; the stoichiometries are summed. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atpDup_c'}; +metsToAdd.metNames={'atp duplicate'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPshared'}; +rxnsToAdd.equations={'atp_c + atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +modelNew=replaceMets(model,'atpDup_c','atp_c',false,true); + +rxnIdx=strcmp(modelNew.rxns,'ATPshared'); +verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),rxnIdx)),-2) +verifyEqual(testCase,numel(modelNew.mets),numel(model.mets)-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +end + +function replaceMetsNamesTest(testCase) +%The replacement metabolite has a higher index than the to-be-replaced +%metabolite, and the removed duplicate reaction is not the last reaction. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atpNew_c'}; +metsToAdd.metNames={'ATPnew'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPM_dup','ATPM_other'}; +rxnsToAdd.equations={'atpNew_c + h2o_c => adp_c + h_c + pi_c','h2o_c => h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +[modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'ATP','ATPnew'); + +verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +verifyEqual(testCase,idxDuplRxns,numel(model.rxns)-1) +verifyEqual(testCase,model.rxns(idxDuplRxns),removedRxns) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp_c'))) +end + +function replaceMetsNamesRepeatedReplacementTest(testCase) +%The replacement name occurs twice in the same compartment; the +%stoichiometry is added once. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atp2_c','atpNew_c'}; +metsToAdd.metNames={'ATP2','ATPnew'}; +metsToAdd.compartments={'c','c'}; +evalc('model=addMets(model,metsToAdd);'); +model.metNames{strcmp(model.mets,'atp2_c')}='ATP'; +rxnsToAdd.rxns={'rNew'}; +rxnsToAdd.equations={'atpNew_c => adp_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +modelNew=replaceMets(model,'ATPnew','ATP'); + +verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),strcmp(modelNew.rxns,'rNew'))),-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpNew_c'))) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp2_c'))) +end