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