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)
   identifiers     true if 'metabolite' and 'replacement' refer to
                   metabolite identifiers instead of metabolite names
                   (optional, default false)
 
 Output:
   model           model structure with selected metabolites replaced
   removedRxns     identifiers of duplicate reactions that were removed
   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,
 then the 'identifiers' flag should be set to true.

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

0001 function [model, removedRxns, idxDuplRxns]=replaceMets(model,metabolite,replacement,verbose,identifiers)
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 %   identifiers     true if 'metabolite' and 'replacement' refer to
0017 %                   metabolite identifiers instead of metabolite names
0018 %                   (optional, default false)
0019 %
0020 % Output:
0021 %   model           model structure with selected metabolites replaced
0022 %   removedRxns     identifiers of duplicate reactions that were removed
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,
0027 % then the 'identifiers' flag should be set to true.
0028 %
0029 % Usage: [model, removedRxns, idxDuplRxns] = replaceMets(model, metabolite, replacement, verbose)
0030 
0031 metabolite=char(metabolite);
0032 replacement=char(replacement);
0033 
0034 if nargin<4 || isempty(verbose)
0035     verbose=false;
0036 end
0037 if nargin<5
0038     identifiers = false;
0039 end
0040 
0041 % Find occurence of replacement metabolites. Annotation will be taken from
0042 % first metabolite found.
0043 if identifiers
0044     repIdx = find(strcmp(replacement,model.mets));
0045 else
0046     repIdx = find(strcmp(replacement,model.metNames));
0047 end
0048 if isempty(repIdx)
0049     error('The replacement metabolite cannot be found in the model.');
0050 end
0051 
0052 % Change name and information from metabolite to replacement metabolite
0053 if identifiers
0054     metIdx = find(strcmp(metabolite,model.mets));
0055 else
0056     metIdx = find(strcmp(metabolite,model.metNames));
0057 end
0058 if isempty(metIdx)
0059     error('The to-be-replaced metabolite cannot be found in the model.');
0060 end
0061 
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