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     index of removedRxns in original 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     index of removedRxns in original 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 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