replaceMets

 replaceMets
   Replaces metabolite names and annotation with replacement metabolite
   that is already in the model. If this results in duplicate metabolites,
   the replacement metabolite will be kept, while the S matrix is updated
   to use the replacement metabolite instead.

   model               a model structure
   metabolite          string with name of metabolite to be replace
   replacement         string with name of replacement metabolite
   verbose             logical whether to print the ids of reactions that
                       involve the replaced metabolite (opt, default
                       false)

   This function is useful when the model contains both 'oxygen' and 'o2'
   as metabolites.

   Usage: model=replaceMets(model,metabolite,replacement,verbose)

0001 function model=replaceMets(model,metabolite,replacement,verbose)
0002 % replaceMets
0003 %   Replaces metabolite names and annotation with replacement metabolite
0004 %   that is already in the model. If this results in duplicate metabolites,
0005 %   the replacement metabolite will be kept, while the S matrix is updated
0006 %   to use the replacement metabolite instead.
0007 %
0008 %   model               a model structure
0009 %   metabolite          string with name of metabolite to be replace
0010 %   replacement         string with name of replacement metabolite
0011 %   verbose             logical whether to print the ids of reactions that
0012 %                       involve the replaced metabolite (opt, default
0013 %                       false)
0014 %
0015 %   This function is useful when the model contains both 'oxygen' and 'o2'
0016 %   as metabolites.
0017 %
0018 %   Usage: model=replaceMets(model,metabolite,replacement,verbose)
0019 
0020 metabolite=char(metabolite);
0021 replacement=char(replacement);
0022 
0023 if nargin<4
0024     verbose=false;
0025 end
0026 
0027 % Find occurence of replacement metabolites. Annotation will be taken from
0028 % first metabolite found. Metabolite ID from replacement will be used where
0029 % possible.
0030 repIdx = find(strcmp(replacement,model.metNames));
0031 if isempty(repIdx)
0032     error('The replacement metabolite name cannot be found in the model.');
0033 end
0034 
0035 % Change name and information from metabolite to replacement metabolite
0036 metIdx = find(strcmp(metabolite,model.metNames));
0037 if isempty(metIdx)
0038     error('The to-be-replaced metabolite name cannot be found in the model.');
0039 end
0040 if verbose==true
0041     fprintf('\n\nThe following reactions contain the replaced metabolite as reactant:\n')
0042     fprintf(strjoin(model.rxns(find(model.S(metIdx,:))),'\n'))
0043     fprintf('\n')
0044 end
0045 model.metNames(metIdx) = model.metNames(repIdx(1));
0046 if isfield(model,'metFormulas')
0047     model.metFormulas(metIdx) = model.metFormulas(repIdx(1));
0048 end
0049 if isfield(model,'metMiriams')
0050     model.metMiriams(metIdx) = model.metMiriams(repIdx(1));
0051 end
0052 if isfield(model,'metCharges')
0053     model.metCharges(metIdx) = model.metCharges(repIdx(1));
0054 end
0055 if isfield(model,'metDeltaG')
0056     model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1));
0057 end
0058 if isfield(model,'inchis')
0059     model.inchis(metIdx) = model.inchis(repIdx(1));
0060 end
0061 if isfield(model,'metSmiles')
0062     model.metSmiles(metIdx) = model.metSmiles(repIdx(1));
0063 end
0064 % Run through replacement metabolites and their compartments. If any of the
0065 % to-be-replaced metabolites is already present (checked by
0066 % metaboliteName[compartment], then the replacement metabolite is kept and
0067 % the to-be-replace metabolite ID deleted.
0068 
0069 % Build list of metaboliteName[compartment]
0070 metCompsN =cellstr(num2str(model.metComps));
0071 map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps);
0072 metCompsN = map.values(metCompsN);
0073 metCompsN = strcat(lower(model.metNames),'[',metCompsN,']');
0074 
0075 idxDelete=[];
0076 for i = 1:length(repIdx)
0077     metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN));
0078     if length(metCompsNidx)>1
0079         for j = 2:length(metCompsNidx)
0080             model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:);
0081             idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete
0082         end
0083     end
0084 end
0085 
0086 if ~isempty(idxDelete)
0087     model.S(idxDelete,:) =[];
0088     model.mets(idxDelete) = [];
0089     model.metNames(idxDelete) = [];
0090     model.metComps(idxDelete) = [];
0091     model.b(idxDelete) = [];
0092     if isfield(model,'metFormulas')
0093         model.metFormulas(idxDelete) = [];
0094     end
0095     if isfield(model,'unconstrained')
0096         model.unconstrained(idxDelete) = [];
0097     end
0098     if isfield(model,'metMiriams')
0099         model.metMiriams(idxDelete) = [];
0100     end
0101     if isfield(model,'metCharges')
0102         model.metCharges(idxDelete) = [];
0103     end
0104     if isfield(model,'metDeltaG')
0105         model.metDeltaG(idxDelete) = [];
0106     end
0107     if isfield(model,'inchis')
0108         model.inchis(idxDelete) = [];
0109     end
0110     if isfield(model,'metSmiles')
0111         model.metSmiles(idxDelete) = [];
0112     end
0113     if isfield(model,'metFrom')
0114         model.metFrom(idxDelete) = [];
0115     end
0116 end
0117 
0118 % This could now have created duplicate reactions. Contract model.
0119 model=contractModel(model);
0120 end