addMets

 addMets
   Adds metabolites to a model

 Input:
   model       a model structure
   metsToAdd   the metabolite structure can have the following fields:
       mets            cell array with unique strings that identifies each
                       metabolite (optional, default IDs of new
                       metabolites are numbered with the prefix defined
                       below)
       metNames        cell array with the names of each metabolite
       compartments    cell array with the compartment of each
                       metabolite. Should match model.comps. If this is a
                       string rather than a cell array it is assumed that
                       all mets are in that compartment
       b               Nx1 or Nx2 matrix with equality constraints for
                       each metabolite (optional, default 0)
       unconstrained   vector describing if each metabolite is an exchange
                       metabolite (1) or not (0) (optional, default 0)
       inchis          cell array with InChI strings (optional, default '')
       metSmiles       cell array with SMILES strings (optional, default '')
       metFormulas     cell array with the formulas  (optional, default '')
       metMiriams      cell array with MIRIAM structures (optional, default [])
       metCharges      metabolite charge (optional, default NaN)
       metDeltaG       Gibbs free energy of formation at biochemical
                       standard condition in kJ/mole (optional, default NaN)
       metNotes        cell array with metabolite notes as strings
                       (optional, default '')
   copyInfo    when adding metabolites to a compartment where it
               previously did not exist, the function will copy any
               available annotation from the metabolite in another
               compartment (optional, default true)
   prefix      when metsToAdd.mets is not specified, new metabolite IDs
               are generated with the prefix specified here. If IDs with
               the prefix are already used in the model then the
               numbering will start from the highest existing integer+1
               (optional, default 'm_')

 Output:
   newModel     an updated model structure

 This function does not make extensive checks about MIRIAM formats,
 forbidden characters or such.

 If multiple metabolites are added at once, the metMiriams cell array
 should be defined as (example with ChEBI and KEGG):

 metsToAdd.metMiriams{1} = struct('name',{{'chebi';'kegg.compound'}},...
     'value',{{'CHEBI:18072';'C11821'}});
 metsToAdd.metMiriams{2} = struct('name',{{'chebi';'kegg.compound'}},...
     'value',{{'CHEBI:31132';'C12248'}});

 Usage: newModel = addMets(model, metsToAdd, copyInfo, prefix)

