addRxnsGenesMets

 addRxnsGenesMets
   Copies reactions from a source model to a new model, including
   (new) metabolites and genes

   model           draft model where reactions should be copied to
   sourceModel     model where reactions and metabolites are sourced from
   rxns            cell array with reaction IDs (from source model). Can also
                   be string if only one reaction is added
   addGene         three options:
                   false   no genes are annotated to the new reactions
                   true    grRules ared copied from the sourceModel and
                           new genes are added when required
                   string or cell array
                           new grRules are specified as string or cell
                           array, and any new genes are added when
                           required
                   (opt, default false)
   rxnNote         cell array with strings explaining why reactions were copied
                   to the model, to be included as newModel.rxnNotes. Can also
                   be string if same rxnNotes should be added for each new
                   reaction, or only one reaction is to be added (opt, default
                   'Added via addRxnsAndMets()')
   confidence      integer specifying confidence score for all reactions.
                   4:  biochemical data: direct evidence from enzymes
                       assays
                   3:  genetic data: knockout/-in or overexpression
                       analysis
                   2:  physiological data: indirect evidence, e.g.
                       secretion products or defined medium requirement
                       sequence data: genome annotation
                   1:  modeling data: required for functional model,
                       hypothetical reaction
                   0:  no evidence
                   following doi:10.1038/nprot.2009.203 (opt, default 0)

   newModel        an updated model structure

     This function only works if the draft model and source model follow
    the same metabolite and compartment naming convention. Metabolites are
    only matched by metaboliteName[compartment]. Useful if one wants to copy
    additional reactions from source to draft after getModelFromHomology was
    used involving the same models.

   Usage: newModel=addRxnsGenesMets(model,sourceModel,rxns,addGene,rxnNote,confidence)

