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 (optional, default true)
   trimWarnings    true if only a maximal of 10 items should be displayed in
                   a given error/warning (optional, 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 (optional, default true)
0009 %   trimWarnings    true if only a maximal of 10 items should be displayed in
0010 %                   a given error/warning (optional, 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     if ~isfield(model,'genes')
0124         EM='If "grRules" field exists, the model should also contain a "genes" field';
0125         dispEM(EM,throwErrors);
0126     else
0127         geneList = getGenesFromGrRules(model.grRules);
0128         geneList = setdiff(unique(geneList),model.genes);
0129         if ~isempty(geneList)
0130             problemGrRules = model.rxns(contains(model.grRules,geneList));
0131             problemGrRules = strjoin(problemGrRules(:),'; ');
0132             EM=['The reaction(s) "' problemGrRules '" contain the following genes in its "grRules" field, but these are not in the "genes" field:'];
0133             dispEM(EM,throwErrors,geneList);
0134         end
0135     end
0136 end
0137 if isfield(model,'rxnComps')
0138     if ~isnumeric(model.rxnComps)
0139         EM='The "rxnComps" field must be of type "double"';
0140         dispEM(EM,throwErrors);
0141     end
0142 end
0143 if isfield(model,'inchis')
0144     if ~iscellstr(model.inchis)
0145         EM='The "inchis" field must be a cell array of strings';
0146         dispEM(EM,throwErrors);
0147     end
0148 end
0149 if isfield(model,'metSmiles')
0150     if ~iscellstr(model.metSmiles)
0151         EM='The "metSmiles" field must be a cell array of strings';
0152         dispEM(EM,throwErrors);
0153     end
0154 end
0155 if isfield(model,'metFormulas')
0156     if ~iscellstr(model.metFormulas)
0157         EM='The "metFormulas" field must be a cell array of strings';
0158         dispEM(EM,throwErrors);
0159     end
0160 end
0161 if isfield(model,'metCharges')
0162     if ~isnumeric(model.metCharges)
0163         EM='The "metCharges" field must be a double';
0164         dispEM(EM,throwErrors);
0165     end
0166 end
0167 if isfield(model,'metDeltaG')
0168     if ~isnumeric(model.metDeltaG)
0169         EM='The "metDeltaG" field must be a double';
0170         dispEM(EM,throwErrors);
0171     end
0172 end
0173 if isfield(model,'subSystems')
0174     for i=1:numel(model.subSystems)
0175         if ~iscell(model.subSystems{i,1})
0176             EM='The "subSystems" field must be a cell array';
0177             dispEM(EM,throwErrors);
0178         end
0179     end
0180 end
0181 if isfield(model,'eccodes')
0182     if ~iscellstr(model.eccodes)
0183         EM='The "eccodes" field must be a cell array of strings';
0184         dispEM(EM,throwErrors);
0185     end
0186 end
0187 if isfield(model,'unconstrained')
0188     if ~isnumeric(model.unconstrained)
0189         EM='The "unconstrained" field must be of type "double"';
0190         dispEM(EM,throwErrors);
0191     end
0192 end
0193 if isfield(model,'rxnNotes')
0194     if ~iscellstr(model.rxnNotes)
0195         EM='The "rxnNotes" field must be a cell array of strings';
0196         dispEM(EM,throwErrors);
0197     end
0198 end
0199 if isfield(model,'rxnReferences')
0200     if ~iscellstr(model.rxnReferences)
0201         EM='The "rxnReferences" field must be a cell array of strings';
0202         dispEM(EM,throwErrors);
0203     end
0204 end
0205 if isfield(model,'rxnConfidenceScores')
0206     if ~isnumeric(model.rxnConfidenceScores)
0207         EM='The "rxnConfidenceScores" field must be a double';
0208         dispEM(EM,throwErrors);
0209     end
0210 end
0211 if isfield(model,'rxnDeltaG')
0212     if ~isnumeric(model.rxnDeltaG)
0213         EM='The "rxnDeltaG" field must be a double';
0214         dispEM(EM,throwErrors);
0215     end
0216 end
0217 
0218 %Empty strings
0219 if isempty(model.id)
0220     EM='The "id" field cannot be empty';
0221     dispEM(EM,throwErrors);
0222 end
0223 if any(cellfun(@isempty,model.rxns))
0224     EM='The model contains empty reaction IDs';
0225     dispEM(EM,throwErrors);
0226 end
0227 if any(cellfun(@isempty,model.mets))
0228     EM='The model contains empty metabolite IDs';
0229     dispEM(EM,throwErrors);
0230 end
0231 if any(cellfun(@isempty,model.comps))
0232     EM='The model contains empty compartment IDs';
0233     dispEM(EM,throwErrors);
0234 end
0235 EM='The following metabolites have empty names:';
0236 dispEM(EM,throwErrors,model.mets(cellfun(@isempty,model.metNames)),trimWarnings);
0237 
0238 if isfield(model,'genes')
0239     if any(cellfun(@isempty,model.genes))
0240         EM='The model contains empty gene IDs';
0241         dispEM(EM,throwErrors);
0242     end
0243 end
0244 
0245 %Validate format of ids
0246 fields      = {'rxns';'mets';'comps';'genes'};
0247 fieldNames  = {'reaction';'metabolite';'compartment';'gene'};
0248 fieldPrefix = {'R_';'M_';'C_';'G_'};
0249 for i=1:numel(fields)
0250     try
0251         numIDs = ~startsWith(model.(fields{i}),regexpPattern('^[a-zA-Z_]'));
0252     catch
0253         numIDs = [];
0254     end
0255     if any(numIDs)
0256         EM = ['The following ' fieldNames{i} ' identifiers do not start '...
0257             'with a letter or _ (conflicting with SBML specifications). '...
0258             'This does not impact RAVEN functionality, but be aware that '...
0259             'exportModel will automatically add ' fieldPrefix{i} ...
0260             ' prefixes to all ' fieldNames{i} ' identifiers:'];
0261         dispEM(EM,false,{model.(fields{i}){numIDs}},trimWarnings);
0262     end
0263 end
0264 
0265 %Duplicates
0266 EM='The following reaction IDs are duplicates:';
0267 dispEM(EM,throwErrors,model.rxns(duplicates(model.rxns)),trimWarnings);
0268 EM='The following metabolite IDs are duplicates:';
0269 dispEM(EM,throwErrors,model.mets(duplicates(model.mets)),trimWarnings);
0270 EM='The following compartment IDs are duplicates:';
0271 dispEM(EM,throwErrors,model.comps(duplicates(model.comps)),trimWarnings);
0272 if isfield(model,'genes')
0273     EM='The following genes are duplicates:';
0274     dispEM(EM,throwErrors,model.genes(duplicates(model.genes)),trimWarnings);
0275 end
0276 metInComp=strcat(model.metNames,'[',model.comps(model.metComps),']');
0277 EM='The following metabolites already exist in the same compartment:';
0278 dispEM(EM,throwErrors,metInComp(duplicates(metInComp)),trimWarnings);
0279 
0280 %Elements never used (print only as warnings
0281 EM='The following reactions are empty (no involved metabolites):';
0282 dispEM(EM,false,model.rxns(~any(model.S,1)),trimWarnings);
0283 EM='The following metabolites are never used in a reaction:';
0284 dispEM(EM,false,model.mets(~any(model.S,2)),trimWarnings);
0285 if isfield(model,'genes')
0286     EM='The following genes are not associated to a reaction:';
0287     dispEM(EM,false,model.genes(~any(model.rxnGeneMat,1)),trimWarnings);
0288 end
0289 I=true(numel(model.comps),1);
0290 I(model.metComps)=false;
0291 EM='The following compartments contain no metabolites:';
0292 dispEM(EM,false,model.comps(I),trimWarnings);
0293 
0294 %Contradicting bounds
0295 EM='The following reactions have contradicting bounds (lower bound is higher than upper bound):';
0296 dispEM(EM,throwErrors,model.rxns(model.lb>model.ub),trimWarnings);
0297 EM='The following reactions have lower and upper bounds that indicate reversibility, but are indicated as irreversible in model.rev:';
0298 dispEM(EM,false,model.rxns(model.lb < 0 & model.ub > 0 & model.rev==0),trimWarnings);
0299 
0300 %Multiple or no objective functions not allowed in SBML L3V1 FBCv2
0301 if numel(find(model.c))>1
0302     EM='Multiple objective functions found. This might be intended, but results in FBCv2 non-compliant SBML file when exported';
0303     dispEM(EM,false,model.rxns(find(model.c)),trimWarnings);
0304 elseif ~any(model.c)
0305     EM='No objective function found. This might be intended, but results in FBCv2 non-compliant SBML file when exported';
0306     dispEM(EM,false);
0307 end
0308 
0309 %Mapping of compartments
0310 if isfield(model,'compOutside')
0311     EM='The following compartments are in "compOutside" but not in "comps":';
0312     dispEM(EM,throwErrors,setdiff(model.compOutside,[{''};model.comps]),trimWarnings);
0313 end
0314 
0315 %Met names which start with number
0316 I=false(numel(model.metNames),1);
0317 for i=1:numel(model.metNames)
0318     index=strfind(model.metNames{i},' ');
0319     if any(index)
0320         if any(str2double(model.metNames{i}(1:index(1)-1)))
0321             I(i)=true;
0322         end
0323     end
0324 end
0325 EM='The following metabolite names begin with a number directly followed by space, which could potentially cause problems:';
0326 dispEM(EM,false,model.metNames(I),trimWarnings);
0327 
0328 %Non-parseable composition
0329 if isfield(model,'metFormulas')
0330     [~, ~, exitFlag]=parseFormulas(model.metFormulas,true,false);
0331     EM='The composition for the following metabolites could not be parsed:';
0332     dispEM(EM,false,model.mets(exitFlag==-1),trimWarnings);
0333 end
0334 
0335 %Check if there are metabolites with different names but the same MIRIAM
0336 %codes
0337 if isfield(model,'metMiriams')
0338     miriams=containers.Map();
0339     for i=1:numel(model.mets)
0340         if ~isempty(model.metMiriams{i})
0341             %Loop through and add for each miriam
0342             for j=1:numel(model.metMiriams{i}.name)
0343                 %Get existing metabolite indexes
0344                 current=strcat(model.metMiriams{i}.name{j},'/',model.metMiriams{i}.value{j});
0345                 if isKey(miriams,current)
0346                     existing=miriams(current);
0347                 else
0348                     existing=[];
0349                 end
0350                 miriams(current)=[existing;i];
0351             end
0352         end
0353     end
0354     
0355     %Get all keys
0356     allMiriams=keys(miriams);
0357     
0358     hasMultiple=false(numel(allMiriams),1);
0359     for i=1:numel(allMiriams)
0360         if numel(miriams(allMiriams{i}))>1
0361             %Check if they all have the same name
0362             if numel(unique(model.metNames(miriams(allMiriams{i}))))>1
0363                 if ~regexp(allMiriams{i},'^sbo\/SBO:') % SBO terms are expected to be multiple
0364                     hasMultiple(i)=true;
0365                 end                
0366             end
0367         end
0368     end
0369     
0370     %Print output
0371     EM='The following MIRIAM strings are associated to more than one unique metabolite name:';
0372     dispEM(EM,false,allMiriams(hasMultiple),trimWarnings);
0373 end
0374 
0375 %Check if there are metabolites with different names but the same InChI
0376 %codes
0377 if isfield(model,'inchis')
0378     inchis=containers.Map();
0379     for i=1:numel(model.mets)
0380         if ~isempty(model.inchis{i})
0381             %Get existing metabolite indexes
0382             if isKey(inchis,model.inchis{i})
0383                 existing=inchis(model.inchis{i});
0384             else
0385                 existing=[];
0386             end
0387             inchis(model.inchis{i})=[existing;i];
0388         end
0389     end
0390     
0391     %Get all keys
0392     allInchis=keys(inchis);
0393     
0394     hasMultiple=false(numel(allInchis),1);
0395     for i=1:numel(allInchis)
0396         if numel(inchis(allInchis{i}))>1
0397             %Check if they all have the same name
0398             if numel(unique(model.metNames(inchis(allInchis{i}))))>1
0399                 hasMultiple(i)=true;
0400             end
0401         end
0402     end
0403     
0404     %Print output
0405     EM='The following InChI strings are associated to more than one unique metabolite name:';
0406     dispEM(EM,false,allInchis(hasMultiple),trimWarnings);
0407 end
0408 
0409 % %Check if there are metabolites with different names but the same SMILES
0410 % if isfield(model,'metSmiles')
0411 %     metSmiles=containers.Map();
0412 %     for i=1:numel(model.mets)
0413 %         if ~isempty(model.metSmiles{i})
0414 %             %Get existing metabolite indexes
0415 %             if isKey(metSmiles,model.metSmiles{i})
0416 %                 existing=metSmiles(model.metSmiles{i});
0417 %             else
0418 %                 existing=[];
0419 %             end
0420 %             metSmiles(model.metSmiles{i})=[existing;i];
0421 %         end
0422 %     end
0423 %
0424 %     %Get all keys
0425 %     allmetSmiles=keys(metSmiles);
0426 %
0427 %     hasMultiple=false(numel(metSmiles),1);
0428 %     for i=1:numel(metSmiles)
0429 %         if numel(metSmiles(metSmiles{i}))>1
0430 %             %Check if they all have the same name
0431 %             if numel(unique(model.metNames(metSmiles(allmetSmiles{i}))))>1
0432 %                 hasMultiple(i)=true;
0433 %             end
0434 %         end
0435 %     end
0436 %
0437 %     %Print output
0438 %     EM='The following metSmiles strings are associated to more than one unique metabolite name:';
0439 %     dispEM(EM,false,allmetSmiles(hasMultiple),trimWarnings);
0440 % end
0441 end
0442 
0443 function I=duplicates(strings)
0444 I=false(numel(strings),1);
0445 [J, K]=unique(strings);
0446 if numel(J)~=numel(strings)
0447     L=1:numel(strings);
0448     L(K)=[];
0449     I(L)=true;
0450 end
0451 end