Home > RAVEN > exportModel.m

exportModel

PURPOSE ^

exportModel

SYNOPSIS ^

function exportModel(model,fileName)

DESCRIPTION ^

 exportModel
   Exports a constraint-based model to a SBML file

   model         a model structure
   fileName      the filename  to export the model to

   The resulting SBML file is in the format defined in the yeast consensus
   model. No checks at all are made to ensure the correctness of the
   model.
   
   NOTE: This function is currently not able to deal with gene complexes
         and saves them only as 'or' relationships. Use
         exportToExcelFormat if you have gene complexes in your model.

   Usage: exportModel(model,fileName)

   Rasmus Agren, 2013-02-06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function exportModel(model,fileName)
0002 % exportModel
0003 %   Exports a constraint-based model to a SBML file
0004 %
0005 %   model         a model structure
0006 %   fileName      the filename  to export the model to
0007 %
0008 %   The resulting SBML file is in the format defined in the yeast consensus
0009 %   model. No checks at all are made to ensure the correctness of the
0010 %   model.
0011 %
0012 %   NOTE: This function is currently not able to deal with gene complexes
0013 %         and saves them only as 'or' relationships. Use
0014 %         exportToExcelFormat if you have gene complexes in your model.
0015 %
0016 %   Usage: exportModel(model,fileName)
0017 %
0018 %   Rasmus Agren, 2013-02-06
0019 %
0020 
0021 %Check if the "unconstrained" field is still present. This shows if
0022 %exchange metabolites have been removed
0023 if ~isfield(model,'unconstrained')
0024    fprintf('WARNING: There is no unconstrained field in the model structure. This means that no metabolites are considered exchange metabolites\n');
0025    model.unconstrained=zeros(numel(model.mets),1);
0026 end
0027 
0028 %Check if genes have associated compartments. This is only needed for the
0029 %new format and has no real impact on anything.
0030 if ~isfield(model,'geneComps') && isfield(model,'genes')
0031     fprintf('WARNING: There are no compartments specified for the genes. All genes will be assigned to a new compartment named "Unspecified" and abbreviated "u".\n');
0032     model.comps=[model.comps;{'u'}];
0033     model.compNames=[model.compNames;{'Unspecified'}];
0034 
0035     if isfield(model,'compOutside')
0036         model.compOutside=[model.compOutside;{''}];
0037     end
0038     if isfield(model,'compMiriams')
0039         model.compMiriams=[model.compMiriams;{[]}];
0040     end
0041     geneComps=cell(numel(model.genes),1);
0042     geneComps(:,:)=cellstr('u');
0043     model.geneComps=geneComps;
0044 end
0045 
0046 %Generate temporary filename
0047 tempFile=tempname;
0048 
0049 %Open a stream
0050 fid = fopen(tempFile,'w');
0051 
0052 %Writes the intro
0053 intro=['<?xml version="1.0" encoding="UTF-8" ?>'...
0054 '<sbml xmlns="http://www.sbml.org/sbml/level2/version3" level="2" version="3">'...
0055 '<model metaid="metaid_' model.id '" id="' model.id '" name="' model.description '">'...
0056 '<notes>'...
0057 '<body xmlns="http://www.w3.org/1999/xhtml">'...
0058 '<p>This file was generated using the exportModel function in RAVEN Toolbox.</p>'...
0059 '</body>'...
0060 '</notes>'...
0061 '<annotation>'...
0062 '<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/">'...
0063 '<rdf:Description rdf:about="#metaid_' model.id '">'...
0064 '<dc:creator rdf:parseType="Resource">'...
0065 '<rdf:Bag>'...
0066 '<rdf:li rdf:parseType="Resource">'...
0067 '<vCard:N rdf:parseType="Resource">'...
0068 '<vCard:Family>Agren</vCard:Family>' ...
0069 '<vCard:Given>Rasmus</vCard:Given>' ...
0070 '</vCard:N>'...
0071 '<vCard:EMAIL>rasmus.agren@chalmers.se</vCard:EMAIL>'... 
0072 '<vCard:ORG>'...
0073 '<vCard:Orgname>Chalmers University of Technology, Gothenburg</vCard:Orgname>'... 
0074 '</vCard:ORG>'...
0075 '</rdf:li>'...
0076 '</rdf:Bag>'...
0077 '</dc:creator>'...
0078 '<dcterms:created rdf:parseType="Resource">'...
0079 '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF>'... 
0080 '</dcterms:created>'...
0081 '<dcterms:modified rdf:parseType="Resource">'...
0082 '<dcterms:W3CDTF>' datestr(now,'yyyy-mm-ddTHH:MM:SSZ') '</dcterms:W3CDTF>'... 
0083 '</dcterms:modified>'...
0084 '</rdf:Description>'...
0085 '</rdf:RDF>'...
0086 '</annotation>'...
0087 '<listOfUnitDefinitions>'...
0088 '<unitDefinition id="mmol_per_gDW_per_hr">'...
0089 '<listOfUnits>'...
0090 '<unit kind="mole" scale="-3"/>'...
0091 '<unit kind="second" multiplier="0.00027778" exponent="-1"/>'...
0092 '</listOfUnits>'...
0093 '</unitDefinition>'...
0094 '</listOfUnitDefinitions>'...
0095 '<listOfCompartments>'];
0096 
0097 %Write intro
0098 fprintf(fid,intro);
0099 
0100 %Write compartments
0101 for i=1:numel(model.comps)
0102     %Check if it's outside anything
0103     if isfield(model, 'compOutside')
0104         if ~isempty(model.compOutside{i})
0105             append=[' outside="C_' model.compOutside{i} '" spatialDimensions="3"'];
0106         else
0107             append=' spatialDimensions="3"';
0108         end
0109     else
0110         append=' spatialDimensions="3"';
0111     end
0112     
0113     fprintf(fid,['<compartment metaid="metaid_C_' model.comps{i} '" id="C_' model.comps{i} '" name="' model.compNames{i} '"' append ' sboTerm="SBO:0000290">']);
0114 
0115     %Print associated Miriam strings
0116     if isfield(model,'compMiriams')
0117         miriamString=getMiriam(model.compMiriams{i});
0118     else
0119         miriamString=[];
0120     end
0121 
0122     if ~isempty(miriamString)
0123         compinfo=['<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='...
0124                 '"http://purl.org/dc/terms/" xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" '...
0125                 'xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#metaid_C_' model.comps{i} '">'...
0126                 '<bqbiol:is><rdf:Bag>' miriamString '</rdf:Bag></bqbiol:is>'...
0127                 '</rdf:Description></rdf:RDF></annotation></compartment>'];
0128     else
0129         compinfo='</compartment>\n';
0130     end
0131     fprintf(fid,compinfo);
0132 end
0133 
0134 intro='</listOfCompartments><listOfSpecies>';
0135 fprintf(fid,intro);
0136 
0137 %Begin writing species
0138 for i=1:numel(model.mets)
0139     if model.unconstrained(i)
0140         unbounded='true';
0141     else
0142         unbounded='false';
0143     end
0144     
0145     toprint=['<species metaid="metaid_M_' model.mets{i} '" id="M_' model.mets{i} '" name="' model.metNames{i} '" compartment="C_' model.comps{model.metComps(i)} '" boundaryCondition="' unbounded '" sboTerm="SBO:0000299">'];
0146     
0147     %Print some stuff if there is a formula for the compound
0148     if isfield(model,'metFormulas')
0149         %Only print formula if there is no InChI. The formula is loaded
0150         %from the InChI-string in importModel. This may have to change!
0151         
0152         hasInchi=false;
0153         if isfield(model,'inchis')
0154             if ~isempty(model.inchis{i})
0155                 hasInchi=true;
0156             end
0157         end
0158         
0159         if ~isempty(model.metFormulas{i}) && hasInchi==false
0160             toprint=[toprint '<notes><body xmlns="http://www.w3.org/1999/xhtml"><p>FORMULA: '  model.metFormulas{i} '</p></body></notes>'];
0161         end
0162     end
0163     
0164     %Only print annotations for metabolites with some miriam link. This is because I don't know how "unknown"
0165     %metabolites should be presented, and it seems unlikely that you will
0166     %have InChI without a database link. This might be temporary....
0167     if isfield(model,'metMiriams')
0168         miriamString=getMiriam(model.metMiriams{i});
0169     else
0170         miriamString=[];
0171     end
0172 
0173     if ~isempty(miriamString)
0174         toprint=[toprint '<annotation>'];
0175 
0176         %Print InChI if available
0177         if isfield(model,'inchis')
0178            if ~isempty(model.inchis{i})
0179                 toprint=[toprint '<in:inchi xmlns:in="http://biomodels.net/inchi" metaid="metaid_M_' model.mets{i} '_inchi">InChI=' model.inchis{i} '</in:inchi>'];
0180                 isInchi=1;
0181            else
0182                 isInchi=0;
0183            end
0184         else
0185             isInchi=0;
0186         end
0187         %Print some more annotation stuff
0188         toprint=[toprint '<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/" '...
0189             '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/">'...
0190             '<rdf:Description rdf:about="#metaid_M_' model.mets{i} '">'...
0191             '<bqbiol:is>'...
0192             '<rdf:Bag>'];
0193         if isInchi==1
0194             toprint=[toprint '<rdf:li rdf:resource="#metaid_M_' model.mets{i} '_inchi" />'];
0195         end
0196 
0197         %Print Miriam and finish up
0198         toprint=[toprint miriamString '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0199     end
0200     toprint=[toprint '</species>\n'];
0201     fprintf(fid,toprint);
0202 end
0203 
0204 %Write genes and complexes
0205 %Write genes
0206 if isfield(model,'genes')
0207     for i=1:numel(model.genes)
0208        toprint=['<species metaid="metaid_E_' int2str(i) '" id="E_' int2str(i) '" name="' model.genes{i} '" compartment="C_' geneComps{i} '" sboTerm="SBO:0000014">'];
0209 
0210        %Print gene name if present
0211        if isfield(model,'geneShortNames')
0212            if ~isempty(model.geneShortNames{i})
0213                 toprint=[toprint '<notes><body xmlns="http://www.w3.org/1999/xhtml"><p>SHORT NAME: '  model.geneShortNames{i} '</p></body></notes>'];
0214            end
0215        end
0216 
0217        if isfield(model,'geneMiriams')
0218             miriamString=getMiriam(model.geneMiriams{i});
0219        else
0220             miriamString=[]; 
0221        end
0222 
0223        if ~isempty(miriamString)
0224            toprint=[toprint '<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:dc="http://purl.org/dc/elements/1.1/" '...
0225                 '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/" '...
0226                 'xmlns:bqmodel="http://biomodels.net/model-qualifiers/"><rdf:Description rdf:about="#metaid_E_' int2str(i) '"><bqbiol:is><rdf:Bag>'...
0227                 miriamString '</rdf:Bag></bqbiol:is></rdf:Description></rdf:RDF></annotation>'];
0228        end
0229 
0230        toprint=[toprint '</species>\n'];
0231        fprintf(fid,toprint);
0232     end
0233 end
0234 %Write complexes (NOT DONE YET)
0235 
0236 %Finish metbolites
0237 fprintf(fid,'</listOfSpecies>');
0238 
0239 %Add reactions
0240 fprintf(fid,'<listOfReactions>');
0241 
0242 for i=1:length(model.rxns)
0243     %Get reversibility
0244     reversible='false';
0245     if model.rev(i)==1
0246         reversible='true';
0247     end
0248 
0249     %TEMP, MUST FIX SBO-TERM
0250     fprintf(fid,['<reaction metaid="metaid_R_' model.rxns{i} '" id="R_' model.rxns{i} '" name="' model.rxnNames{i} '" reversible="' reversible '" sboTerm="SBO:0000176">']);
0251 
0252     if isfield(model,'subSystems')
0253         subsystem=model.subSystems{i};
0254     else
0255         subsystem=[];
0256     end
0257     
0258     if ~isempty(subsystem)
0259         fprintf(fid,['<notes><body xmlns="http://www.w3.org/1999/xhtml"><p>SUBSYSTEM: ' subsystem '</p></body></notes>']);
0260     end
0261     
0262     %Print annotation
0263     if isfield(model,'rxnMiriams')
0264         miriamString=getMiriam(model.rxnMiriams{i});
0265     else
0266         miriamString=[];
0267     end
0268 
0269     %Check to see if there is an associated ec-code
0270     if isfield(model,'eccodes')
0271         if ~isempty(model.eccodes{i})
0272             miriamString=[miriamString  '<rdf:li rdf:resource="urn:miriam:ec-code:' model.eccodes{i} '"/>'];
0273         end
0274     end
0275 
0276     if ~isempty(miriamString)
0277         fprintf(fid,['<annotation><rdf:RDF xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" '...
0278         'xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:dcterms="http://purl.org/dc/terms/" '...
0279         'xmlns:vCard="http://www.w3.org/2001/vcard-rdf/3.0#" xmlns:bqbiol="http://biomodels.net/biology-qualifiers/" '...
0280         'xmlns:bqmodel="http://biomodels.net/model-qualifiers/">'...
0281         '<rdf:Description rdf:about="#metaid_R_' model.rxns{i} '">'...
0282         '<bqbiol:is>'...
0283         '<rdf:Bag>'...
0284         miriamString... 
0285         '</rdf:Bag>'...
0286         '</bqbiol:is>'...
0287         '</rdf:Description>'...
0288         '</rdf:RDF>'...
0289         '</annotation>']);
0290     end
0291     
0292     %The reactants have negative values in the stochiometric matrix
0293     compounds=model.S(:,i);
0294     reactants=find(compounds<0);
0295     products=find(compounds>0);
0296     
0297     if any(reactants)
0298         fprintf(fid,'<listOfReactants>');
0299         for j=1:length(reactants)
0300             tempmetname=model.mets{reactants(j)};
0301             fprintf(fid,['<speciesReference species="M_' tempmetname '" stoichiometry="' num2str(-1*compounds(reactants(j))) '"/>']);
0302         end  
0303 
0304         fprintf(fid,'</listOfReactants>');
0305     end
0306     if any(products)
0307         fprintf(fid,'<listOfProducts>');
0308         for j=1:length(products)
0309            tempmetname=model.mets{products(j)};
0310            fprintf(fid,['<speciesReference species="M_' tempmetname '" stoichiometry="' num2str(compounds(products(j))) '"></speciesReference>']);
0311         end  
0312         fprintf(fid,'</listOfProducts>');
0313     end
0314     
0315     genes=find(model.rxnGeneMat(i,:));
0316     if any(genes)
0317         fprintf(fid,'<listOfModifiers>');
0318         for j=1:numel(genes)
0319             fprintf(fid,['<modifierSpeciesReference species="E_' num2str(genes(j)) '" />']);
0320         end
0321         fprintf(fid,'</listOfModifiers>');
0322     end
0323 
0324     %Print constraints, reversibility, and objective. It's assumed that all
0325     %information is present in the model structure. This should be ok, since
0326     %it should have been set to standard values when imported
0327     
0328     %Note that the order of parameters is hard-coded. It's the same thing
0329     %in importModel
0330 
0331     fprintf(fid,'<kineticLaw><math xmlns="http://www.w3.org/1998/Math/MathML"><ci>FLUX_VALUE</ci></math><listOfParameters>');
0332     fprintf(fid,['<parameter id="LB_R_' model.rxns{i} '" name="LOWER_BOUND" value="' sprintf('%15.8f',model.lb(i)) '" units="mmol_per_gDW_per_hr"/><parameter id="UB_R_' model.rxns{i} '" name="UPPER_BOUND" value="' sprintf('%15.8f',model.ub(i)) '" units="mmol_per_gDW_per_hr"/><parameter id="OBJ_R_' model.rxns{i} '" name="OBJECTIVE_COEFFICIENT" value="' sprintf('%15.8f',model.c(i)) '" units="dimensionless"/><parameter id="FLUX_VALUE" value="0.00000000" units="mmol_per_gDW_per_hr"/>']);
0333     fprintf(fid,'</listOfParameters></kineticLaw>');
0334     fprintf(fid,'</reaction>');
0335 end
0336 
0337 %Write outro
0338 outro='</listOfReactions></model></sbml>';
0339 fprintf(fid,outro);
0340 
0341 fclose(fid);
0342 
0343 %Replace the target file with the temporary file
0344 delete(fileName);
0345 movefile(tempFile,fileName);
0346 end
0347 
0348 function miriamString=getMiriam(miriamStruct)
0349 %Returns a string with list elements for a miriam structure ('<rdf:li
0350 %rdf:resource="urn:miriam:obo.go:GO:0005739"/>' for example). This is just
0351 %to speed ut things since this is done many times during the exporting
0352 
0353 miriamString='';
0354 if isfield(miriamStruct,'name')
0355     for i=1:numel(miriamStruct.name)
0356         miriamString=[miriamString '<rdf:li rdf:resource="urn:miriam:' miriamStruct.name{i} ':' miriamStruct.value{i} '"/>'];
0357     end
0358 end
0359 end

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005