importModel

 importModel
   Import a constraint-based model from an SBML file.

 Input:
   fileName        a SBML file to import. A dialog window will open if
                   no file name is specified.
   removeExcMets   true if exchange metabolites should be removed. This is
                   needed to be able to run simulations, but it could also
                   be done using simplifyModel at a later stage (optional,
                   default true)
   removePrefix    true if identifier prefixes should be removed when
                   loading the model: G_ for genes, R_ for reactions,
                   M_ for metabolites, and C_ for compartments. These are
                   only removed if all identifiers of a certain type
                   contain the prefix. (optional, default true)
   supressWarnings true if warnings regarding the model structure should
                   be supressed (optional, default false)

 Output:
   model
       id               model ID
       name             name of model contents
       annotation       additional information about model
       rxns             reaction ids
       mets             metabolite ids
       S                stoichiometric matrix
       lb               lower bounds
       ub               upper bounds
       rev              reversibility vector
       c                objective coefficients
       b                equality constraints for the metabolite equations
       comps            compartment ids
       compNames        compartment names
       compOutside      the id (as in comps) for the compartment
                        surrounding each of the compartments
       compMiriams      structure with MIRIAM information about the
                        compartments
       rxnNames         reaction description
       rxnComps         compartments for reactions
       grRules          reaction to gene rules in text form
       rxnGeneMat       reaction-to-gene mapping in sparse matrix form
       subSystems       subsystem name for each reaction
       eccodes          EC-codes for the reactions
       rxnMiriams       structure with MIRIAM information about the reactions
       rxnNotes         reaction notes
       rxnReferences    reaction references
       rxnConfidenceScores reaction confidence scores
       genes            list of all genes
       geneComps        compartments for genes
       geneMiriams      structure with MIRIAM information about the genes
       geneShortNames   gene alternative names (e.g. ERG10)
       proteins         protein associated to each gene
       metNames         metabolite description
       metComps         compartments for metabolites
       inchis           InChI-codes for metabolites
       metFormulas      metabolite chemical formula
       metMiriams       structure with MIRIAM information about the metabolites
       metCharges       metabolite charge
       unconstrained    true if the metabolite is an exchange metabolite

 Note: A number of consistency checks are performed in order to ensure that the
 model is valid. Take these warnings seriously and modify the model
 structure to solve them.

 Usage: model = importModel(fileName, removeExcMets, removePrefix, supressWarnings)