0001 function newModel=addMets(model,metsToAdd,copyInfo,prefix)
0002 % addMets
0003 %   Adds metabolites to a model
0004 %
0005 % Input:
0006 %   model       a model structure
0007 %   metsToAdd   the metabolite structure can have the following fields:
0008 %       mets            cell array with unique strings that identifies each
0009 %                       metabolite (optional, default IDs of new
0010 %                       metabolites are numbered with the prefix defined
0011 %                       below)
0012 %       metNames        cell array with the names of each metabolite
0013 %       compartments    cell array with the compartment of each
0014 %                       metabolite. Should match model.comps. If this is a
0015 %                       string rather than a cell array it is assumed that
0016 %                       all mets are in that compartment
0017 %       b               Nx1 or Nx2 matrix with equality constraints for
0018 %                       each metabolite (optional, default 0)
0019 %       unconstrained   vector describing if each metabolite is an exchange
0020 %                       metabolite (1) or not (0) (optional, default 0)
0021 %       inchis          cell array with InChI strings (optional, default '')
0022 %       metSmiles       cell array with SMILES strings (optional, default '')
0023 %       metFormulas     cell array with the formulas  (optional, default '')
0024 %       metMiriams      cell array with MIRIAM structures (optional, default [])
0025 %       metCharges      metabolite charge (optional, default NaN)
0026 %       metDeltaG       Gibbs free energy of formation at biochemical
0027 %                       standard condition in kJ/mole (optional, default NaN)
0028 %       metNotes        cell array with metabolite notes as strings
0029 %                       (optional, default '')
0030 %   copyInfo    when adding metabolites to a compartment where it
0031 %               previously did not exist, the function will copy any
0032 %               available annotation from the metabolite in another
0033 %               compartment (optional, default true)
0034 %   prefix      when metsToAdd.mets is not specified, new metabolite IDs
0035 %               are generated with the prefix specified here. If IDs with
0036 %               the prefix are already used in the model then the
0037 %               numbering will start from the highest existing integer+1
0038 %               (optional, default 'm_')
0039 %
0040 % Output:
0041 %   newModel     an updated model structure
0042 %
0043 % This function does not make extensive checks about MIRIAM formats,
0044 % forbidden characters or such.
0045 %
0046 % If multiple metabolites are added at once, the metMiriams cell array
0047 % should be defined as (example with ChEBI and KEGG):
0048 %
0049 % metsToAdd.metMiriams{1} = struct('name',{{'chebi';'kegg.compound'}},...
0050 %     'value',{{'CHEBI:18072';'C11821'}});
0051 % metsToAdd.metMiriams{2} = struct('name',{{'chebi';'kegg.compound'}},...
0052 %     'value',{{'CHEBI:31132';'C12248'}});
0053 %
0054 % Usage: newModel = addMets(model, metsToAdd, copyInfo, prefix)
0055 
0056 if nargin<3
0057     copyInfo=true;
0058 end
0059 if nargin<4
0060     prefix='m_';
0061 else
0062     prefix=char(prefix);
0063 end
0064 
0065 newModel=model;
0066 
0067 if isempty(metsToAdd)
0068     return;
0069 end
0070 
0071 %Check some stuff regarding the required fields
0072 if ~isfield(metsToAdd,'mets')
0073     metsToAdd.metNames=convertCharArray(metsToAdd.metNames);
0074     metsToAdd.mets=generateNewIds(newModel,'mets',prefix,numel(metsToAdd.metNames));
0075 else
0076     metsToAdd.mets=convertCharArray(metsToAdd.mets);
0077 end
0078 if ~isfield(metsToAdd,'metNames')
0079     metsToAdd.metNames=metsToAdd.mets;
0080 else
0081     metsToAdd.metNames=convertCharArray(metsToAdd.metNames);
0082 end
0083 if ~isfield(metsToAdd,'compartments')
0084     EM='compartments is a required field in metsToAdd';
0085     dispEM(EM);
0086 else
0087     metsToAdd.compartments=convertCharArray(metsToAdd.compartments);
0088     %If only one compartment is given, assume it is for all metabolites
0089     if numel(metsToAdd.compartments)==1 && numel(metsToAdd.mets)>1
0090         temp=cell(numel(metsToAdd.mets),1);
0091         temp(:)=metsToAdd.compartments;
0092         metsToAdd.compartments=temp;
0093     end
0094 end
0095 
0096 %Number of metabolites
0097 nMets=numel(metsToAdd.mets);
0098 nOldMets=numel(model.mets);
0099 filler=cell(nMets,1);
0100 filler(:)={''};
0101 largeFiller=cell(nOldMets,1);
0102 largeFiller(:)={''};
0103 
0104 %Check that no metabolite ids are already present in the model
0105 I=ismember(metsToAdd.mets,model.mets);
0106 if any(I)
0107     error('One or more elements in metsToAdd.mets are already present in model.mets: %s',...
0108         strjoin(metsToAdd.mets(I),', '));
0109 else
0110     newModel.mets=[newModel.mets;metsToAdd.mets(:)];
0111 end
0112 
0113 %Check that all the compartments could be found
0114 [I, compMap]=ismember(metsToAdd.compartments,model.comps);
0115 if ~all(I)
0116     EM='metsToAdd.compartments must match model.comps';
0117     dispEM(EM);
0118 end
0119 
0120 %Check that the metabolite names aren't present in the same compartment.
0121 %Not the neatest way maybe..
0122 t1=strcat(metsToAdd.metNames(:),'***',metsToAdd.compartments(:));
0123 t2=strcat(model.metNames,'***',model.comps(model.metComps));
0124 if any(ismember(t1,t2))
0125     EM='One or more elements in metsToAdd.metNames already exist in the same compartments as the one it is being added to';
0126     dispEM(EM);
0127 end
0128 
0129 %Some more checks and if they pass then add each field to the structure
0130 if numel(metsToAdd.metNames)~=nMets
0131     EM='metsToAdd.metNames must have the same number of elements as metsToAdd.mets';
0132     dispEM(EM);
0133 else
0134     newModel.metNames=[newModel.metNames;metsToAdd.metNames(:)];
0135 end
0136 
0137 if numel(compMap)~=nMets
0138     EM='metsToAdd.compartments must have the same number of elements as metsToAdd.mets';
0139     dispEM(EM);
0140 else
0141     newModel.metComps=[newModel.metComps;compMap(:)];
0142 end
0143 
0144 if isfield(metsToAdd,'b')
0145     if size(metsToAdd.b,1)~=nMets
0146         EM='metsToAdd.b must have the same number of elements as metsToAdd.mets';
0147         dispEM(EM);
0148     else
0149         %Add empty field if it doesn't exist
0150         if ~isfield(newModel,'b')
0151             newModel.b=zeros(nOldMets,1);
0152         end
0153         
0154         %If the original is only one vector
0155         if size(metsToAdd.b,2)>size(newModel.b,2)
0156             newModel.b=[newModel.b newModel.b];
0157         end
0158         %Add the new ones
0159         newModel.b=[newModel.b;metsToAdd.b];
0160     end
0161 else
0162     if isfield(newModel,'b')
0163         %Add the default
0164         newModel.b=[newModel.b;zeros(nMets,size(newModel.b,2))];
0165     end
0166 end
0167 
0168 if isfield(metsToAdd,'unconstrained')
0169     if numel(metsToAdd.unconstrained)~=nMets
0170         EM='metsToAdd.unconstrained must have the same number of elements as metsToAdd.mets';
0171         dispEM(EM);
0172     else
0173         %Add empty field if it doesn't exist
0174         if ~isfield(newModel,'unconstrained')
0175             newModel.unconstrained=zeros(nOldMets,1);
0176         end
0177         
0178         %Add the new ones
0179         newModel.unconstrained=[newModel.unconstrained;metsToAdd.unconstrained(:)];
0180     end
0181 else
0182     if isfield(newModel,'unconstrained')
0183         %Add the default
0184         newModel.unconstrained=[newModel.unconstrained;zeros(nMets,1)];
0185     end
0186 end
0187 
0188 if isfield(metsToAdd,'inchis')
0189     metsToAdd.inchis=convertCharArray(metsToAdd.inchis);
0190     if numel(metsToAdd.inchis)~=nMets
0191         EM='metsToAdd.inchis must have the same number of elements as metsToAdd.mets';
0192         dispEM(EM);
0193     end
0194     %Add empty field if it doesn't exist
0195     if ~isfield(newModel,'inchis')
0196         newModel.inchis=largeFiller;
0197     end
0198     newModel.inchis=[newModel.inchis;metsToAdd.inchis(:)];
0199 else
0200     %Add empty strings if structure is in model
0201     if isfield(newModel,'inchis')
0202         newModel.inchis=[newModel.inchis;filler];
0203     end
0204 end
0205 
0206 
0207 if isfield(metsToAdd,'metSmiles')
0208     metsToAdd.metSmiles=convertCharArray(metsToAdd.metSmiles);
0209     if numel(metsToAdd.metSmiles)~=nMets
0210         EM='metsToAdd.metSmiles must have the same number of elements as metsToAdd.mets';
0211         dispEM(EM);
0212     end
0213     %Add empty field if it doesn't exist
0214     if ~isfield(newModel,'metSmiles')
0215         newModel.metSmiles=largeFiller;
0216     end
0217     newModel.metSmiles=[newModel.metSmiles;metsToAdd.metSmiles(:)];
0218 else
0219     %Add empty strings if structure is in model
0220     if isfield(newModel,'metSmiles')
0221         newModel.metSmiles=[newModel.metSmiles;filler];
0222     end
0223 end
0224 
0225 if isfield(metsToAdd,'metFormulas')
0226     metsToAdd.metFormulas=convertCharArray(metsToAdd.metFormulas);
0227     if numel(metsToAdd.metFormulas)~=nMets
0228         EM='metsToAdd.metFormulas must have the same number of elements as metsToAdd.mets';
0229         dispEM(EM);
0230     end
0231     %Add empty field if it doesn't exist
0232     if ~isfield(newModel,'metFormulas')
0233         newModel.metFormulas=largeFiller;
0234     end
0235     newModel.metFormulas=[newModel.metFormulas;metsToAdd.metFormulas(:)];
0236 else
0237     %Add default
0238     if isfield(newModel,'metFormulas')
0239         newModel.metFormulas=[newModel.metFormulas;filler];
0240     end
0241 end
0242 
0243 if isfield(metsToAdd,'metCharges')
0244     if numel(metsToAdd.metCharges)~=nMets
0245         EM='metsToAdd.metCharges must have the same number of elements as metsToAdd.mets';
0246         dispEM(EM);
0247     end
0248     if ~isnumeric(metsToAdd.metCharges)
0249         EM='metsToAdd.metCharges must be of type "double"';
0250         dispEM(EM);
0251     end
0252     if ~isfield(newModel,'metCharges')
0253         newModel.metCharges=NaN(numel(largeFiller),1);
0254     end
0255     newModel.metCharges=[newModel.metCharges;metsToAdd.metCharges(:)];
0256 else
0257     %Add default
0258     if isfield(newModel,'metCharges')
0259         newModel.metCharges=[newModel.metCharges;NaN(numel(filler),1)];
0260     end
0261 end
0262 
0263 if isfield(metsToAdd,'metDeltaG')
0264     if numel(metsToAdd.metDeltaG)~=nMets
0265         EM='metsToAdd.metDeltaG must have the same number of elements as metsToAdd.mets';
0266         dispEM(EM);
0267     end
0268     if ~isnumeric(metsToAdd.metDeltaG)
0269         EM='metsToAdd.metDeltaG must be of type "double"';
0270         dispEM(EM);
0271     end
0272     if ~isfield(newModel,'metDeltaG')
0273         newModel.metDeltaG=NaN(numel(largeFiller),1);
0274     end
0275     newModel.metDeltaG=[newModel.metDeltaG;metsToAdd.metDeltaG(:)];
0276 else
0277     %Add default
0278     if isfield(newModel,'metDeltaG')
0279         newModel.metDeltaG=[newModel.metDeltaG;NaN(numel(filler),1)];
0280     end
0281 end
0282 
0283 
0284 if isfield(metsToAdd,'metNotes')
0285     metsToAdd.metNotes=convertCharArray(metsToAdd.metNotes);
0286     if numel(metsToAdd.metNotes)==1 && numel(metsToAdd.mets)>1
0287         temp=cell(numel(metsToAdd.mets),1);
0288         temp(:)=metsToAdd.metNotes;
0289         metsToAdd.metNotes=temp;
0290     end
0291     if numel(metsToAdd.metNotes)~=nMets
0292         EM='metsToAdd.metNotes must have the same number of elements as metsToAdd.mets';
0293         dispEM(EM);
0294     end
0295     %Add empty field if it doesn't exist
0296     if ~isfield(newModel,'metNotes')
0297         newModel.metNotes=largeFiller;
0298     end
0299     newModel.metNotes=[newModel.metNotes;metsToAdd.metNotes(:)];
0300 else
0301     %Add empty strings if structure is in model
0302     if isfield(newModel,'metNotes')
0303         newModel.metNotes=[newModel.metNotes;filler];
0304     end
0305 end
0306 
0307 %Don't check the type of metMiriams
0308 if isfield(metsToAdd,'metMiriams')
0309     if numel(metsToAdd.metMiriams)==1 && numel(metsToAdd.mets)>1
0310         temp=cell(numel(metsToAdd.mets),1);
0311         temp(:)={metsToAdd.metMiriams};
0312         metsToAdd.metMiriams=temp;
0313     end
0314     if numel(metsToAdd.metMiriams)~=nMets
0315         EM='metsToAdd.metMiriams must have the same number of elements as metsToAdd.mets';
0316         dispEM(EM);
0317     end
0318     %Add empty field if it doesn't exist
0319     if ~isfield(newModel,'metMiriams')
0320         newModel.metMiriams=cell(nOldMets,1);
0321     end
0322     newModel.metMiriams=[newModel.metMiriams;metsToAdd.metMiriams(:)];
0323 else
0324     if isfield(newModel,'metMiriams')
0325         newModel.metMiriams=[newModel.metMiriams;cell(nMets,1)];
0326     end
0327 end
0328 
0329 if isfield(newModel,'metFrom')
0330     newModel.metFrom=[newModel.metFrom;filler];
0331 end
0332 
0333 %Expand the S matrix
0334 newModel.S=[newModel.S;sparse(nMets,size(newModel.S,2))];
0335 
0336 if copyInfo==true
0337     [I, J]=ismember(metsToAdd.metNames,model.metNames);
0338     J=J(I);
0339     %I is the indexes of the new metabolites for which a metabolite with
0340     %the same name existed
0341     I=find(I)+nOldMets;
0342     %Go through each of the added mets and copy annotation if it doesn't
0343     %exist
0344     for i=1:numel(I)
0345         if isfield(newModel,'inchis')
0346             if isempty(newModel.inchis{I(i)})
0347                 newModel.inchis(I(i))=newModel.inchis(J(i));
0348             end
0349         end
0350         if isfield(newModel,'metSmiles')
0351             if isempty(newModel.metSmiles{I(i)})
0352                 newModel.metSmiles(I(i))=newModel.metSmiles(J(i));
0353             end
0354         end        
0355         if isfield(newModel,'metFormulas')
0356             if isempty(newModel.metFormulas{I(i)})
0357                 newModel.metFormulas(I(i))=newModel.metFormulas(J(i));
0358             end
0359         end
0360         if isfield(newModel,'metMiriams')
0361             if isempty(newModel.metMiriams{I(i)})
0362                 newModel.metMiriams(I(i))=newModel.metMiriams(J(i));
0363             end
0364         end
0365         if isfield(newModel,'metCharges')
0366             newModel.metCharges(I(i))=newModel.metCharges(J(i));
0367         end
0368         if isfield(newModel,'metDeltaG')
0369             newModel.metDeltaG(I(i))=newModel.metDeltaG(J(i));
0370         end
0371     end
0372 end
0373 end