checkModelStruct

 checkModelStruct
   Performs a number of checks to ensure that a model structure is ok

   model           a model structure
   throwErrors     true if the function should throw errors if
                   inconsistencies are found. The alternative is to
                   print warnings for all types of issues (opt, default true)
   trimWarnings    true if only a maximal of 10 items should be displayed in
                   a given error/warning (opt, default true)

   NOTE: This is performed after importing a model from Excel or before
   attempting to export a model to SBML format.

   Usage: checkModelStruct(model,throwErrors,trimWarnings)

0001 function checkModelStruct(model,throwErrors,trimWarnings)
0002 % checkModelStruct
0003 %   Performs a number of checks to ensure that a model structure is ok
0004 %
0005 %   model           a model structure
0006 %   throwErrors     true if the function should throw errors if
0007 %                   inconsistencies are found. The alternative is to
0008 %                   print warnings for all types of issues (opt, default true)
0009 %   trimWarnings    true if only a maximal of 10 items should be displayed in
0010 %                   a given error/warning (opt, default true)
0011 %
0012 %   NOTE: This is performed after importing a model from Excel or before
0013 %   attempting to export a model to SBML format.
0014 %
0015 %   Usage: checkModelStruct(model,throwErrors,trimWarnings)
0016 
0017 if nargin<2
0018     throwErrors=true;
0019 end
0020 if nargin<3
0021     trimWarnings=true;
0022 end
0023 
0024 %Missing elements
0025 fields={'id';'name';'rxns';'mets';'S';'lb';'ub';'rev';'c';'b';'comps';'metComps'};
0026 for i=1:numel(fields)
0027     if ~isfield(model,fields{i})
0028         EM=['The model is missing the "' fields{i} '" field'];
0029         dispEM(EM,throwErrors);
0030     end
0031 end
0032 
0033 %Type check
0034 if ~ischar(model.id)
0035     EM='The "id" field must be a string';
0036     dispEM(EM,throwErrors);
0037 end
0038 if ~ischar(model.name)
0039     EM='The "name" field must be a string';
0040     dispEM(EM,throwErrors);
0041 end
0042 if ~iscellstr(model.rxns)
0043     EM='The "rxns" field must be a cell array of strings';
0044     dispEM(EM,throwErrors);
0045 end
0046 if ~iscellstr(model.mets)
0047     EM='The "mets" field must be a cell array of strings';
0048     dispEM(EM,throwErrors);
0049 end
0050 if ~isnumeric(model.S)
0051     EM='The "S" field must be of type "double"';
0052     dispEM(EM,throwErrors);
0053 end
0054 if ~isnumeric(model.lb)
0055     EM='The "lb" field must be of type "double"';
0056     dispEM(EM,throwErrors);
0057 end
0058 if ~isnumeric(model.ub)
0059     EM='The "ub" field must be of type "double"';
0060     dispEM(EM,throwErrors);
0061 end
0062 if ~isnumeric(model.rev)
0063     EM='The "rev" field must be of type "double"';
0064     dispEM(EM,throwErrors);
0065 end
0066 if ~isnumeric(model.c)
0067     EM='The "c" field must be of type "double"';
0068     dispEM(EM,throwErrors);
0069 end
0070 if ~isnumeric(model.b)
0071     EM='The "b" field must be of type "double"';
0072     dispEM(EM,throwErrors);
0073 end
0074 if ~iscellstr(model.comps)
0075     EM='The "comps" field must be a cell array of strings';
0076     dispEM(EM,throwErrors);
0077 end
0078 if ~isnumeric(model.metComps)
0079     EM='The "metComps" field must be of type "double"';
0080     dispEM(EM,throwErrors);
0081 end
0082 if isfield(model,'compNames')
0083     if ~iscellstr(model.compNames)
0084         EM='The "compNames" field must be a cell array of strings';
0085         dispEM(EM,throwErrors);
0086     end
0087 end
0088 if isfield(model,'compOutside')
0089     if ~iscellstr(model.compOutside)
0090         EM='The "compOutside" field must be a cell array of strings';
0091         dispEM(EM,throwErrors);
0092     end
0093 end
0094 if isfield(model,'rxnNames')
0095     if ~iscellstr(model.rxnNames)
0096         EM='The "rxnNames" field must be a cell array of strings';
0097         dispEM(EM,throwErrors);
0098     end
0099 end
0100 if isfield(model,'metNames')
0101     if ~iscellstr(model.metNames)
0102         EM='The "metNames" field must be a cell array of strings';
0103         dispEM(EM,throwErrors);
0104     end
0105 end
0106 if isfield(model,'genes')
0107     if ~iscellstr(model.genes)
0108         EM='The "genes" field must be a cell array of strings';
0109         dispEM(EM,throwErrors);
0110     end
0111 end
0112 if isfield(model,'rxnGeneMat')
0113     if ~isnumeric(model.rxnGeneMat)
0114         EM='The "rxnGeneMat" field must be of type "double"';
0115         dispEM(EM,throwErrors);
0116     end
0117 end
0118 if isfield(model,'grRules')
0119     if ~iscellstr(model.grRules)
0120         EM='The "grRules" field must be a cell array of strings';
0121         dispEM(EM,throwErrors);
0122     end
0123 end
0124 if isfield(model,'rxnComps')
0125     if ~isnumeric(model.rxnComps)
0126         EM='The "rxnComps" field must be of type "double"';
0127         dispEM(EM,throwErrors);
0128     end
0129 end
0130 if isfield(model,'inchis')
0131     if ~iscellstr(model.inchis)
0132         EM='The "inchis" field must be a cell array of strings';
0133         dispEM(EM,throwErrors);
0134     end
0135 end
0136 if isfield(model,'metSmiles')
0137     if ~iscellstr(model.metSmiles)
0138         EM='The "metSmiles" field must be a cell array of strings';
0139         dispEM(EM,throwErrors);
0140     end
0141 end
0142 if isfield(model,'metFormulas')
0143     if ~iscellstr(model.metFormulas)
0144         EM='The "metFormulas" field must be a cell array of strings';
0145         dispEM(EM,throwErrors);
0146     end
0147 end
0148 if isfield(model,'metCharges')
0149     if ~isnumeric(model.metCharges)
0150         EM='The "metCharges" field must be a double';
0151         dispEM(EM,throwErrors);
0152     end
0153 end
0154 if isfield(model,'metDeltaG')
0155     if ~isnumeric(model.metDeltaG)
0156         EM='The "metDeltaG" field must be a double';
0157         dispEM(EM,throwErrors);
0158     end
0159 end
0160 if isfield(model,'subSystems')
0161     for i=1:numel(model.subSystems)
0162         if ~iscell(model.subSystems{i,1})
0163             EM='The "subSystems" field must be a cell array';
0164             dispEM(EM,throwErrors);
0165         end
0166     end
0167 end
0168 if isfield(model,'eccodes')
0169     if ~iscellstr(model.eccodes)
0170         EM='The "eccodes" field must be a cell array of strings';
0171         dispEM(EM,throwErrors);
0172     end
0173 end
0174 if isfield(model,'unconstrained')
0175     if ~isnumeric(model.unconstrained)
0176         EM='The "unconstrained" field must be of type "double"';
0177         dispEM(EM,throwErrors);
0178     end
0179 end
0180 if isfield(model,'rxnNotes')
0181     if ~iscellstr(model.rxnNotes)
0182         EM='The "rxnNotes" field must be a cell array of strings';
0183         dispEM(EM,throwErrors);
0184     end
0185 end
0186 if isfield(model,'rxnReferences')
0187     if ~iscellstr(model.rxnReferences)
0188         EM='The "rxnReferences" field must be a cell array of strings';
0189         dispEM(EM,throwErrors);
0190     end
0191 end
0192 if isfield(model,'rxnConfidenceScores')
0193     if ~isnumeric(model.rxnConfidenceScores)
0194         EM='The "rxnConfidenceScores" field must be a double';
0195         dispEM(EM,throwErrors);
0196     end
0197 end
0198 if isfield(model,'rxnDeltaG')
0199     if ~isnumeric(model.rxnDeltaG)
0200         EM='The "rxnDeltaG" field must be a double';
0201         dispEM(EM,throwErrors);
0202     end
0203 end
0204 
0205 %Empty strings
0206 if isempty(model.id)
0207     EM='The "id" field cannot be empty';
0208     dispEM(EM,throwErrors);
0209 end
0210 if any(cellfun(@isempty,model.rxns))
0211     EM='The model contains empty reaction IDs';
0212     dispEM(EM,throwErrors);
0213 end
0214 if any(cellfun(@isempty,model.mets))
0215     EM='The model contains empty metabolite IDs';
0216     dispEM(EM,throwErrors);
0217 end
0218 if any(cellfun(@isempty,model.comps))
0219     EM='The model contains empty compartment IDs';
0220     dispEM(EM,throwErrors);
0221 end
0222 EM='The following metabolites have empty names:';
0223 dispEM(EM,throwErrors,model.mets(cellfun(@isempty,model.metNames)),trimWarnings);
0224 
0225 if isfield(model,'genes')
0226     if any(cellfun(@isempty,model.genes))
0227         EM='The model contains empty gene IDs';
0228         dispEM(EM,throwErrors);
0229     end
0230 end
0231 
0232 %Duplicates
0233 EM='The following reaction IDs are duplicates:';
0234 dispEM(EM,throwErrors,model.rxns(duplicates(model.rxns)),trimWarnings);
0235 EM='The following metabolite IDs are duplicates:';
0236 dispEM(EM,throwErrors,model.mets(duplicates(model.mets)),trimWarnings);
0237 EM='The following compartment IDs are duplicates:';
0238 dispEM(EM,throwErrors,model.comps(duplicates(model.comps)),trimWarnings);
0239 if isfield(model,'genes')
0240     EM='The following genes are duplicates:';
0241     dispEM(EM,throwErrors,model.genes(duplicates(model.genes)),trimWarnings);
0242 end
0243 metInComp=strcat(model.metNames,'[',model.comps(model.metComps),']');
0244 EM='The following metabolites already exist in the same compartment:';
0245 dispEM(EM,throwErrors,metInComp(duplicates(metInComp)),trimWarnings);
0246 
0247 %Elements never used (print only as warnings
0248 EM='The following reactions are empty (no involved metabolites):';
0249 dispEM(EM,false,model.rxns(~any(model.S,1)),trimWarnings);
0250 EM='The following metabolites are never used in a reaction:';
0251 dispEM(EM,false,model.mets(~any(model.S,2)),trimWarnings);
0252 if isfield(model,'genes')
0253     EM='The following genes are not associated to a reaction:';
0254     dispEM(EM,false,model.genes(~any(model.rxnGeneMat,1)),trimWarnings);
0255 end
0256 I=true(numel(model.comps),1);
0257 I(model.metComps)=false;
0258 EM='The following compartments contain no metabolites:';
0259 dispEM(EM,false,model.comps(I),trimWarnings);
0260 
0261 %Contradicting bounds
0262 EM='The following reactions have contradicting bounds:';
0263 dispEM(EM,throwErrors,model.rxns(model.lb>model.ub),trimWarnings);
0264 EM='The following reactions have bounds contradicting their reversibility:';
0265 dispEM(EM,throwErrors,model.rxns(model.lb<0 & model.rev==0),trimWarnings);
0266 
0267 %Multiple or no objective functions not allowed in SBML L3V1 FBCv2
0268 if numel(find(model.c))>1
0269     EM='Multiple objective functions found. This might be intended, but results in FBCv2 non-compliant SBML file when exported';
0270     dispEM(EM,false,model.rxns(find(model.c)),trimWarnings);
0271 elseif ~any(model.c)
0272     EM='No objective function found. This might be intended, but results in FBCv2 non-compliant SBML file when exported';
0273     dispEM(EM,false);
0274 end
0275     
0276 EM='The following reactions have contradicting bounds:';
0277 dispEM(EM,throwErrors,model.rxns(model.lb>model.ub),trimWarnings);
0278 
0279 %Mapping of compartments
0280 if isfield(model,'compOutside')
0281     EM='The following compartments are in "compOutside" but not in "comps":';
0282     dispEM(EM,throwErrors,setdiff(model.compOutside,[{''};model.comps]),trimWarnings);
0283 end
0284 
0285 %Met names which start with number
0286 I=false(numel(model.metNames),1);
0287 for i=1:numel(model.metNames)
0288     index=strfind(model.metNames{i},' ');
0289     if any(index)
0290         if any(str2double(model.metNames{i}(1:index(1)-1)))
0291             I(i)=true;
0292         end
0293     end
0294 end
0295 EM='The following metabolite IDs begin with a number directly followed by space:';
0296 dispEM(EM,throwErrors,model.mets(I),trimWarnings);
0297 
0298 %Non-parseable composition
0299 if isfield(model,'metFormulas')
0300     [~, ~, exitFlag]=parseFormulas(model.metFormulas,true,false);
0301     EM='The composition for the following metabolites could not be parsed:';
0302     dispEM(EM,false,model.mets(exitFlag==-1),trimWarnings);
0303 end
0304 
0305 %Check if there are metabolites with different names but the same MIRIAM
0306 %codes
0307 if isfield(model,'metMiriams')
0308     miriams=containers.Map();
0309     for i=1:numel(model.mets)
0310         if ~isempty(model.metMiriams{i})
0311             %Loop through and add for each miriam
0312             for j=1:numel(model.metMiriams{i}.name)
0313                 %Get existing metabolite indexes
0314                 current=strcat(model.metMiriams{i}.name{j},'/',model.metMiriams{i}.value{j});
0315                 if isKey(miriams,current)
0316                     existing=miriams(current);
0317                 else
0318                     existing=[];
0319                 end
0320                 miriams(current)=[existing;i];
0321             end
0322         end
0323     end
0324     
0325     %Get all keys
0326     allMiriams=keys(miriams);
0327     
0328     hasMultiple=false(numel(allMiriams),1);
0329     for i=1:numel(allMiriams)
0330         if numel(miriams(allMiriams{i}))>1
0331             %Check if they all have the same name
0332             if numel(unique(model.metNames(miriams(allMiriams{i}))))>1
0333                 if ~regexp(allMiriams{i},'^sbo\/SBO:') % SBO terms are expected to be multiple
0334                     hasMultiple(i)=true;
0335                 end                
0336             end
0337         end
0338     end
0339     
0340     %Print output
0341     EM='The following MIRIAM strings are associated to more than one unique metabolite name:';
0342     dispEM(EM,false,allMiriams(hasMultiple),trimWarnings);
0343 end
0344 
0345 %Check if there are metabolites with different names but the same InChI
0346 %codes
0347 if isfield(model,'inchis')
0348     inchis=containers.Map();
0349     for i=1:numel(model.mets)
0350         if ~isempty(model.inchis{i})
0351             %Get existing metabolite indexes
0352             if isKey(inchis,model.inchis{i})
0353                 existing=inchis(model.inchis{i});
0354             else
0355                 existing=[];
0356             end
0357             inchis(model.inchis{i})=[existing;i];
0358         end
0359     end
0360     
0361     %Get all keys
0362     allInchis=keys(inchis);
0363     
0364     hasMultiple=false(numel(allInchis),1);
0365     for i=1:numel(allInchis)
0366         if numel(inchis(allInchis{i}))>1
0367             %Check if they all have the same name
0368             if numel(unique(model.metNames(inchis(allInchis{i}))))>1
0369                 hasMultiple(i)=true;
0370             end
0371         end
0372     end
0373     
0374     %Print output
0375     EM='The following InChI strings are associated to more than one unique metabolite name:';
0376     dispEM(EM,false,allInchis(hasMultiple),trimWarnings);
0377 end
0378 
0379 % %Check if there are metabolites with different names but the same SMILES
0380 % if isfield(model,'metSmiles')
0381 %     metSmiles=containers.Map();
0382 %     for i=1:numel(model.mets)
0383 %         if ~isempty(model.metSmiles{i})
0384 %             %Get existing metabolite indexes
0385 %             if isKey(metSmiles,model.metSmiles{i})
0386 %                 existing=metSmiles(model.metSmiles{i});
0387 %             else
0388 %                 existing=[];
0389 %             end
0390 %             metSmiles(model.metSmiles{i})=[existing;i];
0391 %         end
0392 %     end
0393 %
0394 %     %Get all keys
0395 %     allmetSmiles=keys(metSmiles);
0396 %
0397 %     hasMultiple=false(numel(metSmiles),1);
0398 %     for i=1:numel(metSmiles)
0399 %         if numel(metSmiles(metSmiles{i}))>1
0400 %             %Check if they all have the same name
0401 %             if numel(unique(model.metNames(metSmiles(allmetSmiles{i}))))>1
0402 %                 hasMultiple(i)=true;
0403 %             end
0404 %         end
0405 %     end
0406 %
0407 %     %Print output
0408 %     EM='The following metSmiles strings are associated to more than one unique metabolite name:';
0409 %     dispEM(EM,false,allmetSmiles(hasMultiple),trimWarnings);
0410 % end
0411 end
0412 
0413 function I=duplicates(strings)
0414 I=false(numel(strings),1);
0415 [J, K]=unique(strings);
0416 if numel(J)~=numel(strings)
0417     L=1:numel(strings);
0418     L(K)=[];
0419     I(L)=true;
0420 end
0421 end