Tutorial 5: Reconstruct a GEM from KEGG#
The exercise creates a model from KEGG on the basis of protein sequences in a FASTA file and runs functionality checks on the result. The example organism is Saccharomyces cerevisiae.
Unlike Tutorials 1 to 4, this tutorial is a demonstration. It provides a scaffold that can be adapted to reconstruct a GEM for any organism.
The de novo reconstruction searches large HMM databases. Building the model takes 20 to 35 minutes on macOS and Unix systems and 40 to 55 minutes on Windows, depending on the hardware and the size of the proteome of the target organism. gapReport can take several to many hours, depending on the number of gaps in the model.
Step by step#
1. Reconstructing from KEGG#
getKEGGModelForOrganism needs trained hidden Markov models (HMMs) for eukaryotes. When the library file <dataDir>.hmm or its .gz archive exists, for example after a manual download, it is used. Otherwise it is downloaded from the matching release of toolbox-assets. Here the archive kegg118_eukaryotes is selected for automatic download.
The function then creates a model for S. cerevisiae from the FASTA file of the whole proteome (sce.fa). The parameters exclude general or unclear reactions and reactions with undefined stoichiometry. help getKEGGModelForOrganism describes the parameters.
model=getKEGGModelForOrganism('sce','sce.fa','kegg118_eukaryotes','output',false,false,false,false,10^-30,0.8,0.3,-1);
disp(model);
The resulting model should contain around 1589 reactions, 1600 metabolites and 836 genes. Small variations are possible because the algorithm is heuristic, and different KEGG versions give slightly different results.
2. Removing reactions that make something from nothing#
A first control is that the model cannot produce any metabolite without the uptake of some metabolites. A violation often arises when metabolites have a different meaning in different reactions. removeBadRxns finds and removes such reactions automatically, and help removeBadRxns gives the details.
[newModel, removedRxns]=removeBadRxns(model);
The output reports that H⁺ can be made even when no reaction was unbalanced. Protons are problematic because the pH at which the formulas are written is arbitrary. For this analysis the protons can be ignored and corrected later, so repeat the call with the production of H⁺ allowed.
[newModel, removedRxns]=removeBadRxns(model,1,{'H+'},true);
disp(removedRxns);
One reaction was removed because it enabled the model to produce something from nothing. As it is a single reaction, examine it in more detail.
3. Investigating the problematic reaction#
According to KEGG, the removed reaction is a general polymer reaction. findLeakMetabolite shows the flux distributions in more detail and shows whether a better alternative exists for deletion.
[fluxes, metabolite]=findLeakMetabolite(model,'produce',{'H+'},true);
model.metNames(metabolite)
printFluxes(model, fluxes, false, [], [],'%rxnID (%rxnName):\n\t%eqn: %flux\n')
The model could produce H₂O, and the output contains more fluxes than can be read at once. Excluding the elementally balanced reactions shortens the list. Water was produced, so the reactions that are unbalanced for oxygen (column 6 of the elemental balance) are of interest.
balanceStructure=getElementalBalance(model);
goodOnes=balanceStructure.leftComp(:,6)==balanceStructure.rightComp(:,6);
printFluxes(removeReactions(model,goodOnes), fluxes(~goodOnes), false, [], [],'%rxnID (%rxnName):\n\t%eqn: %flux\n');
Many reactions remain. The reactions that involve amylose or starch are kept, as one of the problematic reactions from step 2 involves them.
printFluxes(model, fluxes, false, [], [],'%rxnID (%rxnName):\n\t%eqn: %flux\n',{'Amylose';'Starch'});
Two elementally unbalanced reactions remain, including the reaction identified by removeBadRxns. They contradict each other: the first shows that amylose and starch are interconvertible, and the second shows that amylose contains one glucose unit fewer than starch. General reactions of this kind are fixed manually. The result of removeBadRxns is trusted, and R02110 is deleted.
model=removeReactions(model,'R02110');
4. Checking consumption and adding uptakes#
After the removal of R02110, the model cannot make something from nothing. The next check is whether it can consume something without any output.
[solution, metabolite]=findLeakMetabolite(model,'consume',{'H+'},true);
model.metNames(metabolite)
Nothing is consumed without output. Add some uptakes and check the products of the model.
[~, J]=ismember({'D-Glucose';'H2O';'Orthophosphate';'Oxygen';'NH3';'Sulfate'},model.metNames);
[model, addedRxns]=addExchangeRxns(model,'in',J);
canExchange reports the metabolites that can be produced with these uptakes. It allows the output of all metabolites, which does not happen in a real cell, and it serves functionality testing.
I=canExchange(model,'produce');
fprintf('%d%%\n', round(sum(I)/numel(model.mets)*100));
About 31 % of the metabolites can be synthesised. The result does not show whether this is high or low. Many metabolites are not expected to be synthesised from glucose, water, orthophosphate, oxygen, NH3 and sulfate alone.
5. Filling gaps with the full KEGG model#
Gap-filling with the full KEGG model tests whether it gives a higher fraction of synthesisable metabolites. getModelFromKEGG retrieves the full KEGG model. It is associated with more than 6,400,000 genes that gap-filling does not use, and they are removed to speed up the step.
keggModel=getModelFromKEGG(false,false,false,false);
keggModel=rmfield(keggModel,'genes');
keggModel=rmfield(keggModel,'rxnGeneMat');
Some reactions in KEGG are known to be unbalanced, and the balanced reactions are used for gap-filling.
balanceStructure=getElementalBalance(keggModel);
keggModel=removeReactions(keggModel,balanceStructure.balanceStatus~=1,true,true);
With these settings, fillGaps tries to include reactions that give flux through all reactions in the model. The first flag allows the production of all metabolites.
params.relGap=0.6; %Lower number for a more exhaustive search
params.printReport=true;
[newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,keggModel,true,false,false,[],params);
fillGaps connected around 29 reactions (newConnected) by including around 41 reactions from the KEGG model (addedRxns). These reactions should be checked manually to confirm that they exist in yeast, and here they are assumed to.
6. Reporting connectivity#
The connectivity of the model is improved further by identifying metabolites that should be connected. gapReport summarises how connected the model is, together with the supporting detail.
[noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, minToConnect,...
neededForProductionMat]=gapReport(newModel);
532 of 1634 reactions cannot carry flux. Output is allowed for all metabolites, so the number is 532 in both cases. 543 of 1617 metabolites cannot be synthesised from the supplied precursors. The model has 6 subnetworks, and 1605 of 1617 metabolites belong to the first.
gapReport also prints the metabolites to connect. Net production of all metabolites requires the connection of 322 metabolites, and the top-ranked one is Acyl-CoA, which connects 162 metabolites. The top 10 of the list connect 271 of 543 (50 %) of them. The list shows where the gap-filling has the greatest effect.
7. Adding minimal-media uptakes and iterating#
The list shows that yeast cannot grow on the substrates defined so far and requires other precursors for cofactor synthesis. Add uptake reactions for the minimal-medium constituents that yeast needs to grow.
[~, J]=ismember({'4-Aminobenzoate';'Riboflavin';'Thiamine';'Biotin';'Folate';'Nicotinate';'Zymosterol';'Choline'},newModel.metNames);
[newModel, addedRxns]=addExchangeRxns(newModel,'in',J);
Run gapReport again, and use its output to direct the gap-filling. gapReport prints part of the information, and the output structures hold the rest. Repeat these steps until the model meets the quality criteria of the project, for example a fraction of synthesisable metabolites.
Full script#
% tutorial5
% This exercise is about creating a model from KEGG, based on protein
% sequences in a FASTA file, and doing some functionality checks on the
% model. The example case is for the yeast Saccharomyces cerevisiae. This
% tutorial is more of a showcase than the previous four, and its main
% purpose is to serve as a scaffold to reconstruct a GEM for any
% organism.
% This refers to Tutorial 5 on the RAVEN wiki:
% https://raven-docs.readthedocs.io/en/latest/tutorials/
%
% Start by downloading trained Hidden Markov Models for eukaryotes. This can
% be done automatically or manually from the RAVEN Wiki in its GitHub
% repository. In this tutorial, the archive "kegg118_eukaryotes" is picked for
% the automatic download. See the documentation in the RAVEN Wiki for more
% information regarding preparation of such archive.
%
% This creates a model for S. cerevisiae. The parameters are set to exclude
% general or unclear reactions and reactions with undefined stoichiometry.
% Type "help getKEGGModelForOrganism" to see what the different parameters
% are for. This process takes up to 20-35 minutes in macOS, Unix systems and
% 40-55 minutes in Windows, depending on your hardware and the size of
% target organism proteome
model=getKEGGModelForOrganism('sce','sce.fa','kegg118_eukaryotes','output',false,false,false,false,10^-30,0.8,0.3,-1);
% The resulting model should contain around 1589 reactions, 1600
% metabolites and 836 genes. Small variations are possible since it is an
% heuristic algorithm and different KEGG versions will give slightly
% different results.
disp(model);
% A first control is that the model should not be able to produce any
% metabolites without uptake of some metabolites. This commonly happens when
% metabolites have a different meaning in different reactions. The best way
% to find such reactions is to run findLeakMetabolite and analyze the resulting
% solution for bad reactions. An automated approach is to use removeBadRxns,
% which tries to do the same thing in an automated manner. Type
% "help removeBadRxns" for details.
[newModel, removedRxns]=removeBadRxns(model);
% One can see an error about that H+ can be made even if no reactions were
% unbalanced. Protons are particularly problematic since it is rather
% arbitary at which pH the formulas are written for. For the purpose of this
% analysis, the protons can be ignored and fixed later.
[newModel, removedRxns]=removeBadRxns(model,1,{'H+'},true);
% Only one reaction was removed because it enabled the model to produce
% something from nothing. Since it is only one reaction, it might be
% worthwhile to look into this in more detail.
disp(removedRxns);
% According to the information in KEGG about this reaction it is a general
% polymer reaction. One might want to look at the flux distributions
% in more detail to try to find out if there is any better alternative to
% delete. Use findLeakMetabolite to do this.
[fluxes, metabolite]=findLeakMetabolite(model,'produce',{'H+'},true);
model.metNames(metabolite)
% The model could produce H2O using the following reactions
printFluxes(model, fluxes, false, [], [],'%rxnID (%rxnName):\n\t%eqn: %flux\n')
% That resulted in quite a lot of fluxes to look through. Maybe it is easier
% if the elementally balanced reactions are excluded.
balanceStructure=getElementalBalance(model);
% The hydrogen balancing is a bit tricky, so this time it seems more
% relevant to only look at the ones unbalanced for oxygen (since water was
% produced)
goodOnes=balanceStructure.leftComp(:,6)==balanceStructure.rightComp(:,6);
printFluxes(removeReactions(model,goodOnes), fluxes(~goodOnes), false, [], [],'%rxnID (%rxnName):\n\t%eqn: %flux\n');
% There is still a good number of reactions. Leave only the reactions which
% involve amylose or starch, from one of the problematic reactions that was
% identified earlier.
printFluxes(model, fluxes, false, [], [],'%rxnID (%rxnName):\n\t%eqn: %flux\n',{'Amylose';'Starch'});
% There are two elementally unbalanced reactions, including, the reaction
% which was identified by removeBadRxns. When looking to these reactions
% closer, one can notice the contradiction between the reactions. The first
% one shows that amylose and starch are interconvertible, while the second
% reaction shows that amylose contains one less glucose unit than starch.
% This type of general reactions are problematic and should be fixed
% manually. It is hereby chosen to trust removeBadRxns and delete R02110.
model=removeReactions(model,'R02110');
% The model can no longer make something from nothing. Can it consume
% something without any output?
[solution, metabolite]=findLeakMetabolite(model,'consume',{'H+'},true);
model.metNames(metabolite)
% Nope, so that was good. Add some uptakes and see what it can produce.
[~, J]=ismember({'D-Glucose';'H2O';'Orthophosphate';'Oxygen';'NH3';'Sulfate'},model.metNames);
[model, addedRxns]=addExchangeRxns(model,'in',J);
% Check which metabolites can be produced given these uptakes. The
% canExchange function allows for output of all metabolites. This will not
% happen in the real cell, but it is very useful for functionality testing
% of the model. Once it is functional, the excretion reactions based on
% evidence can be added as well.
I=canExchange(model,'produce');
fprintf('%d%%\n', round(sum(I)/numel(model.mets)*100));
% It seems that around 31% of the metabolites could be synthesized. It is
% not directly clear whether this is a high or low number, many metabolites
% should not be possible to synthesize from those simple precursors.
% Try to fill gaps using the full KEGG model to see if that gives a
% significantly higher number
keggModel=getModelFromKEGG(false,false,false,false);
% The KEGG model is associated to more than 6,400,000 genes. They will not
% be used for the gapfilling, so they are removed to make this a little
% faster
keggModel=rmfield(keggModel,'genes');
keggModel=rmfield(keggModel,'rxnGeneMat');
% It is already known from the previous commands that there are some
% unbalanced reactions in KEGG. Only use the balanced ones for the gap
% filling
balanceStructure=getElementalBalance(keggModel);
keggModel=removeReactions(keggModel,balanceStructure.balanceStatus~=1,true,true);
% The function fillGaps with these settings will try to include reactions in
% order to have flux through all reactions in the model. There are other
% settings as well. The first flag says that production of all metabolites
% should be allowed.
params.relGap=0.6; %Lower number for a more exhaustive search
params.printReport=true;
[newConnected, cannotConnect, addedRxns, newModel, exitFlag]=fillGaps(model,keggModel,true,false,false,[],params);
% The results show that fillGaps could connect around 29 reactions
% (newConnected) by including around 41 reactions from the KEGG model
% (addedRxns). Those should of course be checked manually to confirm that
% they exist in yeast, but in this tutorial it is assumed that all these
% reactions indeed occur in yeast.
% Continue to improve the connectivity of the model by identifying
% metabolites that should be connected. A convenient way to get an overview
% of how connected the model is, and at the same time getting a lot of
% useful data, is to use gapReport. Note that it can take several to many
% hours to run, depending on the number of gaps in the model.
[noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, minToConnect,...
neededForProductionMat]=gapReport(newModel);
% The results show that 532 from 1634 reactions cannot carry flux. Remember
% that an output is allowed for all metabolites, that is why it calculates
% 532 in both cases. 543 from 1617 metabolites cannot be synthesized from
% the precursors which are supplied. There are 6 subnetworks in the model,
% of which 1605 from 1617 metabolites belong to the first one.
% It will also print something similar to:
% To enable net production of all metabolites, a total of 322 metabolites
% must be connected
% Top 10 metabolites to connect:
% 1. Acyl-CoA[s] (connects 162 metabolites)
% 2. Thiamin diphosphate[s] (connects 17 metabolites)
% 3. G00003[s] (connects 14 metabolites)
% 4. Androstenediol[s] (connects 13 metabolites)
% 5. [Dihydrolipoyllysine-residue succinyltransferase] S-glutaryldihydrolipoyllysine[s] (connects 13 metabolites)
% 6. Selenomethionyl-tRNA(Met)[s] (connects 12 metabolites)
% 7. Retinyl ester[s] (connects 11 metabolites)
% 8. 1-Alkyl-2-acylglycerol[s] (connects 10 metabolites)
% 9. G00146[s] (connects 10 metabolites)
% 10. Progesterone[s] (connects 9 metabolites)
% This is a very useful way of directing the gap-filling tasks to where they
% are of the greatest use. In this case it says that to have net
% production of all metabolites in the model from the simple inputs which
% have been defined (543 metabolites can currently not have net production)
% one needs to connect 322 unconnected metabolites in. However, by
% connecting only the top 10 in the list 271/543 (50%) of them will now be
% connected. When looking at the list one can see that the top ranking one
% is Acyl-CoA, which represents many different fatty acid chains. The second
% metabolite, thiamine disphosphate, is an enzyme co-factor. It turns out
% that yeast cannot grow on only the substrates that had been defined, but
% that it requires some other precursors for co-factor synthesis as well.
% Add uptake reactions for the minimal media constituents needed for yeast
% to grow.
[~, J]=ismember({'4-Aminobenzoate';'Riboflavin';'Thiamine';'Biotin';'Folate';'Nicotinate';'Zymosterol';'Choline'},newModel.metNames);
[newModel, addedRxns]=addExchangeRxns(newModel,'in',J);
% Rerun gapReport and use the output for targeting the gap-filling efforts.
% Note that only some info is printed; most of it is available in the output
% structures. Work like this in an iterative manner until the model is of
% sufficient quality.