0001 function model=addRxnsGenesMets(model,sourceModel,rxns,addGene,rxnNote,confidence)
0002 % addRxnsGenesMets
0003 %   Copies reactions from a source model to a new model, including
0004 %   (new) metabolites and genes
0005 %
0006 %   model           draft model where reactions should be copied to
0007 %   sourceModel     model where reactions and metabolites are sourced from
0008 %   rxns            cell array with reaction IDs (from source model). Can also
0009 %                   be string if only one reaction is added
0010 %   addGene         three options:
0011 %                   false   no genes are annotated to the new reactions
0012 %                   true    grRules ared copied from the sourceModel and
0013 %                           new genes are added when required
0014 %                   string or cell array
0015 %                           new grRules are specified as string or cell
0016 %                           array, and any new genes are added when
0017 %                           required
0018 %                   (opt, default false)
0019 %   rxnNote         cell array with strings explaining why reactions were copied
0020 %                   to the model, to be included as newModel.rxnNotes. Can also
0021 %                   be string if same rxnNotes should be added for each new
0022 %                   reaction, or only one reaction is to be added (opt, default
0023 %                   'Added via addRxnsAndMets()')
0024 %   confidence      integer specifying confidence score for all reactions.
0025 %                   4:  biochemical data: direct evidence from enzymes
0026 %                       assays
0027 %                   3:  genetic data: knockout/-in or overexpression
0028 %                       analysis
0029 %                   2:  physiological data: indirect evidence, e.g.
0030 %                       secretion products or defined medium requirement
0031 %                       sequence data: genome annotation
0032 %                   1:  modeling data: required for functional model,
0033 %                       hypothetical reaction
0034 %                   0:  no evidence
0035 %                   following doi:10.1038/nprot.2009.203 (opt, default 0)
0036 %
0037 %   newModel        an updated model structure
0038 %
0039 %     This function only works if the draft model and source model follow
0040 %    the same metabolite and compartment naming convention. Metabolites are
0041 %    only matched by metaboliteName[compartment]. Useful if one wants to copy
0042 %    additional reactions from source to draft after getModelFromHomology was
0043 %    used involving the same models.
0044 %
0045 %   Usage: newModel=addRxnsGenesMets(model,sourceModel,rxns,addGene,rxnNote,confidence)
0046 
0047 if nargin<6
0048     confidence=0;
0049 end
0050 rxns=convertCharArray(rxns);
0051 if nargin<5
0052     rxnNote={'Added via addRxnsGenesMets()'};
0053 else
0054     rxnNote=convertCharArray(rxnNote);
0055 end
0056 if numel(rxnNote)==1 && numel(rxns)>1
0057     rxnNoteArray=cell(1,numel(rxns));
0058     rxnNoteArray(:)=rxnNote;
0059     rxnNote=rxnNoteArray;
0060 end
0061 if nargin<4
0062     addGene=false;
0063 end
0064 
0065 % Obtain indexes of reactions in source model
0066 [notNewRxn,oldRxn]=ismember(rxns,model.rxns);
0067 rxns=rxns(~notNewRxn);
0068 if isempty(rxns)
0069     error('All reactions are already in the model.');
0070 elseif ~isempty(notNewRxn)
0071     fprintf('The following reactions were already present in the model and will not be added:\n')
0072     fprintf([strjoin(model.rxns(oldRxn(find(oldRxn))),'\n') '\n'])
0073 end
0074 
0075 [~, rxnIdx]=ismember(rxns,sourceModel.rxns); % Get rxnIDs
0076 if any(rxnIdx==0)
0077     dispEM('The following reaction IDs could not be found in the source model:',true,...
0078         rxns(rxnIdx==0));
0079 end
0080 
0081 % Add new metabolites
0082 metIdx=find(any(sourceModel.S(:,rxnIdx),2)); % Get metabolite IDs
0083 % Many of the metabolites in are already in the draft model, so only add
0084 % the new metabolites
0085 
0086 % Match by metNames[metComps]. First make these structures for each model.
0087 metCompsN =cellstr(num2str(model.metComps));
0088 map=containers.Map(cellstr(num2str(transpose([1:length(model.comps)]))),model.comps);
0089 metCompsN = map.values(metCompsN);
0090 metCompsN = strcat(model.metNames,'[',metCompsN,']');
0091 
0092 sourcemetCompsN=cellstr(num2str(sourceModel.metComps));
0093 map=containers.Map(cellstr(num2str(transpose([1:length(sourceModel.comps)]))),sourceModel.comps);
0094 sourcemetCompsN = map.values(sourcemetCompsN);
0095 sourcemetCompsN=strcat(sourceModel.metNames,'[',sourcemetCompsN,']');
0096 
0097 newMetCompsN=sourcemetCompsN(metIdx);
0098 notNewMet=newMetCompsN(ismember(newMetCompsN,metCompsN));
0099 
0100 if ~isempty(notNewMet)
0101     fprintf('\nThe following metabolites were already present in the model and will not be added:\n')
0102     fprintf([strjoin(transpose(notNewMet),'\n') '\n'])
0103 end
0104 
0105 metIdx=metIdx(~ismember(sourcemetCompsN(metIdx),metCompsN));
0106 
0107 if ~isempty(metIdx)
0108     fprintf('\nThe following metabolites will be added to the model:\n')
0109     fprintf([strjoin(transpose(sourcemetCompsN(metIdx)),'\n') '\n'])
0110     
0111     if isfield(sourceModel,'mets')
0112         metsToAdd.mets=sourceModel.mets(metIdx);
0113     end
0114     if isfield(sourceModel,'metNames')
0115         metsToAdd.metNames=sourceModel.metNames(metIdx);
0116     end
0117     if isfield(sourceModel,'metFormulas')
0118         metsToAdd.metFormulas=sourceModel.metFormulas(metIdx);
0119     end
0120     if isfield(sourceModel,'metCharges')
0121         metsToAdd.metCharges=sourceModel.metCharges(metIdx);
0122     end
0123     if isfield(sourceModel,'metMiriams')
0124         metsToAdd.metMiriams=sourceModel.metMiriams(metIdx);
0125     end
0126     if isfield(sourceModel,'metNotes')
0127         metsToAdd.metNotes=sourceModel.metNotes(metIdx);
0128     end
0129     if isfield(sourceModel,'inchis')
0130         metsToAdd.inchis=sourceModel.inchis(metIdx);
0131     end
0132     if isfield(sourceModel,'metSmiles')
0133         metsToAdd.metSmiles=sourceModel.metSmiles(metIdx);
0134     end
0135     if isfield(sourceModel,'metDeltaG')
0136         metsToAdd.metDeltaG=sourceModel.metDeltaG(metIdx);
0137     end
0138     
0139     metsToAdd.compartments=strtrim(cellstr(num2str(sourceModel.metComps(metIdx)))); % Convert from compartment string to compartment number
0140     [~,idx]=ismember(metsToAdd.compartments,strsplit(num2str(1:length(sourceModel.comps)))); % Match compartment number to compartment abbreviation
0141     metsToAdd.compartments=sourceModel.comps(idx); % Fill in compartment abbreviations
0142     
0143     model=addMets(model,metsToAdd);
0144 end
0145 fprintf('\nNumber of metabolites added to the model:\n')
0146 fprintf([num2str(numel(metIdx)),'\n'])
0147 
0148 % Add new genes
0149 if ~islogical(addGene) | addGene ~= false
0150     if ischar(addGene)
0151         rxnToAdd.grRules={addGene};
0152     elseif iscell(addGene)
0153         rxnToAdd.grRules=addGene(~notNewRxn); 
0154     else
0155         rxnToAdd.grRules=sourceModel.grRules(rxnIdx); % Get the relevant grRules
0156     end
0157 end
0158 % Add new reactions
0159 rxnToAdd.equations=constructEquations(sourceModel,rxnIdx);
0160 rxnToAdd.rxnNames=sourceModel.rxnNames(rxnIdx);
0161 rxnToAdd.rxns=sourceModel.rxns(rxnIdx);
0162 rxnToAdd.lb=sourceModel.lb(rxnIdx);
0163 rxnToAdd.ub=sourceModel.ub(rxnIdx);
0164 rxnToAdd.rxnNotes(:)=rxnNote(~notNewRxn);
0165 rxnToAdd.rxnConfidenceScores=NaN(1,numel(rxnToAdd.rxns));
0166 if ~isnumeric(confidence)
0167     EM='confidence score must be numeric';
0168     dispEM(EM, true);
0169 end
0170 rxnToAdd.rxnConfidenceScores(:)=confidence;
0171 if isfield(sourceModel,'rxnDeltaG')
0172     rxnToAdd.rxnDeltaG=sourceModel.rxnDeltaG(rxnIdx);
0173 end
0174 if isfield(sourceModel,'subSystems')
0175     rxnToAdd.subSystems=sourceModel.subSystems(rxnIdx);
0176 end
0177 if isfield(sourceModel,'eccodes')
0178     rxnToAdd.eccodes=sourceModel.eccodes(rxnIdx);
0179 end
0180 model=addRxns(model,rxnToAdd,3,'',false,true);
0181 
0182 fprintf('\nNumber of reactions added to the model:\n')
0183 fprintf([num2str(numel(rxnIdx)),'\n'])
0184 end