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.
   exportGeneComplexes true if gene complexes (all gene sets linked with
                       AND relationship) should be recognised and exported
                       (opt, default false)
   supressWarnings     true if warnings should be supressed (opt, default
                       false)
   sortIds             logical whether metabolites, reactions and genes
                       should be sorted alphabetically by their
                       identifiers (opt, default false)

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

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