|
1 | | -function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose) |
| 1 | +function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose,identifiers) |
2 | 2 | % replaceMets |
3 | 3 | % Replaces metabolite names and annotation with replacement metabolite |
4 | 4 | % that is already in the model. If this results in duplicate metabolites, |
|
13 | 13 | % verbose logical whether to print the ids of reactions that |
14 | 14 | % involve the replaced metabolite (optional, default |
15 | 15 | % false) |
| 16 | +% identifiers true if 'metabolite' and 'replacement' refer to |
| 17 | +% metabolite identifiers instead of metabolite names |
| 18 | +% (optional, default false) |
16 | 19 | % |
17 | 20 | % Output: |
18 | 21 | % model model structure with selected metabolites replaced |
19 | 22 | % removedRxns identifiers of duplicate reactions that were removed |
20 | 23 | % idxDuplRxns index of removedRxns in original model |
21 | 24 | % |
22 | 25 | % Note: This function is useful when the model contains both 'oxygen' and |
23 | | -% 'o2' as metabolites. |
| 26 | +% 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, |
| 27 | +% then the 'identifiers' flag should be set to true. |
24 | 28 | % |
25 | 29 | % Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose) |
26 | 30 |
|
27 | 31 | metabolite=char(metabolite); |
28 | 32 | replacement=char(replacement); |
29 | 33 |
|
30 | | -if nargin<4 |
| 34 | +if nargin<4 || isempty(verbose) |
31 | 35 | verbose=false; |
32 | 36 | end |
| 37 | +if nargin<5 |
| 38 | + identifiers = false; |
| 39 | +end |
33 | 40 |
|
34 | 41 | % Find occurence of replacement metabolites. Annotation will be taken from |
35 | | -% first metabolite found. Metabolite ID from replacement will be used where |
36 | | -% possible. |
37 | | -repIdx = find(strcmp(replacement,model.metNames)); |
| 42 | +% first metabolite found. |
| 43 | +if identifiers |
| 44 | + repIdx = find(strcmp(replacement,model.mets)); |
| 45 | +else |
| 46 | + repIdx = find(strcmp(replacement,model.metNames)); |
| 47 | +end |
38 | 48 | if isempty(repIdx) |
39 | | - error('The replacement metabolite name cannot be found in the model.'); |
| 49 | + error('The replacement metabolite cannot be found in the model.'); |
40 | 50 | end |
41 | 51 |
|
42 | 52 | % Change name and information from metabolite to replacement metabolite |
43 | | -metIdx = find(strcmp(metabolite,model.metNames)); |
| 53 | +if identifiers |
| 54 | + metIdx = find(strcmp(metabolite,model.mets)); |
| 55 | +else |
| 56 | + metIdx = find(strcmp(metabolite,model.metNames)); |
| 57 | +end |
44 | 58 | if isempty(metIdx) |
45 | | - error('The to-be-replaced metabolite name cannot be found in the model.'); |
| 59 | + error('The to-be-replaced metabolite cannot be found in the model.'); |
46 | 60 | end |
| 61 | + |
| 62 | +rxnsWithMet = find(model.S(metIdx,:)); |
47 | 63 | if verbose==true |
48 | | - fprintf('\n\nThe following reactions contain the replaced metabolite as reactant:\n') |
49 | | - fprintf(strjoin(model.rxns(find(model.S(metIdx,:))),'\n')) |
| 64 | + fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') |
| 65 | + fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) |
50 | 66 | fprintf('\n') |
51 | 67 | end |
| 68 | + |
52 | 69 | model.metNames(metIdx) = model.metNames(repIdx(1)); |
53 | 70 | if isfield(model,'metFormulas') |
54 | 71 | model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); |
|
68 | 85 | if isfield(model,'metSmiles') |
69 | 86 | model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); |
70 | 87 | end |
71 | | -% Run through replacement metabolites and their compartments. If any of the |
72 | | -% to-be-replaced metabolites is already present (checked by |
73 | | -% metaboliteName[compartment], then the replacement metabolite is kept and |
74 | | -% the to-be-replace metabolite ID deleted. |
75 | | - |
76 | | -% Build list of metaboliteName[compartment] |
77 | | -metCompsN =cellstr(num2str(model.metComps)); |
78 | | -map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); |
79 | | -metCompsN = map.values(metCompsN); |
80 | | -metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); |
81 | 88 |
|
82 | 89 | idxDelete=[]; |
83 | | -for i = 1:length(repIdx) |
84 | | - metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); |
85 | | - if length(metCompsNidx)>1 |
86 | | - for j = 2:length(metCompsNidx) |
87 | | - model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); |
88 | | - idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete |
| 90 | +if identifiers |
| 91 | + originalStoch = model.S(metIdx,rxnsWithMet); |
| 92 | + model.S(repIdx,rxnsWithMet) = originalStoch; |
| 93 | + model.S(metIdx,rxnsWithMet) = 0; |
| 94 | + idxDelete = metIdx; |
| 95 | +else |
| 96 | + % Run through replacement metabolites and their compartments. If any of the |
| 97 | + % to-be-replaced metabolites is already present (checked by |
| 98 | + % metaboliteName[compartment], then the replacement metabolite is kept and |
| 99 | + % the to-be-replace metabolite ID deleted. |
| 100 | + |
| 101 | + % Build list of metaboliteName[compartment] |
| 102 | + metCompsN =cellstr(num2str(model.metComps)); |
| 103 | + map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); |
| 104 | + metCompsN = map.values(metCompsN); |
| 105 | + metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); |
| 106 | + |
| 107 | + for i = 1:length(repIdx) |
| 108 | + metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); |
| 109 | + if length(metCompsNidx)>1 |
| 110 | + for j = 2:length(metCompsNidx) |
| 111 | + model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); |
| 112 | + idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete |
| 113 | + end |
89 | 114 | end |
90 | 115 | end |
91 | 116 | end |
|
0 commit comments