exportModel

 exportModel
   Exports a constraint-based model to an SBML file (L3V1 FBCv2)

 Input:
   model               a model structure
   fileName            filename to export the model to. A dialog window
                       will open if no file name is specified.
   neverPrefix         true if prefixes are never added to identifiers,
                       even if start with e.g. digits. This might result
                       in invalid SBML files (optional, default false)
   supressWarnings     true if warnings should be supressed. This might
                       results in invalid SBML files, as no checks are
                       performed (optional, default false)
   sortIds             logical whether metabolites, reactions and genes
                       should be sorted alphabetically by their
                       identifiers (optional, default false)

 Usage: exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)

0001 function exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
0002 % exportModel
0003 %   Exports a constraint-based model to an SBML file (L3V1 FBCv2)
0004 %
0005 % Input:
0006 %   model               a model structure
0007 %   fileName            filename to export the model to. A dialog window
0008 %                       will open if no file name is specified.
0009 %   neverPrefix         true if prefixes are never added to identifiers,
0010 %                       even if start with e.g. digits. This might result
0011 %                       in invalid SBML files (optional, default false)
0012 %   supressWarnings     true if warnings should be supressed. This might
0013 %                       results in invalid SBML files, as no checks are
0014 %                       performed (optional, default false)
0015 %   sortIds             logical whether metabolites, reactions and genes
0016 %                       should be sorted alphabetically by their
0017 %                       identifiers (optional, default false)
0018 %
0019 % Usage: exportModel(model,fileName,neverPrefix,supressWarnings,sortIds)
0020 
0021 if nargin<2 || isempty(fileName)
0022     [fileName, pathName] = uiputfile({'*.xml;*.sbml'}, 'Select file for model export',[model.id '.xml']);
0023     if fileName == 0
0024         error('You should provide a file location')
0025     else
0026         fileName = fullfile(pathName,fileName);
0027     end
0028 end
0029 fileName=char(fileName);
0030 if nargin<3 || isempty(neverPrefix)
0031     neverPrefix=false;
0032 end
0033 if nargin<4 || isempty(supressWarnings)
0034     supressWarnings=false;
0035 end
0036 if nargin<5 || isempty(sortIds)
0037     sortIds=false;
0038 end
0039 if sortIds==true
0040     model=sortIdentifiers(model);
0041 end
0042 
0043 if isfield(model,'ec')
0044     warning("exportModel does not store information from the 'model.ec' structure. Use 'writeYAMLmodel(model)' to export all content from a GECKO model.")
0045 end
0046 
0047 %If no subSystems are defined, then no need to use groups package
0048 if isfield(model,'subSystems')
0049     modelHasSubsystems=true;
0050 else
0051     modelHasSubsystems=false;
0052 end
0053 
0054 %The default SBML format settings, which are used as input for appropriate
0055 %libSBML functions to generate the blank SBML model structure before using
0056 %exporting in with OutputSBML to xml file
0057 sbmlLevel=3;
0058 sbmlVersion=1;
0059 sbmlPackages={'fbc'};
0060 sbmlPackageVersions=2;
0061 if modelHasSubsystems
0062     sbmlPackages={sbmlPackages,'groups'};
0063     sbmlPackageVersions=[sbmlPackageVersions,1];
0064 end
0065 
0066 %Check if the "unconstrained" field is still present. This shows if
0067 %exchange metabolites have been removed
0068 if ~isfield(model,'unconstrained')
0069     model.unconstrained=zeros(numel(model.mets),1);
0070 end
0071 
0072 %If model id and name do not exist, make sure that default
0073 %strings are included
0074 if ~isfield(model,'id')
0075     fprintf('WARNING: The model is missing the "id" field. Uses "blankID". \n');
0076     model.id='blankID';
0077 end
0078 if ~isfield(model,'name')
0079     fprintf('WARNING: The model is missing the "name" field. Uses "blankName". \n');
0080     model.name='blankName';
0081 end
0082 
0083 % Add prefixes if required
0084 if ~neverPrefix
0085     [model,hasChanged] = addIdentifierPrefix(model);
0086     dispEM(['The following fields have one or more entries that do not start '...
0087         'with a letter or _ (conflicting with SBML specifications). Prefixes '...
0088         'are added to all entries in those fields:'],false,hasChanged)
0089 end
0090 
0091 %Check the model structure
0092 if supressWarnings==false
0093     checkModelStruct(model);
0094 end
0095 
0096 %Add several blank fields, if they do not exist already. This is to reduce
0097 %the number of conditions below
0098 if ~isfield(model,'compMiriams')
0099     model.compMiriams=cell(numel(model.comps),1);
0100 end
0101 if ~isfield(model,'inchis')
0102     model.inchis=cell(numel(model.mets),1);
0103 end
0104 if ~isfield(model,'metFormulas')
0105     model.metFormulas=cell(numel(model.mets),1);
0106 end
0107 if ~isfield(model,'metMiriams')
0108     model.metMiriams=cell(numel(model.mets),1);
0109 end
0110 if ~isfield(model,'geneMiriams') && isfield(model,'genes')
0111     model.geneMiriams=cell(numel(model.genes),1);
0112 end
0113 if ~isfield(model,'geneShortNames') && isfield(model,'genes')
0114     model.geneShortNames=cell(numel(model.genes),1);
0115 end
0116 if ~isfield(model,'proteins') && isfield(model,'genes')
0117     model.proteins=cell(numel(model.genes),1);
0118 end
0119 if ~isfield(model,'subSystems')
0120     model.subSystems=cell(numel(model.rxns),1);
0121 end
0122 if ~isfield(model,'eccodes')
0123     model.eccodes=cell(numel(model.rxns),1);
0124 end
0125 if ~isfield(model,'rxnReferences')
0126     model.rxnReferences=cell(numel(model.rxns),1);
0127 end
0128 if ~isfield(model,'rxnConfidenceScores')
0129     model.rxnConfidenceScores=NaN(numel(model.rxns),1);
0130 end
0131 if ~isfield(model,'rxnNotes')
0132     model.rxnNotes=cell(numel(model.rxns),1);
0133 end
0134 if ~isfield(model,'rxnMiriams')
0135     model.rxnMiriams=cell(numel(model.rxns),1);
0136 end
0137 
0138 if sbmlLevel<3
0139     %Check if genes have associated compartments
0140     if ~isfield(model,'geneComps') && isfield(model,'genes')
0141         if supressWarnings==false
0142             EM='There are no compartments specified for genes. All genes will be assigned to the first compartment. This is because the SBML structure requires all elements to be assigned to a compartment';
0143             dispEM(EM,false);
0144         end
0145         model.geneComps=ones(numel(model.genes),1);
0146     end
0147 end
0148 
0149 %Convert ids to SBML-convenient format. This is to avoid the data loss when
0150 %unsupported characters are included in ids. Here we are using part from
0151 %convertSBMLID, originating from the COBRA Toolbox
0152 model.rxns=regexprep(model.rxns,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0153 model.mets=regexprep(model.mets,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0154 model.comps=regexprep(model.comps,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0155 if isfield(model,'genes')
0156     problemGenes=find(~cellfun('isempty',regexp(model.genes,'([^0-9_a-zA-Z])')));
0157     originalGenes=model.genes(problemGenes);
0158     replacedGenes=regexprep(model.genes(problemGenes),'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0159     model.genes(problemGenes)=replacedGenes;
0160     for i=1:numel(problemGenes)
0161         model.grRules = regexprep(model.grRules, ['(^|\s|\()' originalGenes{i} '($|\s|\))'], ['$1' replacedGenes{i} '$2']);
0162     end
0163 end
0164 
0165 %Generate an empty SBML structure
0166 modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0167 modelSBML.metaid=model.id;
0168 modelSBML.id=regexprep(model.id,'([^0-9_a-zA-Z])','__${num2str($1+0)}__');
0169 modelSBML.name=model.name;
0170 
0171 if isfield(model,'annotation')
0172     if isfield(model.annotation,'note')
0173         modelSBML.notes=['<notes><body xmlns="http://www.w3.org/1999/xhtml"><p>',regexprep(model.annotation.note,'<p>|</p>',''),'</p></body></notes>'];
0174     end
0175 else
0176     modelSBML.notes='<notes><body xmlns="http://www.w3.org/1999/xhtml"><p>This file was generated using the exportModel function in RAVEN Toolbox 2 and OutputSBML in libSBML </p></body></notes>';
0177 end
0178 
0179 if isfield(model,'annotation')
0180     nameString='';
0181     if isfield(model.annotation,'familyName')
0182         if ~isempty(model.annotation.familyName)
0183             nameString=['<vCard:Family>' model.annotation.familyName '</vCard:Family>'];
0184         end
0185     end
0186     if isfield(model.annotation,'givenName')
0187         if ~isempty(model.annotation.givenName)
0188             nameString=[nameString '<vCard:Given>' model.annotation.givenName '</vCard:Given>'];
0189         end
0190     end
0191     email='';
0192     if isfield(model.annotation,'email')
0193         if ~isempty(model.annotation.email)
0194             email=['<vCard:EMAIL>' model.annotation.email '</vCard:EMAIL>'];
0195         end
0196     end
0197     org='';
0198     if isfield(model.annotation,'organization')
0199         if ~isempty(model.annotation.organization)
0200             org=['<vCard:ORG rdf:parseType="Resource"><vCard:Orgname>' model.annotation.organization '</vCard:Orgname></vCard:ORG>'];
0201         end
0202     end
0203     if ~isempty(nameString) || ~isempty(email) || ~isempty(org) % Only fill .annotation if ownership data is provided
0204         modelSBML.annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.id '">'];
0205         modelSBML.annotation=[modelSBML.annotation '<dc:creator><rdf:Bag><rdf:li rdf:parseType="Resource">'];
0206         if ~isempty(nameString)
0207             modelSBML.annotation=[modelSBML.annotation '<vCard:N rdf:parseType="Resource">' nameString '</vCard:N>'];
0208         end
0209         modelSBML.annotation=[modelSBML.annotation email org '</rdf:li></rdf:Bag></dc:creator>'];
0210         modelSBML.annotation=[modelSBML.annotation '<dcterms:created rdf:parseType="Resource">'...
0211             '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF></dcterms:created><dcterms:modified rdf:parseType="Resource">'...
0212             '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF></dcterms:modified>'];
0213         if isfield(model.annotation,'taxonomy')
0214             modelSBML.annotation=[modelSBML.annotation '<bqbiol:is><rdf:Bag><rdf:li rdf:resource="https://identifiers.org/taxonomy/' regexprep(model.annotation.taxonomy,'taxonomy/','') '"/></rdf:Bag></bqbiol:is>'];
0215         end
0216         modelSBML.annotation=[modelSBML.annotation '</rdf:Description></rdf:RDF></annotation>'];
0217     end
0218 end
0219 
0220 %Prepare compartments
0221 for i=1:numel(model.comps)
0222     %Add the default values, as these will be the same in all entries
0223     if i==1
0224         if isfield(modelSBML.compartment, 'sboTerm')
0225             modelSBML.compartment(i).sboTerm=290;
0226         end
0227         if isfield(modelSBML.compartment, 'spatialDimensions')
0228             modelSBML.compartment(i).spatialDimensions=3;
0229         end
0230         if isfield(modelSBML.compartment, 'size')
0231             modelSBML.compartment(i).size=1;
0232         end
0233         if isfield(modelSBML.compartment, 'constant')
0234             modelSBML.compartment(i).constant=1;
0235         end
0236         if isfield(modelSBML.compartment, 'isSetSize')
0237             modelSBML.compartment(i).isSetSize=1;
0238         end
0239         if isfield(modelSBML.compartment, 'isSetSpatialDimensions')
0240             modelSBML.compartment(i).isSetSpatialDimensions=1;
0241         end
0242     end
0243     %Copy the default values to the next entry as long as it is not the
0244     %last one
0245     if i<numel(model.comps)
0246         modelSBML.compartment(i+1)=modelSBML.compartment(i);
0247     end
0248     
0249     if isfield(modelSBML.compartment,'metaid')
0250         modelSBML.compartment(i).metaid=model.comps{i};
0251     end
0252     %Prepare Miriam strings
0253     if ~isempty(model.compMiriams{i})
0254         [~,sbo_ind] = ismember('sbo',model.compMiriams{i}.name);
0255         if sbo_ind > 0
0256             modelSBML.compartment(i).sboTerm=str2double(regexprep(model.compMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0257             % remove the SBO term from compMiriams so the information is
0258             % not duplicated in the "annotation" field later on
0259             model.compMiriams{i}.name(sbo_ind) = [];
0260             model.compMiriams{i}.value(sbo_ind) = [];
0261         end
0262     end
0263     if ~isempty(model.compMiriams{i}) && isfield(modelSBML.compartment(i),'annotation')
0264         modelSBML.compartment(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.comps{i} '">'];
0265         modelSBML.compartment(i).annotation=[modelSBML.compartment(i).annotation '<bqbiol:is><rdf:Bag>'];
0266         modelSBML.compartment(i).annotation=[modelSBML.compartment(i).annotation getMiriam(model.compMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0267     end
0268     if isfield(modelSBML.compartment, 'name')
0269         modelSBML.compartment(i).name=model.compNames{i};
0270     end
0271     if isfield(modelSBML.compartment, 'id')
0272         modelSBML.compartment(i).id=model.comps{i};
0273     end
0274     
0275 end
0276 
0277 %Begin writing species
0278 for i=1:numel(model.mets)
0279     %Add the default values, as these will be the same in all entries
0280     if i==1
0281         if isfield(modelSBML.species, 'sboTerm')
0282             modelSBML.species(i).sboTerm=247;
0283         end
0284         if isfield(modelSBML.species, 'initialAmount')
0285             modelSBML.species(i).initialAmount=1;
0286         end
0287         if isfield(modelSBML.species, 'initialConcentration')
0288             modelSBML.species(i).initialConcentration=0;
0289         end
0290         if isfield(modelSBML.species, 'isSetInitialAmount')
0291             modelSBML.species(i).isSetInitialAmount=1;
0292         end
0293         if isfield(modelSBML.species, 'isSetInitialConcentration')
0294             modelSBML.species(i).isSetInitialConcentration=1;
0295         end
0296     end
0297     %Copy the default values to the next entry as long as it is not the
0298     %last one
0299     if i<numel(model.mets)
0300         modelSBML.species(i+1)=modelSBML.species(i);
0301     end
0302     
0303     if isfield(modelSBML.species,'metaid')
0304         modelSBML.species(i).metaid=model.mets{i};
0305     end
0306     if isfield(modelSBML.species, 'name')
0307         modelSBML.species(i).name=model.metNames{i};
0308     end
0309     if isfield(modelSBML.species, 'id')
0310         modelSBML.species(i).id=model.mets{i};
0311     end
0312     if isfield(modelSBML.species, 'compartment')
0313         modelSBML.species(i).compartment=model.comps{model.metComps(i)};
0314     end
0315     if isfield(model,'unconstrained')
0316         if model.unconstrained(i)
0317             modelSBML.species(i).boundaryCondition=1;
0318         end
0319     end
0320     if isfield(modelSBML.species, 'fbc_charge') && isfield(model,'metCharges')
0321         if ~isnan(model.metCharges(i))
0322             modelSBML.species(i).fbc_charge=model.metCharges(i);
0323             modelSBML.species(i).isSetfbc_charge=1;
0324         else
0325             modelSBML.species(i).isSetfbc_charge=0;
0326         end
0327     end
0328     if ~isempty(model.metMiriams{i})
0329         [~,sbo_ind] = ismember('sbo',model.metMiriams{i}.name);
0330         if sbo_ind > 0
0331             modelSBML.species(i).sboTerm=str2double(regexprep(model.metMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0332             % remove the SBO term from metMiriams so the information is
0333             % not duplicated in the "annotation" field later on
0334             model.metMiriams{i}.name(sbo_ind) = [];
0335             model.metMiriams{i}.value(sbo_ind) = [];
0336         end
0337     end
0338     if isfield(modelSBML.species,'annotation')
0339         if ~isempty(model.metMiriams{i}) || ~isempty(model.metFormulas{i})
0340             hasInchi=false;
0341             if ~isempty(model.metFormulas{i})
0342                 %Only export formula if there is no InChI. This is because
0343                 %the metFormulas field is populated by InChIs if available
0344                 if ~isempty(model.inchis{i})
0345                     hasInchi=true;
0346                 end
0347                 if hasInchi==false
0348                     modelSBML.species(i).fbc_chemicalFormula=model.metFormulas{i};
0349                 end
0350             end
0351             if ~isempty(model.metMiriams{i}) || hasInchi==true
0352                 modelSBML.species(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.mets{i} '">'];
0353                 modelSBML.species(i).annotation=[modelSBML.species(i).annotation '<bqbiol:is><rdf:Bag>'];
0354                 if ~isempty(model.metMiriams{i})
0355                     modelSBML.species(i).annotation=[modelSBML.species(i).annotation getMiriam(model.metMiriams{i})];
0356                 end
0357                 if hasInchi==true
0358                     modelSBML.species(i).annotation=[modelSBML.species(i).annotation '<rdf:li rdf:resource="https://identifiers.org/inchi/InChI=' regexprep(model.inchis{i},'^InChI=','') '"/>'];
0359                     modelSBML.species(i).fbc_chemicalFormula=char(regexp(model.inchis{i}, '/(\w+)/', 'tokens', 'once'));
0360                 end
0361                 modelSBML.species(i).annotation=[modelSBML.species(i).annotation '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0362             end
0363         end
0364     end
0365 end
0366 
0367 if isfield(model,'genes')
0368     for i=1:numel(model.genes)
0369         %Add the default values, as these will be the same in all entries
0370         if i==1
0371             if isfield(modelSBML.fbc_geneProduct, 'sboTerm')
0372                 modelSBML.fbc_geneProduct(i).sboTerm=243;
0373             end
0374         end
0375         %Copy the default values to the next index as long as it is not the
0376         %last one
0377         if i<numel(model.genes)
0378             modelSBML.fbc_geneProduct(i+1)=modelSBML.fbc_geneProduct(i);
0379         end
0380         
0381         if isfield(modelSBML.fbc_geneProduct,'metaid')
0382             modelSBML.fbc_geneProduct(i).metaid=model.genes{i};
0383         end
0384         if ~isempty(model.geneMiriams{i})
0385             [~,sbo_ind] = ismember('sbo',model.geneMiriams{i}.name);
0386             if sbo_ind > 0
0387                 modelSBML.fbc_geneProduct(i).sboTerm=str2double(regexprep(model.geneMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0388                 % remove the SBO term from compMiriams so the information is
0389                 % not duplicated in the "annotation" field later on
0390                 model.geneMiriams{i}.name(sbo_ind) = [];
0391                 model.geneMiriams{i}.value(sbo_ind) = [];
0392             end
0393         end
0394         if ~isempty(model.geneMiriams{i}) && isfield(modelSBML.fbc_geneProduct(i),'annotation')
0395             modelSBML.fbc_geneProduct(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.genes{i} '">'];
0396             modelSBML.fbc_geneProduct(i).annotation=[modelSBML.fbc_geneProduct(i).annotation '<bqbiol:is><rdf:Bag>'];
0397             modelSBML.fbc_geneProduct(i).annotation=[modelSBML.fbc_geneProduct(i).annotation getMiriam(model.geneMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0398         end
0399         if isfield(modelSBML.fbc_geneProduct, 'fbc_id')
0400             modelSBML.fbc_geneProduct(i).fbc_id=model.genes{i};
0401         end
0402         if isfield(modelSBML.fbc_geneProduct, 'fbc_label') && isfield(model,'geneShortNames')
0403             if isempty(model.geneShortNames{i})
0404                 modelSBML.fbc_geneProduct(i).fbc_label=model.genes{i};
0405             else
0406                 modelSBML.fbc_geneProduct(i).fbc_label=model.geneShortNames{i};
0407             end
0408         end
0409         if isfield(modelSBML.fbc_geneProduct, 'fbc_name') && isfield(model,'proteins')
0410             if ~isempty(model.proteins{i})
0411                 modelSBML.fbc_geneProduct(i).fbc_name=model.proteins{i};
0412             end
0413         end
0414     end
0415 end
0416 
0417 %Generate a list of unique fbc_bound names
0418 totalValues=[model.lb; model.ub];
0419 totalNames=cell(size(totalValues,1),1);
0420 
0421 listUniqueValues=unique(totalValues);
0422 
0423 for i=1:length(listUniqueValues)
0424     listUniqueNames{i,1}=['FB',num2str(i),'N',num2str(abs(round(listUniqueValues(i))))]; % create unique flux bound IDs.
0425     ind=find(ismember(totalValues,listUniqueValues(i)));
0426     totalNames(ind)=listUniqueNames(i,1);
0427 end
0428 
0429 for i=1:length(listUniqueNames)
0430     %Add the default values, as these will be the same in all entries
0431     if i==1
0432         if isfield(modelSBML.parameter, 'constant')
0433             modelSBML.parameter(i).constant=1;
0434         end
0435         if isfield(modelSBML.parameter, 'isSetValue')
0436             modelSBML.parameter(i).isSetValue=1;
0437         end
0438     end
0439     %Copy the default values to the next index as long as it is not the
0440     %last one
0441     if i<numel(listUniqueNames)
0442         modelSBML.parameter(i+1)=modelSBML.parameter(i);
0443     end
0444     modelSBML.parameter(i).id=listUniqueNames{i};
0445     modelSBML.parameter(i).value=listUniqueValues(i);
0446 end
0447 
0448 for i=1:numel(model.rxns)
0449     %Add the default values, as these will be the same in all entries
0450     if i==1
0451         if isfield(modelSBML.reaction, 'sboTerm')
0452             modelSBML.reaction(i).sboTerm=176;
0453         end
0454         if isfield(modelSBML.reaction, 'isSetFast')
0455             modelSBML.reaction(i).isSetFast=1;
0456         end
0457     end
0458     %Copy the default values to the next index as long as it is not the
0459     %last one
0460     if i<numel(model.rxns)
0461         modelSBML.reaction(i+1)=modelSBML.reaction(i);
0462     end
0463     
0464     if isfield(modelSBML.reaction,'metaid')
0465         modelSBML.reaction(i).metaid=model.rxns{i};
0466     end
0467     
0468     %Export notes information
0469     if (~isnan(model.rxnConfidenceScores(i)) || ~isempty(model.rxnReferences{i}) || ~isempty(model.rxnNotes{i}))
0470         modelSBML.reaction(i).notes='<notes><body xmlns="http://www.w3.org/1999/xhtml">';
0471         if ~isnan(model.rxnConfidenceScores(i))
0472             modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '<p>Confidence Level: ' num2str(model.rxnConfidenceScores(i)) '</p>'];
0473         end
0474         if ~isempty(model.rxnReferences{i})
0475             modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '<p>AUTHORS: ' model.rxnReferences{i} '</p>'];
0476         end
0477         if ~isempty(model.rxnNotes{i})
0478             modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '<p>NOTES: ' model.rxnNotes{i} '</p>'];
0479         end
0480         modelSBML.reaction(i).notes=[modelSBML.reaction(i).notes '</body></notes>'];
0481     end
0482     
0483     % Export SBO terms from rxnMiriams
0484     if ~isempty(model.rxnMiriams{i})
0485         [~,sbo_ind] = ismember('sbo',model.rxnMiriams{i}.name);
0486         if sbo_ind > 0
0487             modelSBML.reaction(i).sboTerm=str2double(regexprep(model.rxnMiriams{i}.value{sbo_ind},'SBO:','','ignorecase'));
0488             % remove the SBO term from rxnMiriams so the information is not
0489             % duplicated in the "annotation" field later on
0490             model.rxnMiriams{i}.name(sbo_ind) = [];
0491             model.rxnMiriams{i}.value(sbo_ind) = [];
0492         end
0493     end
0494     
0495     %Export annotation information from rxnMiriams
0496     if (~isempty(model.rxnMiriams{i}) && isfield(modelSBML.reaction(i),'annotation')) || ~isempty(model.eccodes{i})
0497         modelSBML.reaction(i).annotation=['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#meta_' model.rxns{i} '">'];
0498         modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation '<bqbiol:is><rdf:Bag>'];
0499         if ~isempty(model.eccodes{i})
0500             eccodes=regexp(model.eccodes{i},';','split');
0501             for j=1:numel(eccodes)
0502                 modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation  '<rdf:li rdf:resource="https://identifiers.org/ec-code/' regexprep(eccodes{j},'ec-code/|EC','') '"/>'];
0503             end
0504         end
0505         modelSBML.reaction(i).annotation=[modelSBML.reaction(i).annotation getMiriam(model.rxnMiriams{i}) '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0506     end
0507     
0508     if isfield(modelSBML.reaction, 'name')
0509         modelSBML.reaction(i).name=model.rxnNames{i};
0510     end
0511     if isfield(modelSBML.reaction, 'id')
0512         modelSBML.reaction(i).id=model.rxns{i};
0513     end
0514     
0515     %Add the information about reactants and products
0516     involvedMets=addReactantsProducts(model,modelSBML,i);
0517     for j=1:numel(involvedMets.reactant)
0518         if j<numel(involvedMets.reactant)
0519             modelSBML.reaction(i).reactant(j+1)=modelSBML.reaction(i).reactant(j);
0520         end
0521         modelSBML.reaction(i).reactant(j).species=involvedMets.reactant(j).species;
0522         modelSBML.reaction(i).reactant(j).stoichiometry=involvedMets.reactant(j).stoichiometry;
0523         modelSBML.reaction(i).reactant(j).isSetStoichiometry=involvedMets.reactant(j).isSetStoichiometry;
0524         modelSBML.reaction(i).reactant(j).constant=involvedMets.reactant(j).constant;
0525     end
0526     if numel(involvedMets.reactant)==0
0527         modelSBML.reaction(i).reactant='';
0528     end
0529     for j=1:numel(involvedMets.product)
0530         if j<numel(involvedMets.product)
0531             modelSBML.reaction(i).product(j+1)=modelSBML.reaction(i).product(j);
0532         end
0533         modelSBML.reaction(i).product(j).species=involvedMets.product(j).species;
0534         modelSBML.reaction(i).product(j).stoichiometry=involvedMets.product(j).stoichiometry;
0535         modelSBML.reaction(i).product(j).isSetStoichiometry=involvedMets.product(j).isSetStoichiometry;
0536         modelSBML.reaction(i).product(j).constant=involvedMets.product(j).constant;
0537     end
0538     if numel(involvedMets.product)==0
0539         modelSBML.reaction(i).product='';
0540     end
0541     %Export reversibility information. Reactions are irreversible by
0542     %default
0543     if model.rev(i)==1
0544         modelSBML.reaction(i).reversible=1;
0545     end
0546     if isfield(model, 'rxnComps')
0547         modelSBML.reaction(i).compartment=model.comps{model.rxnComps(i)};
0548     end
0549     if isfield(model, 'grRules')
0550         modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association=model.grRules{i};
0551     end
0552     modelSBML.reaction(i).fbc_lowerFluxBound=totalNames{i};
0553     modelSBML.reaction(i).fbc_upperFluxBound=totalNames{length(model.lb)+i};
0554 end
0555 
0556 %Prepare subSystems Code taken from COBRA functions getModelSubSystems,
0557 %writeSBML, findRxnsFromSubSystem under GNU General Public License v3.0,
0558 %license file in readme/GPL.MD. Code modified for RAVEN
0559 if modelHasSubsystems
0560     modelSBML.groups_group.groups_kind = 'partonomy';
0561     modelSBML.groups_group.sboTerm = 633;
0562     tmpStruct=modelSBML.groups_group;
0563 
0564     rxns=model.rxns;
0565     if ~any(cellfun(@iscell,model.subSystems))
0566         if ~any(~cellfun(@isempty,model.subSystems))
0567             subSystems = {};
0568         else
0569             subSystems = setdiff(model.subSystems,'');
0570         end
0571     else
0572         orderedSubs = cellfun(@(x) columnVector(x),model.subSystems,'UniformOUtput',false);
0573         subSystems = setdiff(vertcat(orderedSubs{:}),'');
0574     end
0575     if isempty(subSystems)
0576         subSystems = {};
0577     end
0578     if ~isempty(subSystems)
0579         %Build the groups for the group package
0580         groupIDs = strcat('group',cellfun(@num2str, num2cell(1:length(subSystems)),'UniformOutput',false));
0581         for i = 1:length(subSystems)
0582             cgroup = tmpStruct;
0583             if ~any(cellfun(@iscell,model.subSystems))
0584                 present = ismember(model.subSystems,subSystems{i});
0585             else
0586                 present = cellfun(@(x) any(ismember(x,subSystems{i})),model.subSystems);
0587             end
0588             groupMembers = rxns(present);
0589             for j = 1:numel(groupMembers)
0590                 cMember = tmpStruct.groups_member;
0591                 cMember.groups_idRef = groupMembers{j};
0592                 if j == 1
0593                     cgroup.groups_member = cMember;
0594                 else
0595                     cgroup.groups_member(j) = cMember;
0596                 end
0597             end
0598             cgroup.groups_id = groupIDs{i};
0599             cgroup.groups_name = subSystems{i};
0600             if i == 1
0601                 modelSBML.groups_group = cgroup;
0602             else
0603                 modelSBML.groups_group(i) = cgroup;
0604             end
0605         end
0606     end
0607 end
0608 
0609 %Prepare fbc_objective subfield
0610 
0611 modelSBML.fbc_objective.fbc_type='maximize';
0612 modelSBML.fbc_objective.fbc_id='obj';
0613 
0614 ind=find(model.c);
0615 
0616 if isempty(ind)
0617     modelSBML.fbc_objective.fbc_fluxObjective.fbc_coefficient=0;
0618 else
0619     for i=1:length(ind)
0620         %Copy the default values to the next index as long as it is not the
0621         %last one
0622         if i<numel(ind)
0623             modelSBML.reaction(i+1)=modelSBML.reaction(i);
0624         end
0625         values=model.c(model.c~=0);
0626         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction=modelSBML.reaction(ind(i)).id;
0627         modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient=values(i);
0628         modelSBML.fbc_objective(i).fbc_fluxObjective.isSetfbc_coefficient=1;
0629     end
0630 end
0631 
0632 modelSBML.fbc_activeObjective=modelSBML.fbc_objective.fbc_id;
0633 
0634 fbcStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/fbc/version',num2str(sbmlPackageVersions(1))];
0635 if modelHasSubsystems
0636     groupStr=['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/groups/version',num2str(sbmlPackageVersions(2))];
0637     modelSBML.namespaces=struct('prefix',{'','fbc','groups'},...
0638     'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
0639     fbcStr,groupStr});
0640 else
0641     modelSBML.namespaces=struct('prefix',{'','fbc'},...
0642     'uri',{['http://www.sbml.org/sbml/level', num2str(sbmlLevel), '/version', num2str(sbmlVersion), '/core'],...
0643     fbcStr});
0644 end
0645 
0646 if sbmlPackageVersions(1) == 2
0647     modelSBML.fbc_strict=1;
0648     modelSBML.isSetfbc_strict = 1;
0649 end
0650 
0651 modelSBML.rule=[];
0652 modelSBML.constraint=[];
0653 
0654 [ravenDir,prevDir]=findRAVENroot();
0655 fileName=checkFileExistence(fileName,1,true,false);
0656 
0657 OutputSBML_RAVEN(modelSBML,fileName,1,0,[1,0]);
0658 end
0659 
0660 
0661 function modelSBML=getSBMLStructure(sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions)
0662 %Returns the blank SBML model structure by using appropriate libSBML
0663 %functions. This creates structure by considering three levels
0664 
0665 sbmlFieldNames=getStructureFieldnames('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0666 sbmlDefaultValues=getDefaultValues('model',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0667 
0668 for i=1:numel(sbmlFieldNames)
0669     modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
0670     sbmlSubfieldNames=getStructureFieldnames(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0671     sbmlSubfieldValues=getDefaultValues(sbmlFieldNames{1,i},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0672     if ~strcmp(sbmlFieldNames{1,i},'event') && ~strcmp(sbmlFieldNames{1,i},'functionDefinition') && ~strcmp(sbmlFieldNames{1,i},'initialAssignment')
0673         for j=1:numel(sbmlSubfieldNames)
0674             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j})=sbmlSubfieldValues{1,j};
0675             sbmlSubsubfieldNames=getStructureFieldnames(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0676             sbmlSubsubfieldValues=getDefaultValues(sbmlSubfieldNames{1,j},sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0677             if ~strcmp(sbmlSubfieldNames{1,j},'modifier') && ~strcmp(sbmlSubfieldNames{1,j},'kineticLaw')
0678                 for k=1:numel(sbmlSubsubfieldNames)
0679                     %'compartment' and 'species' fields are not supposed to
0680                     %have their standalone structures if they are subfields
0681                     %or subsubfields
0682                     if ~strcmp(sbmlSubfieldNames{1,j},'compartment') && ~strcmp(sbmlSubfieldNames{1,j},'species')
0683                         modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k})=sbmlSubsubfieldValues{1,k};
0684                     end
0685                     %If it is fbc_association in the third level, we need
0686                     %to establish the fourth level, since libSBML requires
0687                     %it
0688                     if strcmp(sbmlSubsubfieldNames{1,k},'fbc_association')
0689                         fbc_associationFieldNames=getStructureFieldnames('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0690                         fbc_associationFieldValues=getDefaultValues('fbc_association',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0691                         for l=1:numel(fbc_associationFieldNames)
0692                             modelSBML.(sbmlFieldNames{1,i}).(sbmlSubfieldNames{1,j}).(sbmlSubsubfieldNames{1,k}).(fbc_associationFieldNames{1,l})=fbc_associationFieldValues{1,l};
0693                         end
0694                     end
0695                 end
0696             end
0697         end
0698     end
0699     if ~isstruct(modelSBML.(sbmlFieldNames{1,i}))
0700         modelSBML.(sbmlFieldNames{1,i})=sbmlDefaultValues{1,i};
0701     end
0702 end
0703 
0704 modelSBML.unitDefinition.id='mmol_per_gDW_per_hr';
0705 
0706 unitFieldNames=getStructureFieldnames('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0707 unitDefaultValues=getDefaultValues('unit',sbmlLevel,sbmlVersion,sbmlPackages,sbmlPackageVersions);
0708 
0709 kinds={'mole','gram','second'};
0710 exponents=[1 -1 -1];
0711 scales=[-3 0 0];
0712 multipliers=[1 1 1*60*60];
0713 
0714 for i=1:numel(unitFieldNames)
0715     modelSBML.unitDefinition.unit(1).(unitFieldNames{1,i})=unitDefaultValues{1,i};
0716     for j=1:3
0717         modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=unitDefaultValues{1,i};
0718         if strcmp(unitFieldNames{1,i},'kind')
0719             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=kinds{j};
0720         elseif strcmp(unitFieldNames{1,i},'exponent')
0721             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=exponents(j);
0722         elseif strcmp(unitFieldNames{1,i},'scale')
0723             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=scales(j);
0724         elseif strcmp(unitFieldNames{1,i},'multiplier')
0725             modelSBML.unitDefinition.unit(j).(unitFieldNames{1,i})=multipliers(j);
0726         end
0727     end
0728 end
0729 end
0730 
0731 function miriamString=getMiriam(miriamStruct)
0732 %Returns a string with list elements for a miriam structure ('<rdf:li
0733 %rdf:resource="https://identifiers.org/go/GO:0005739"/>' for example). This
0734 %is just to speed up things since this is done many times during the
0735 %exporting
0736 
0737 miriamString='';
0738 if isfield(miriamStruct,'name')
0739     for i=1:numel(miriamStruct.name)
0740         miriamString=[miriamString '<rdf:li rdf:resource="https://identifiers.org/' miriamStruct.name{i} '/' miriamStruct.value{i} '"/>'];
0741     end
0742 end
0743 end
0744 
0745 function [tmp_Rxn]=addReactantsProducts(model,sbmlModel,i)
0746 %This function provides reactants and products for particular reaction. The
0747 %function was 'borrowed' from writeSBML in COBRA toolbox, lines 663-679
0748 
0749 met_idx = find(model.S(:, i));
0750 tmp_Rxn.product=[];
0751 tmp_Rxn.reactant=[];
0752 for j_met=1:size(met_idx,1)
0753     tmp_idx = met_idx(j_met,1);
0754     sbml_tmp_species_ref.species = sbmlModel.species(tmp_idx).id;
0755     met_stoich = model.S(tmp_idx, i);
0756     sbml_tmp_species_ref.stoichiometry = abs(met_stoich);
0757     sbml_tmp_species_ref.isSetStoichiometry=1;
0758     sbml_tmp_species_ref.constant=1;
0759     if (met_stoich > 0)
0760         tmp_Rxn.product = [ tmp_Rxn.product, sbml_tmp_species_ref ];
0761     else
0762         tmp_Rxn.reactant = [ tmp_Rxn.reactant, sbml_tmp_species_ref];
0763     end
0764 end
0765 end
0766 
0767 function vecT = columnVector(vec)
0768 % Code below taken from COBRA Toolbox under GNU General Public License v3.0
0769 % license file in readme/GPL.MD.
0770 %
0771 % Converts a vector to a column vector
0772 %
0773 % USAGE:
0774 %
0775 %   vecT = columnVector(vec)
0776 %
0777 % INPUT:
0778 %   vec:     a vector
0779 %
0780 % OUTPUT:
0781 %   vecT:    a column vector
0782 
0783 [n, m] = size(vec);
0784 
0785 if n < m
0786     vecT = vec';
0787 else
0788     vecT = vec;
0789 end
0790 end