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