0001 function model=importModel(fileName,removeExcMets,removePrefix,supressWarnings)
0002 % importModel
0003 %   Import a constraint-based model from an SBML file.
0004 %
0005 % Input:
0006 %   fileName        a SBML file to import. A dialog window will open if
0007 %                   no file name is specified.
0008 %   removeExcMets   true if exchange metabolites should be removed. This is
0009 %                   needed to be able to run simulations, but it could also
0010 %                   be done using simplifyModel at a later stage (optional,
0011 %                   default true)
0012 %   removePrefix    true if identifier prefixes should be removed when
0013 %                   loading the model: G_ for genes, R_ for reactions,
0014 %                   M_ for metabolites, and C_ for compartments. These are
0015 %                   only removed if all identifiers of a certain type
0016 %                   contain the prefix. (optional, default true)
0017 %   supressWarnings true if warnings regarding the model structure should
0018 %                   be supressed (optional, default false)
0019 %
0020 % Output:
0021 %   model
0022 %       id               model ID
0023 %       name             name of model contents
0024 %       annotation       additional information about model
0025 %       rxns             reaction ids
0026 %       mets             metabolite ids
0027 %       S                stoichiometric matrix
0028 %       lb               lower bounds
0029 %       ub               upper bounds
0030 %       rev              reversibility vector
0031 %       c                objective coefficients
0032 %       b                equality constraints for the metabolite equations
0033 %       comps            compartment ids
0034 %       compNames        compartment names
0035 %       compOutside      the id (as in comps) for the compartment
0036 %                        surrounding each of the compartments
0037 %       compMiriams      structure with MIRIAM information about the
0038 %                        compartments
0039 %       rxnNames         reaction description
0040 %       rxnComps         compartments for reactions
0041 %       grRules          reaction to gene rules in text form
0042 %       rxnGeneMat       reaction-to-gene mapping in sparse matrix form
0043 %       subSystems       subsystem name for each reaction
0044 %       eccodes          EC-codes for the reactions
0045 %       rxnMiriams       structure with MIRIAM information about the reactions
0046 %       rxnNotes         reaction notes
0047 %       rxnReferences    reaction references
0048 %       rxnConfidenceScores reaction confidence scores
0049 %       genes            list of all genes
0050 %       geneComps        compartments for genes
0051 %       geneMiriams      structure with MIRIAM information about the genes
0052 %       geneShortNames   gene alternative names (e.g. ERG10)
0053 %       proteins         protein associated to each gene
0054 %       metNames         metabolite description
0055 %       metComps         compartments for metabolites
0056 %       inchis           InChI-codes for metabolites
0057 %       metFormulas      metabolite chemical formula
0058 %       metMiriams       structure with MIRIAM information about the metabolites
0059 %       metCharges       metabolite charge
0060 %       unconstrained    true if the metabolite is an exchange metabolite
0061 %
0062 % Note: A number of consistency checks are performed in order to ensure that the
0063 % model is valid. Take these warnings seriously and modify the model
0064 % structure to solve them.
0065 %
0066 % Usage: model = importModel(fileName, removeExcMets, removePrefix, supressWarnings)
0067 
0068 if nargin<1 || isempty(fileName)
0069     [fileName, pathName] = uigetfile({'*.xml;*.sbml'}, 'Please select the model file');
0070     if fileName == 0
0071         error('You should select a model file')
0072     else
0073         fileName = fullfile(pathName,fileName);
0074     end
0075 end
0076 fileName=char(fileName);
0077 if nargin<2 || isempty(removeExcMets)
0078     removeExcMets=true;
0079 end
0080 
0081 if nargin<3 || isempty(removePrefix)
0082     removePrefix=true;
0083 end
0084 
0085 if nargin<4
0086     supressWarnings=false;
0087 end
0088 
0089 fileName=checkFileExistence(fileName,1);
0090 % If path contains non-ASCII characters, copy file to tempdir first, as
0091 % libSBML is known to have problems with this on Windows:
0092 % https://sbml.org/software/libsbml/libsbml-docs/known-pitfalls/#matlab-on-windows-has-issues-with-unicode-filenames
0093 if ispc && any(double(fileName)>128)
0094     [~,originalFile,ext] = fileparts(fileName);
0095     tempFile = fullfile(tempdir,[originalFile ext]);
0096     copyfile(fileName,tempFile);
0097     fileName = tempFile;
0098 end
0099 
0100 %This is to match the order of the fields to those you get from importing
0101 %from Excel
0102 model=[];
0103 model.id=[];
0104 model.name=[];
0105 model.annotation=[];
0106 model.rxns={};
0107 model.mets={};
0108 model.S=[];
0109 model.lb=[];
0110 model.ub=[];
0111 model.rev=[];
0112 model.c=[];
0113 model.b=[];
0114 model.comps={};
0115 model.compNames={};
0116 model.compOutside={};
0117 model.compMiriams={};
0118 model.rxnNames={};
0119 model.rxnComps=[];
0120 model.grRules={};
0121 model.rxnGeneMat=[];
0122 model.subSystems={};
0123 model.eccodes={};
0124 model.rxnMiriams={};
0125 model.rxnNotes={};
0126 model.rxnReferences={};
0127 model.rxnConfidenceScores=[];
0128 model.genes={};
0129 model.geneComps=[];
0130 model.geneMiriams={};
0131 model.geneShortNames={};
0132 model.proteins={};
0133 model.metNames={};
0134 model.metComps=[];
0135 model.inchis={};
0136 model.metFormulas={};
0137 model.metMiriams={};
0138 model.metCharges=[];
0139 model.unconstrained=[];
0140 
0141 %Load the model using libSBML
0142 [modelSBML,errorMsg] = TranslateSBML_RAVEN(fileName,0,0,[1 1]);
0143 if exist('tempFile','var')
0144     delete(tempFile)
0145 end
0146 
0147 if isempty(modelSBML)
0148     EM=['There is a problem with the SBML file. Try using the SBML Validator at http://sbml.org/Facilities/Validator.\nlibSBML reports: ', errorMsg.message];
0149     dispEM(EM);
0150 end
0151 
0152 %Retrieve compartment names and IDs
0153 compartmentNames=cell(numel(modelSBML.compartment),1);
0154 compartmentIDs=cell(numel(modelSBML.compartment),1);
0155 compartmentOutside=cell(numel(modelSBML.compartment),1);
0156 compartmentMiriams=cell(numel(modelSBML.compartment),1);
0157 
0158 if isfield(modelSBML.compartment,'sboTerm') && isscalar(unique([modelSBML.compartment.sboTerm]))
0159     %If all the SBO terms are identical, don't add them to compMiriams
0160     modelSBML.compartment = rmfield(modelSBML.compartment,'sboTerm');
0161 end
0162 
0163 for i=1:numel(modelSBML.compartment)
0164     compartmentNames{i}=modelSBML.compartment(i).name;
0165     compartmentIDs{i}=modelSBML.compartment(i).id;
0166     if isfield(modelSBML.compartment(i),'outside')
0167         if ~isempty(modelSBML.compartment(i).outside)
0168             compartmentOutside{i}=modelSBML.compartment(i).outside;
0169         else
0170             compartmentOutside{i}='';
0171         end
0172     else
0173         compartmentOutside{i}=[];
0174     end
0175 
0176     if isfield(modelSBML.compartment(i),'annotation')
0177         compartmentMiriams{i}=parseMiriam(modelSBML.compartment(i).annotation);
0178     else
0179         compartmentMiriams{i}=[];
0180     end
0181 
0182     if isfield(modelSBML.compartment(i),'sboTerm') && ~(modelSBML.compartment(i).sboTerm==-1)
0183         compartmentMiriams{i} = addSBOtoMiriam(compartmentMiriams{i},modelSBML.compartment(i).sboTerm);
0184     end
0185 end
0186 
0187 %If there are no compartment names then use compartment id as name
0188 if all(cellfun(@isempty,compartmentNames))
0189     compartmentNames=compartmentIDs;
0190 end
0191 
0192 %Retrieve info on metabolites, genes, complexes
0193 metaboliteNames={};
0194 metaboliteIDs={};
0195 metaboliteCompartments={};
0196 metaboliteUnconstrained=[];
0197 metaboliteFormula={};
0198 metaboliteInChI={};
0199 metaboliteMiriams={};
0200 metaboliteCharges=[];
0201 
0202 geneNames={};
0203 geneIDs={};
0204 geneMiriams={};
0205 geneShortNames={};
0206 proteins={};
0207 geneCompartments={};
0208 complexIDs={};
0209 complexNames={};
0210 
0211 %If the file is not a COBRA Toolbox model. According to the format
0212 %specified in the yeast consensus model both metabolites and genes are a
0213 %type of 'species'. The metabolites have names starting with 'M_' and genes
0214 %with 'E_'
0215 geneSBOs = [];
0216 metSBOs = [];
0217 %Regex of compartment names, later to be used to remove from metabolite
0218 %names if present as suffix.
0219 regexCompNames = ['\s?\[((' strjoin({modelSBML.compartment.name},')|(') '))\]$'];
0220 for i=1:numel(modelSBML.species)
0221     if length(modelSBML.species(i).id)>=2 && strcmpi(modelSBML.species(i).id(1:2),'E_')
0222         geneNames{numel(geneNames)+1,1}=modelSBML.species(i).name;
0223 
0224         %The "E_" is included in the ID. This is because it's only used
0225         %internally in this file and it makes the matching a little
0226         %smoother
0227         geneIDs{numel(geneIDs)+1,1}=modelSBML.species(i).id;
0228         geneCompartments{numel(geneCompartments)+1,1}=modelSBML.species(i).compartment;
0229 
0230         %Get Miriam structure
0231         if isfield(modelSBML.species(i),'annotation')
0232             %Get Miriam info
0233             geneMiriam=parseMiriam(modelSBML.species(i).annotation);
0234             geneMiriams{numel(geneMiriams)+1,1}=geneMiriam;
0235         else
0236             geneMiriams{numel(geneMiriams)+1,1}=[];
0237         end
0238 
0239         %Protein short names (for example ERG10) are saved as SHORT
0240         %NAME: NAME in the notes-section of metabolites for SBML Level
0241         %2 and as PROTEIN_ASSOCIATION for each reaction in SBML Level 2
0242         %COBRA Toolbox format. For now only the SHORT NAME is loaded
0243         %and no mapping takes place
0244         if isfield(modelSBML.species(i),'notes')
0245             geneShortNames{numel(geneShortNames)+1,1}=parseNote(modelSBML.species(i).notes,'SHORT NAME');
0246         else
0247             geneShortNames{numel(geneShortNames)+1,1}='';
0248         end
0249 
0250         %Get SBO term
0251         if isfield(modelSBML.species(i),'sboTerm') && ~(modelSBML.species(i).sboTerm==-1)
0252             geneSBOs(end+1,1) = modelSBML.species(i).sboTerm;
0253         end
0254     elseif length(modelSBML.species(i).id)>=2 && strcmpi(modelSBML.species(i).id(1:3),'Cx_')
0255         %If it's a complex keep the ID and name
0256         complexIDs=[complexIDs;modelSBML.species(i).id];
0257         complexNames=[complexNames;modelSBML.species(i).name];
0258     else
0259         %If it is not gene or complex, then it must be a metabolite
0260         metaboliteNames{numel(metaboliteNames)+1,1}=modelSBML.species(i).name;
0261         metaboliteIDs{numel(metaboliteIDs)+1,1}=modelSBML.species(i).id;
0262         metaboliteCompartments{numel(metaboliteCompartments)+1,1}=modelSBML.species(i).compartment;
0263         metaboliteUnconstrained(numel(metaboliteUnconstrained)+1,1)=modelSBML.species(i).boundaryCondition;
0264 
0265         %For each metabolite retrieve the formula and the InChI code if
0266         %available First add the InChI code and the formula from the
0267         %InChI. This allows for overwriting the formula by setting the
0268         %actual formula field
0269         if ~isempty(modelSBML.species(i).annotation)
0270             %Get the formula if available
0271             startString='>InChI=';
0272             endString='</in:inchi>';
0273             formStart=strfind(modelSBML.species(i).annotation,startString);
0274             if isempty(formStart)
0275                 startString='InChI=';
0276                 endString='"/>';
0277             end
0278             formStart=strfind(modelSBML.species(i).annotation,startString);
0279             if ~isempty(formStart)
0280                 formEnd=strfind(modelSBML.species(i).annotation,endString);
0281                 formEndIndex=find(formEnd>formStart, 1 );
0282                 formula=modelSBML.species(i).annotation(formStart+numel(startString):formEnd(formEndIndex)-1);
0283                 metaboliteInChI{numel(metaboliteInChI)+1,1}=formula;
0284 
0285                 %The composition is most often present between the
0286                 %first and second "/" in the model. In some simple
0287                 %molecules, such as salts, there is no second "/". The
0288                 %formula is then assumed to be to the end of the string
0289                 compositionIndexes=strfind(formula,'/');
0290                 if numel(compositionIndexes)>1
0291                     metaboliteFormula{numel(metaboliteFormula)+1,1}=...
0292                         formula(compositionIndexes(1)+1:compositionIndexes(2)-1);
0293                 else
0294                     if isscalar(compositionIndexes)
0295                         %Probably a simple molecule which can have only
0296                         %one conformation
0297                         metaboliteFormula{numel(metaboliteFormula)+1,1}=...
0298                             formula(compositionIndexes(1)+1:numel(formula));
0299                     else
0300                         metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0301                     end
0302                 end
0303             elseif isfield(modelSBML.species(i),'fbc_chemicalFormula')
0304                 metaboliteInChI{numel(metaboliteInChI)+1,1}='';
0305                 if ~isempty(modelSBML.species(i).fbc_chemicalFormula)
0306                     %Cannot extract InChi from formula, so remains
0307                     %empty
0308                     metaboliteFormula{numel(metaboliteFormula)+1,1}=modelSBML.species(i).fbc_chemicalFormula;
0309                 else
0310                     metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0311                 end
0312             else
0313                 metaboliteInChI{numel(metaboliteInChI)+1,1}='';
0314                 metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0315             end
0316 
0317             %Get Miriam info
0318             metMiriam=parseMiriam(modelSBML.species(i).annotation);
0319             metaboliteMiriams{numel(metaboliteMiriams)+1,1}=metMiriam;
0320         else
0321             metaboliteInChI{numel(metaboliteInChI)+1,1}='';
0322             if isfield(modelSBML.species(i),'notes')
0323                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
0324             else
0325                 metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0326             end
0327             metaboliteMiriams{numel(metaboliteMiriams)+1,1}=[];
0328         end
0329         if ~isempty(modelSBML.species(i).notes)
0330             if ~isfield(modelSBML.species(i),'annotation')
0331                 metaboliteFormula{numel(metaboliteFormula)+1,1}=parseNote(modelSBML.species(i).notes,'FORMULA');
0332             end
0333         elseif ~isfield(modelSBML.species(i),'annotation')
0334             metaboliteFormula{numel(metaboliteFormula)+1,1}='';
0335         end
0336         %Get SBO term
0337         if isfield(modelSBML.species(i),'sboTerm') && ~(modelSBML.species(i).sboTerm==-1)
0338             metSBOs(end+1,1) = modelSBML.species(i).sboTerm;
0339         end
0340     end
0341 
0342     %The following lines are executed regardless isSBML2COBRA setting
0343     if isempty(modelSBML.species(i).id) || ~strcmpi(modelSBML.species(i).id(1:2),'E_')
0344         if isempty(modelSBML.species(i).id) || ~strcmpi(modelSBML.species(i).id(1:3),'Cx_')
0345             %Remove trailing [compartment] from metabolite name if present
0346             metaboliteNames{end,1}=regexprep(metaboliteNames{end,1},regexCompNames,'');
0347             metaboliteNames{end,1}=metaboliteNames{end,1};
0348             if isfield(modelSBML.species(i),'fbc_charge')
0349                 if ~isempty(modelSBML.species(i).fbc_charge) && modelSBML.species(i).isSetfbc_charge
0350                     metaboliteCharges(numel(metaboliteCharges)+1,1)=double(modelSBML.species(i).fbc_charge);
0351                 else
0352                     if isfield(modelSBML.species(i),'notes')
0353                         if strfind(modelSBML.species(i).notes,'CHARGE')
0354                             metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
0355                         else
0356                             metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0357                         end
0358                     else
0359                         metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0360                     end
0361                 end
0362             elseif isfield(modelSBML.species(i),'notes')
0363                 if strfind(modelSBML.species(i).notes,'CHARGE')
0364                     metaboliteCharges(numel(metaboliteCharges)+1,1)=str2double(parseNote(modelSBML.species(i).notes,'CHARGE'));
0365                 else
0366                     metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0367                 end
0368             else
0369                 metaboliteCharges(numel(metaboliteCharges)+1,1)=NaN;
0370             end
0371             %Additional information from FBC format Chemical formula
0372             if isfield(modelSBML.species(i),'fbc_chemicalFormula')
0373                 if ~isempty(modelSBML.species(i).fbc_chemicalFormula)
0374                     metaboliteFormula{numel(metaboliteFormula),1}=modelSBML.species(i).fbc_chemicalFormula;
0375                 end
0376             end
0377         end
0378     end
0379 end
0380 
0381 %Add SBO terms to gene and metabolite miriam fields
0382 if numel(unique(geneSBOs)) > 1  % don't add if they're all identical
0383     for i = 1:numel(geneNames)
0384         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},geneSBOs(i));
0385     end
0386 end
0387 if numel(unique(metSBOs)) > 1
0388     for i = 1:numel(metaboliteNames)
0389         metaboliteMiriams{i} = addSBOtoMiriam(metaboliteMiriams{i},metSBOs(i));
0390     end
0391 end
0392 
0393 %Retrieve info on reactions
0394 reactionNames=cell(numel(modelSBML.reaction),1);
0395 reactionIDs=cell(numel(modelSBML.reaction),1);
0396 subsystems=cell(numel(modelSBML.reaction),1);
0397 eccodes=cell(numel(modelSBML.reaction),1);
0398 eccodes(:,:)=cellstr('');
0399 rxnconfidencescores=NaN(numel(modelSBML.reaction),1);
0400 rxnreferences=cell(numel(modelSBML.reaction),1);
0401 rxnreferences(:,:)=cellstr('');
0402 rxnnotes=cell(numel(modelSBML.reaction),1);
0403 rxnnotes(:,:)=cellstr('');
0404 grRules=cell(numel(modelSBML.reaction),1);
0405 grRules(:,:)=cellstr('');
0406 grRulesFromModifier=grRules;
0407 rxnComps=zeros(numel(modelSBML.reaction),1);
0408 rxnMiriams=cell(numel(modelSBML.reaction),1);
0409 reactionReversibility=zeros(numel(modelSBML.reaction),1);
0410 reactionUB=zeros(numel(modelSBML.reaction),1);
0411 reactionLB=zeros(numel(modelSBML.reaction),1);
0412 reactionObjective=zeros(numel(modelSBML.reaction),1);
0413 
0414 %Construct the stoichiometric matrix while the reaction info is read
0415 S=zeros(numel(metaboliteIDs),numel(modelSBML.reaction));
0416 
0417 counter=0;
0418 %If FBC, then bounds have parameter ids defined for the whole model
0419 if isfield(modelSBML,'parameter')
0420     parameter.name=cell(numel(modelSBML.parameter),1);
0421     parameter.name={modelSBML.parameter(:).id}';
0422     parameter.value={modelSBML.parameter(:).value}';
0423 end
0424 
0425 if isfield(modelSBML.reaction,'sboTerm') && isscalar(unique([modelSBML.reaction.sboTerm]))
0426     %If all the SBO terms are identical, don't add them to rxnMiriams
0427     modelSBML.reaction = rmfield(modelSBML.reaction,'sboTerm');
0428 end
0429 
0430 for i=1:numel(modelSBML.reaction)
0431 
0432     %Check that the reaction doesn't produce a complex and nothing else. If
0433     %so, then jump to the next reaction. This is because I get the genes
0434     %for complexes from the names and not from the reactions that create
0435     %them. This only applies to the non-COBRA format
0436     if isscalar(modelSBML.reaction(i).product)
0437         if length(modelSBML.reaction(i).product(1).species)>=3
0438             if strcmp(modelSBML.reaction(i).product(1).species(1:3),'Cx_')==true
0439                 continue;
0440             end
0441         end
0442     end
0443 
0444     %It didn't look like a gene complex-forming reaction
0445     counter=counter+1;
0446 
0447     reactionNames{counter}=modelSBML.reaction(i).name;
0448 
0449     reactionIDs{counter}=modelSBML.reaction(i).id;
0450     reactionReversibility(counter)=modelSBML.reaction(i).reversible;
0451 
0452     %If model is FBC, first get parameter of bound and then replace it with
0453     %the correct value. Probably faster with replace(), but this was only
0454     %introduced in Matlab R2016b
0455     if isfield(modelSBML.reaction(i),'fbc_lowerFluxBound')
0456         lb=modelSBML.reaction(i).fbc_lowerFluxBound;
0457         ub=modelSBML.reaction(i).fbc_upperFluxBound;
0458         [~,fluxBoundIdx] = ismember({lb,ub},parameter.name);
0459         reactionLB(counter) = parameter.value{fluxBoundIdx(1)};
0460         reactionUB(counter) = parameter.value{fluxBoundIdx(2)};
0461     elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter')
0462         reactionLB(counter)=modelSBML.reaction(i).kineticLaw.parameter(1).value;
0463         reactionUB(counter)=modelSBML.reaction(i).kineticLaw.parameter(2).value;
0464         reactionObjective(counter)=modelSBML.reaction(i).kineticLaw.parameter(3).value;
0465     else
0466         if reactionReversibility(counter)==true
0467             reactionLB(counter)=-inf;
0468         else
0469             reactionLB(counter)=0;
0470         end
0471         reactionUB(counter)=inf;
0472         reactionObjective(counter)=0;
0473     end
0474 
0475     %Find the associated gene if available
0476     %If FBC, get gene association data from corresponding fields
0477     if isfield(modelSBML.reaction(i),'fbc_geneProductAssociation')
0478         if ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation) && ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association)
0479             grRules{counter}=modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association;
0480         end
0481     elseif isfield(modelSBML.reaction(i),'notes')
0482         %This section was previously executed only if isSBML2COBRA is true. Now
0483         %it will be executed, if 'GENE_ASSOCIATION' is found in
0484         %modelSBML.reaction(i).notes
0485         if strfind(modelSBML.reaction(i).notes,'GENE_ASSOCIATION')
0486             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION');
0487         elseif strfind(modelSBML.reaction(i).notes,'GENE ASSOCIATION')
0488             geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION');
0489         else
0490             geneAssociation='';
0491         end
0492         if ~isempty(geneAssociation)
0493             %This adds the grRules. The gene list and rxnGeneMat are created
0494             %later
0495             grRules{counter}=geneAssociation;
0496         end
0497     end
0498     if isempty(grRules{counter}) && ~isempty(modelSBML.reaction(i).modifier)
0499         rules='';
0500         for j=1:numel(modelSBML.reaction(i).modifier)
0501             modifier=modelSBML.reaction(i).modifier(j).species;
0502             if ~isempty(modifier)
0503                 if strcmpi(modifier(1:2),'E_')
0504                     index=find(strcmp(modifier,geneIDs));
0505                     %This should be unique and in the geneIDs list,
0506                     %otherwise something is wrong
0507                     if numel(index)~=1
0508                         EM=['Could not get the gene association data from reaction ' reactionIDs{i}];
0509                         dispEM(EM);
0510                     end
0511                     if ~isempty(rules)
0512                         rules=[rules ' or (' geneNames{index} ')'];
0513                     else
0514                         rules=['(' geneNames{index} ')'];
0515                     end
0516                 elseif strcmp(modifier(1:2),'s_')
0517                     index=find(strcmp(modifier,metaboliteIDs));
0518                     %This should be unique and in the geneIDs list,
0519                     %otherwise something is wrong
0520                     if numel(index)~=1
0521                         EM=['Could not get the gene association data from reaction ' reactionIDs{i}];
0522                         dispEM(EM);
0523                     end
0524                     if ~isempty(rules)
0525                         rules=[rules ' or (' metaboliteIDs{index} ')'];
0526                     else
0527                         rules=['(' metaboliteIDs{index} ')'];
0528                     end
0529                 else
0530                     %It seems to be a complex. Add the corresponding
0531                     %genes from the name of the complex (not the
0532                     %reaction that creates it)
0533                     index=find(strcmp(modifier,complexIDs));
0534                     if isscalar(index)
0535                         if ~isempty(rules)
0536                             rules=[rules ' or (' strrep(complexNames{index},':',' and ') ')'];
0537                         else
0538                             rules=['(' strrep(complexNames{index},':',' and ') ')'];
0539                         end
0540                     else
0541                         %Could not find a complex
0542                         EM=['Could not get the gene association data from reaction ' reactionIDs{i}];
0543                         dispEM(EM);
0544                     end
0545                 end
0546             end
0547         end
0548         grRules{counter}=rules;
0549         grRulesFromModifier{counter}=rules;%Backup copy for grRules, useful to parse Yeast 7.6
0550     end
0551 
0552     %Add reaction compartment
0553     if isfield(modelSBML.reaction(i),'compartment')
0554         if ~isempty(modelSBML.reaction(i).compartment)
0555             rxnComp=modelSBML.reaction(i).compartment;
0556         else
0557             rxnComp='';
0558         end
0559     elseif isfield(modelSBML.reaction(i),'notes')
0560         rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT');
0561     end
0562     if ~isempty(rxnComp)
0563         %Find it in the compartment list
0564         [~, J]=ismember(rxnComp,compartmentIDs);
0565         rxnComps(counter)=J;
0566     end
0567 
0568 
0569     miriamStruct=parseMiriam(modelSBML.reaction(i).annotation);
0570     rxnMiriams{counter}=miriamStruct;
0571     if isfield(modelSBML.reaction(i),'notes')
0572         subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM'));
0573         subsystems{counter,1}(cellfun('isempty',subsystems{counter,1})) = [];
0574         if strfind(modelSBML.reaction(i).notes,'Confidence Level')
0575             confScore = parseNote(modelSBML.reaction(i).notes,'Confidence Level');
0576             if isempty(confScore)
0577                 confScore = 0;
0578             end
0579             rxnconfidencescores(counter)=str2double(confScore);
0580         end
0581         rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS');
0582         rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES');
0583     end
0584 
0585     %Get SBO terms
0586     if isfield(modelSBML.reaction(i),'sboTerm') && ~(modelSBML.reaction(i).sboTerm==-1)
0587         rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm);
0588     end
0589 
0590     %Get ec-codes
0591     eccode='';
0592     if ~isempty(modelSBML.reaction(i).annotation)
0593         if strfind(modelSBML.reaction(i).annotation,'urn:miriam:ec-code')
0594             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code');
0595         elseif strfind(modelSBML.reaction(i).annotation,'http://identifiers.org/ec-code')
0596             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code');
0597         elseif strfind(modelSBML.reaction(i).annotation,'https://identifiers.org/ec-code')
0598             eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code');
0599         end
0600     elseif isfield(modelSBML.reaction(i),'notes')
0601         if strfind(modelSBML.reaction(i).notes,'EC Number')
0602             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')];
0603         elseif strfind(modelSBML.reaction(i).notes,'PROTEIN_CLASS')
0604             eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')];
0605         end
0606     end
0607     eccodes{counter}=eccode;
0608 
0609     %Add all reactants
0610     for j=1:numel(modelSBML.reaction(i).reactant)
0611         %Get the index of the metabolite in metaboliteIDs. External
0612         %metabolites will be removed at a later stage
0613         metIndex=find(strcmp(modelSBML.reaction(i).reactant(j).species,metaboliteIDs),1);
0614         if isempty(metIndex)
0615             EM=['Could not find metabolite ' modelSBML.reaction(i).reactant(j).species ' in reaction ' reactionIDs{counter}];
0616             dispEM(EM);
0617         end
0618         S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).reactant(j).stoichiometry*-1;
0619     end
0620 
0621     %Add all products
0622     for j=1:numel(modelSBML.reaction(i).product)
0623         %Get the index of the metabolite in metaboliteIDs.
0624         metIndex=find(strcmp(modelSBML.reaction(i).product(j).species,metaboliteIDs),1);
0625         if isempty(metIndex)
0626             EM=['Could not find metabolite ' modelSBML.reaction(i).product(j).species ' in reaction ' reactionIDs{counter}];
0627             dispEM(EM);
0628         end
0629         S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).product(j).stoichiometry;
0630     end
0631 end
0632 
0633 %if FBC, objective function is separately defined. Multiple objective
0634 %functions can be defined, one is set as active
0635 if isfield(modelSBML, 'fbc_activeObjective')
0636     obj=modelSBML.fbc_activeObjective;
0637     for i=1:numel(modelSBML.fbc_objective)
0638         if strcmp(obj,modelSBML.fbc_objective(i).fbc_id)
0639             if ~isempty(modelSBML.fbc_objective(i).fbc_fluxObjective)
0640                 rxn=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction;
0641                 idx=ismember(reactionIDs,rxn);
0642                 reactionObjective(idx)=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient;
0643             end
0644         end
0645     end
0646 end
0647 
0648 %subSystems can be stored as groups instead of in annotations
0649 if isfield(modelSBML,'groups_group')
0650     for i=1:numel(modelSBML.groups_group)
0651         groupreactions={modelSBML.groups_group(i).groups_member(:).groups_idRef};
0652         [~, idx] = ismember(groupreactions, reactionIDs);
0653         if any(idx)
0654             for j=1:numel(idx)
0655                 if isempty(subsystems{idx(j)}) % First subsystem
0656                     subsystems{idx(j)} = {modelSBML.groups_group(i).groups_name};
0657                 else % Consecutive subsystems: concatenate
0658                     subsystems{idx(j)} = horzcat(subsystems{idx(j)}, modelSBML.groups_group(i).groups_name);
0659                 end
0660             end
0661         end
0662     end
0663 end
0664 
0665 %Shrink the structures if complex-forming reactions had to be skipped
0666 reactionNames=reactionNames(1:counter);
0667 reactionIDs=reactionIDs(1:counter);
0668 subsystems=subsystems(1:counter);
0669 eccodes=eccodes(1:counter);
0670 rxnconfidencescores=rxnconfidencescores(1:counter);
0671 rxnreferences=rxnreferences(1:counter);
0672 rxnnotes=rxnnotes(1:counter);
0673 grRules=grRules(1:counter);
0674 rxnMiriams=rxnMiriams(1:counter);
0675 reactionReversibility=reactionReversibility(1:counter);
0676 reactionUB=reactionUB(1:counter);
0677 reactionLB=reactionLB(1:counter);
0678 reactionObjective=reactionObjective(1:counter);
0679 S=S(:,1:counter);
0680 
0681 model.name=modelSBML.name;
0682 model.id=modelSBML.id;
0683 model.rxns=reactionIDs;
0684 model.mets=metaboliteIDs;
0685 model.S=sparse(S);
0686 model.lb=reactionLB;
0687 model.ub=reactionUB;
0688 model.rev=reactionReversibility;
0689 model.c=reactionObjective;
0690 model.b=zeros(numel(metaboliteIDs),1);
0691 model.comps=compartmentIDs;
0692 model.compNames=compartmentNames;
0693 model.rxnConfidenceScores=rxnconfidencescores;
0694 model.rxnReferences=rxnreferences;
0695 model.rxnNotes=rxnnotes;
0696 
0697 %Load annotation if available. If there are several authors, only the first
0698 %author credentials are imported
0699 if isfield(modelSBML,'annotation')
0700     endString='</';
0701     I=strfind(modelSBML.annotation,endString);
0702     J=strfind(modelSBML.annotation,'<vCard:Family>');
0703     if any(J)
0704         model.annotation.familyName=modelSBML.annotation(J(1)+14:I(find(I>J(1),1))-1);
0705     end
0706     J=strfind(modelSBML.annotation,'<vCard:Given>');
0707     if any(J)
0708         model.annotation.givenName=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1);
0709     end
0710     J=strfind(modelSBML.annotation,'<vCard:EMAIL>');
0711     if any(J)
0712         model.annotation.email=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1);
0713     end
0714     J=strfind(modelSBML.annotation,'<vCard:Orgname>');
0715     if any(J)
0716         model.annotation.organization=modelSBML.annotation(J(1)+15:I(find(I>J(1),1))-1);
0717     end
0718     endString='"/>';
0719     I=strfind(modelSBML.annotation,endString);
0720     if strfind(modelSBML.annotation,'"urn:miriam:')
0721         J=strfind(modelSBML.annotation,'"urn:miriam:');
0722         if any(J)
0723             model.annotation.taxonomy=modelSBML.annotation(J+12:I(find(I>J,1))-1);
0724         end
0725     else
0726         J=strfind(modelSBML.annotation,'"http://identifiers.org/');
0727         if any(J)
0728             I = I(find(I>J,1))-1;
0729             J = J(find(J<I,1))+24;
0730             model.annotation.taxonomy=modelSBML.annotation(J:I);
0731         else
0732             J=strfind(modelSBML.annotation,'"https://identifiers.org/');
0733             if any(J)
0734                 I = I(find(I>J,1))-1;
0735                 J = J(find(J<I,1))+25;
0736                 model.annotation.taxonomy=modelSBML.annotation(J:I);
0737             end
0738         end
0739     end
0740 end
0741 if isfield(modelSBML,'notes')
0742     %Two shapes reach here. RAVEN's own exportModel writes the note as bare
0743     %text inside a <body>; a writer that serialises a notes dictionary emits
0744     %one '<p>key: value</p>' per entry and no <body> at all, carrying the note
0745     %under its own key. Take the labelled entry when there is one, so the
0746     %other entries (version, and whatever else the dictionary held) stay out
0747     %of the note, and fall back to the <body> text otherwise.
0748     note=parseNote(modelSBML.notes,'note');
0749     if isempty(note)
0750         startString=strfind(modelSBML.notes,'xhtml">');
0751         endString=strfind(modelSBML.notes,'</body>');
0752         if any(startString) && any(endString)
0753             note=modelSBML.notes(startString(1)+7:endString-1);
0754         end
0755     end
0756     if ~isempty(note)
0757         model.annotation.note=regexprep(note,'<p.*?>|</p.*?>','');
0758         model.annotation.note=strtrim(model.annotation.note);
0759         if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML')
0760             model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel
0761         end
0762     end
0763 end
0764 
0765 if any(~cellfun(@isempty,compartmentOutside))
0766     model.compOutside=compartmentOutside;
0767 end
0768 
0769 model.rxnNames=reactionNames;
0770 model.metNames=metaboliteNames;
0771 
0772 %Match the compartments for metabolites
0773 [~, J]=ismember(metaboliteCompartments,model.comps);
0774 model.metComps=J;
0775 
0776 %If any genes have been loaded (only for the new format)
0777 if ~isempty(geneNames)
0778     %In some rare cases geneNames may not necessarily be used in grRules.
0779     %That is true for Yeast 7.6. It's therefore important to change gene
0780     %systematic names to geneIDs in sophisticated way. Gene systematic
0781     %names are not unique, since exactly the same name may be in different
0782     %compartments
0783     if all(cellfun(@isempty,strfind(grRules,geneNames{1})))
0784         geneShortNames=geneNames;
0785         %geneShortNames contain compartments as well, so these are removed
0786         geneShortNames=regexprep(geneShortNames,' \[.+$','');
0787         %grRules obtained from modifier fields contain geneNames. These are
0788         %changed into geneIDs. grRulesFromModifier is a good way to have
0789         %geneIDs and rxns association when it's important to resolve
0790         %systematic name ambiguities
0791         grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs);
0792         grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs);
0793 
0794         %Yeast 7.6 contains several metabolites, which were used in gene
0795         %associations. For that reason, the list of species ID is created
0796         %and we then check whether any of them have kegg.genes annotation
0797         %thereby obtaining systematic gene names
0798         geneShortNames=vertcat(geneShortNames,metaboliteNames);
0799         geneIDs=vertcat(geneIDs,metaboliteIDs);
0800         geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes');
0801         geneCompartments=vertcat(geneCompartments,metaboliteCompartments);
0802         geneMiriams=vertcat(geneMiriams,metaboliteMiriams);
0803 
0804         %Now we retain information for only these entries, which have
0805         %kegg.genes annotation
0806         geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames));
0807         geneIDs=geneIDs(~cellfun('isempty',geneSystNames));
0808         geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames));
0809         geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames));
0810         geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames));
0811         %Now we reorder geneIDs and geneSystNames by geneSystNames string
0812         %length
0813         geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later
0814         [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend');
0815         geneIDs = geneIDs(Indx);
0816         geneSystNames = geneSystNames(Indx);
0817         for i=1:numel(geneSystNames)
0818             for j=1:numel(grRules)
0819                 if strfind(grRules{j},geneSystNames{i})
0820                     if ~isempty(grRules{j})
0821                         if sum(ismember(geneSystNames,geneSystNames{i}))==1
0822                             grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i});
0823                         elseif sum(ismember(geneSystNames,geneSystNames{i}))>1
0824                             counter=0;
0825                             ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i}));
0826                             for k=1:numel(ovrlpIDs)
0827                                 if strfind(grRulesFromModifier{j},ovrlpIDs{k})
0828                                     counter=counter+1;
0829                                     grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k});
0830                                 end
0831                                 if counter>1
0832                                     EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id];
0833                                     dispEM(EM);
0834                                 end
0835                             end
0836                         end
0837                     end
0838                 end
0839             end
0840         end
0841     end
0842     model.genes=geneNames;
0843     model.grRules=grRules;
0844     [grRules,rxnGeneMat] = standardizeGrRules(model,true);
0845     model.grRules = grRules;
0846     model.rxnGeneMat = rxnGeneMat;
0847 
0848     %Match the compartments for genes
0849     [~, J]=ismember(geneCompartments,model.comps);
0850     model.geneComps=J;
0851 else
0852     if ~all(cellfun(@isempty,grRules))
0853         %If fbc_geneProduct exists, follow the specified gene order, such
0854         %that matching geneShortNames in function below will work
0855         if isfield(modelSBML,'fbc_geneProduct')
0856             genes={modelSBML.fbc_geneProduct.fbc_id};
0857 
0858             %Get gene Miriams if they were not retrieved above (this occurs
0859             %when genes are stored as fbc_geneProduct instead of species)
0860             if isempty(geneMiriams)
0861                 geneMiriams = cell(numel(genes),1);
0862                 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && isscalar(unique([modelSBML.fbc_geneProduct.sboTerm]))
0863                     %If all the SBO terms are identical, don't add them to geneMiriams
0864                     modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm');
0865                 end
0866                 for i = 1:numel(genes)
0867                     geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation);
0868                     if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1)
0869                         geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm);
0870                     end
0871                 end
0872             end
0873             proteins={modelSBML.fbc_geneProduct.fbc_name};
0874         else
0875             genes=getGenesFromGrRules(grRules);
0876         end
0877         model.genes=genes;
0878         model.grRules=grRules;
0879         [grRules,rxnGeneMat] = standardizeGrRules(model,true);
0880         model.grRules = grRules;
0881         model.rxnGeneMat = rxnGeneMat;
0882     end
0883 end
0884 
0885 if all(cellfun(@isempty,geneShortNames))
0886     if isfield(modelSBML,'fbc_geneProduct')
0887         for i=1:numel(genes)
0888             if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label)
0889                 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label;
0890             elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name)
0891                 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name;
0892             else
0893                 geneShortNames{i,1}='';
0894             end
0895         end
0896     end
0897 end
0898 
0899 %If any InChIs have been loaded
0900 if any(~cellfun(@isempty,metaboliteInChI))
0901     model.inchis=metaboliteInChI;
0902 end
0903 
0904 %If any formulas have been loaded
0905 if any(~cellfun(@isempty,metaboliteFormula))
0906     model.metFormulas=metaboliteFormula;
0907 end
0908 
0909 %If any charges have been loaded
0910 if ~isempty(metaboliteCharges)
0911     model.metCharges=metaboliteCharges;
0912 end
0913 
0914 %If any gene short names have been loaded
0915 if any(~cellfun(@isempty,geneShortNames))
0916     model.geneShortNames=geneShortNames;
0917 end
0918 
0919 %If any Miriam strings for compartments have been loaded
0920 if any(~cellfun(@isempty,compartmentMiriams))
0921     model.compMiriams=compartmentMiriams;
0922 end
0923 
0924 %If any Miriam strings for metabolites have been loaded
0925 if any(~cellfun(@isempty,metaboliteMiriams))
0926     model.metMiriams=metaboliteMiriams;
0927 end
0928 
0929 %If any subsystems have been loaded
0930 if any(~cellfun(@isempty,subsystems))
0931     model.subSystems=subsystems;
0932 end
0933 if any(rxnComps)
0934     if all(rxnComps)
0935         model.rxnComps=rxnComps;
0936     else
0937         if supressWarnings==false
0938             EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information';
0939             dispEM(EM,false,model.rxns(rxnComps==0));
0940         end
0941     end
0942 end
0943 
0944 %If any ec-codes have been loaded
0945 if any(~cellfun(@isempty,eccodes))
0946     model.eccodes=eccodes;
0947 end
0948 
0949 %If any Miriam strings for reactions have been loaded
0950 if any(~cellfun(@isempty,rxnMiriams))
0951     model.rxnMiriams=rxnMiriams;
0952 end
0953 
0954 %If any Miriam strings for genes have been loaded
0955 if any(~cellfun(@isempty,geneMiriams))
0956     model.geneMiriams=geneMiriams;
0957 end
0958 
0959 %If any protein strings have been loaded
0960 if any(~cellfun(@isempty,proteins))
0961     proteins = reshape(proteins,[],1);
0962     model.proteins=proteins;
0963 end
0964 
0965 model.unconstrained=metaboliteUnconstrained;
0966 
0967 %Convert SBML IDs back into their original strings. Here we are using part
0968 %from convertSBMLID, originating from the COBRA Toolbox
0969 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}');
0970 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}');
0971 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}');
0972 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}');
0973 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}');
0974 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}');
0975 
0976 if removePrefix
0977     [model, hasChanged]=removeIdentifierPrefix(model);
0978     dispEM(['The following fields have prefixes removed from all entries. '...
0979     'If this is undesired, run importModel with removePrefix as false. Example: '...
0980     'importModel(''filename.xml'',[],false);'],false,hasChanged)
0981 end
0982 
0983 %Remove unused fields
0984 if isempty(model.annotation)
0985     model=rmfield(model,'annotation');
0986 end
0987 if isempty(model.compOutside)
0988     model=rmfield(model,'compOutside');
0989 end
0990 if isempty(model.compMiriams)
0991     model=rmfield(model,'compMiriams');
0992 end
0993 if isempty(model.rxnComps)
0994     model=rmfield(model,'rxnComps');
0995 end
0996 if isempty(model.grRules)
0997     model=rmfield(model,'grRules');
0998 end
0999 if isempty(model.rxnGeneMat)
1000     model=rmfield(model,'rxnGeneMat');
1001 end
1002 if isempty(model.subSystems)
1003     model=rmfield(model,'subSystems');
1004 else
1005     model.subSystems(cellfun(@isempty,subsystems))={{''}};
1006     % If all subSystems have single entries, then unnest them
1007     if all(cellfun(@(x) iscell(x) && isscalar(x), model.subSystems))
1008         model.subSystems = transpose([model.subSystems{:}]);
1009     end
1010 end
1011 if isempty(model.eccodes)
1012     model=rmfield(model,'eccodes');
1013 end
1014 if isempty(model.rxnMiriams)
1015     model=rmfield(model,'rxnMiriams');
1016 end
1017 if cellfun(@isempty,model.rxnNotes)
1018     model=rmfield(model,'rxnNotes');
1019 end
1020 if cellfun(@isempty,model.rxnReferences)
1021     model=rmfield(model,'rxnReferences');
1022 end
1023 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores))
1024     model=rmfield(model,'rxnConfidenceScores');
1025 end
1026 if isempty(model.genes)
1027     model=rmfield(model,'genes');
1028 elseif isrow(model.genes)
1029     model.genes=transpose(model.genes);
1030 end
1031 if isempty(model.geneComps)
1032     model=rmfield(model,'geneComps');
1033 end
1034 if isempty(model.geneMiriams)
1035     model=rmfield(model,'geneMiriams');
1036 end
1037 if isempty(model.geneShortNames)
1038     model=rmfield(model,'geneShortNames');
1039 end
1040 if isempty(model.proteins)
1041     model=rmfield(model,'proteins');
1042 end
1043 if isempty(model.inchis)
1044     model=rmfield(model,'inchis');
1045 end
1046 if isempty(model.metFormulas)
1047     model=rmfield(model,'metFormulas');
1048 end
1049 if isempty(model.metMiriams)
1050     model=rmfield(model,'metMiriams');
1051 end
1052 if ~any(model.metCharges)
1053     model=rmfield(model,'metCharges');
1054 end
1055 
1056 %This just removes the grRules if no genes have been loaded
1057 if ~isfield(model,'genes') && isfield(model,'grRules')
1058     model=rmfield(model,'grRules');
1059 end
1060 
1061 %Print warnings about bad structure
1062 if supressWarnings==false
1063     checkModelStruct(model,false);
1064 end
1065 
1066 if removeExcMets==true
1067     model=simplifyModel(model);
1068 end
1069 end
1070 
1071 function fieldContent=parseNote(searchString,fieldName)
1072 %The function obtains the particular information from 'notes' field, using
1073 %fieldName as the dummy string
1074 
1075 fieldContent='';
1076 
1077 if strfind(searchString,fieldName)
1078     [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match');
1079     targetString=regexprep(targetString,'<p>|</p>','');
1080     targetString=regexprep(targetString,[fieldName, ':'],'');
1081     for i=1:numel(targetString)
1082         fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
1083     end
1084     fieldContent=regexprep(fieldContent,'^;|;$','');
1085 else
1086     fieldContent='';
1087 end
1088 end
1089 
1090 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName)
1091 
1092 fieldContent='';
1093 
1094 %Removing whitespace characters from the ending strings, which may occur in
1095 %several cases
1096 searchString=regexprep(searchString,'" />','"/>');
1097 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match');
1098 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
1099 targetString=regexprep(targetString,startString,'');
1100 targetString=regexprep(targetString,[fieldName midString],'');
1101 
1102 for i=1:numel(targetString)
1103     fieldContent=[fieldContent ';' strtrim(targetString{1,i})];
1104 end
1105 
1106 fieldContent=regexprep(fieldContent,'^;|;$','');
1107 end
1108 
1109 function miriamStruct=parseMiriam(searchString)
1110 %Generates miriam structure from annotation field
1111 
1112 %Finding whether miriams are written in the old or the new way
1113 if strfind(searchString,'urn:miriam:')
1114     startString='urn:miriam:';
1115 elseif strfind(searchString,'http://identifiers.org/')
1116     startString='http://identifiers.org/';
1117 elseif strfind(searchString,'https://identifiers.org/')
1118     startString='https://identifiers.org/';
1119 else
1120     miriamStruct=[];
1121     return;
1122 end
1123 
1124 miriamStruct=[];
1125 
1126 searchString=regexprep(searchString,'" />','"/>');
1127 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match');
1128 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>','');
1129 targetString=regexprep(targetString,startString,'');
1130 
1131 fwdslash  = contains(targetString,'/');
1132 midString = cell(numel(targetString),1);
1133 midString(fwdslash) = {'/'};
1134 midString(~fwdslash) = {':'};
1135 
1136 counter=0;
1137 for i=1:numel(targetString)
1138     if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once'))
1139         counter=counter+1;
1140         miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once');
1141         miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once');
1142         miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.','');
1143     end
1144 end
1145 end
1146 
1147 function miriam = addSBOtoMiriam(miriam,sboTerm)
1148 %Appends SBO term to miriam structure
1149 
1150 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]};  % convert to proper format
1151 if isempty(miriam)
1152     miriam.name = {'sbo'};
1153     miriam.value = sboTerm;
1154 elseif any(strcmp('sbo',miriam.name))
1155     currSbo = strcmp('sbo',miriam.name);
1156     miriam.value(currSbo) = sboTerm;
1157 else
1158     miriam.name(end+1) = {'sbo'};
1159     miriam.value(end+1) = sboTerm;
1160 end
1161 end