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. At the end, contractModel is
   run to remove any duplicate reactions that might have occured.

 Input:
   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 (optional, default
                   false)
 
 Output:
   model           model structure with selected metabolites replaced
   removedRxns     identifiers of duplicate reactions that were removed
   idxDuplRxns     index of removedRxns in original model

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

 Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose)

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