SysBioChalmers / Human-GEM

The generic genome-scale metabolic model of Homo sapiens
https://sysbiochalmers.github.io/Human-GEM-guide/
Creative Commons Attribution 4.0 International
98 stars 41 forks source link

Pseudoreactions in humanGEM #28

Closed haowang-bioinfo closed 5 years ago

haowang-bioinfo commented 6 years ago

Description of the issue:

Expected feature/value/output:

Current feature/value/output:

Format COBRA RAVEN
Exchange 0 461+1639-439
Sink 96 (removed) 1
Demand 157 (removed) 0

I hereby confirm that I have:

haowang-bioinfo commented 6 years ago

Sink and demand (DM) reactions are both Pseudoreactions that are formulated for modeling purposes, i.e. consumption and biosynthesis of required metabolites.

By definition, Sink reaction are always reversible while demand reactions are usually irreversible. Given its irreversible direction, should sink_citr[c] be renamed to DM_citr[c]?

haowang-bioinfo commented 6 years ago

In HMR2, pseudoreactions are categorized into "Exchange reactions" subSystem, even though there was one sink reaction HMR_9736 (cholesterol-ester pool[x] <=> cholesterol-ester pool[l]). Recon3D classifies pseudoreactions into "Exchange/demand reaction" subSystem, that is better than just covering exchange ones.

@pecholleyc @pinarkocabas @JonathanRob Please provide your comments regarding this. My view is to apply the latter as subSystem name of pseudoreactions in humanGEM.

haowang-bioinfo commented 6 years ago

Another reaction 'DM_atp_c_' (ATP[c] + H2O[c] => ADP[c] + H+[c] + Pi[c]) appears to be problematic. It doesn't look like a demand reaction. Instead it seems to be a reaction catalyzed by ATPase. Since there is already ATPase reaction 'ATPasel' (ATP[c] + H2O[c] => ADP[c] + H+[l] + Pi[c]), the equation of 'DM_atp_c_' should be wrong. And the correct equation suppose to be (ATP[c] ⇌), according to BiGG database.

Their reaction names (Demand for ATP, Cytosolic and V-Type ATPase, H+ Transporting, Lysosomal) supported this deduction.

pinarkocabas commented 6 years ago

I think we can also use "Exchange/demand reaction" subSystem.

haowang-bioinfo commented 6 years ago

The 253 sink and demand reactions have been converted from COBRA to RAVEN format, by adding corresponding boundary metabolites through using addBoundaryMets.m function in commit 7fec699. Now all pseudoreactions in humanGEM are unified into the same format.

JonathanRob commented 6 years ago

This is nice, I agree with the changes you have made/suggested so far.

To take this a step further (or in a slightly different direction), what do you think of simply deleting all demand and sink reactions from humanGEM? Since they are artificial, and HMR2 didn't contain any such reactions, I don't think we need them. In addition, they may lead to erroneous results in other situations, such as checking metabolic tasks and running tINIT.

haowang-bioinfo commented 6 years ago

HMR2 didn't contain any such reactions

This is not accurate. HMR2 does contain such pseudoreactions, including 461 exchange type and 1 sink type HMR_9736.

In addition, there are 1639 exchange reactions that are included in Recon3D and then introduced into humanGEM v0.2.0. Later in v0.3.0, a total of 439 Recon3D exchange reactions were detected as duplications and removed .

haowang-bioinfo commented 6 years ago

they may lead to erroneous results in other situations, such as checking metabolic tasks and running tINIT.

You probably are referring to sink/DM reactions here. It would be good to have some evidence to exemplify the problem-causing cases.

JonathanRob commented 6 years ago

I was a bit unclear in my previous comment.

HMR2 didn't contain any such reactions

Here, I was referring specifically to sink/demand reactions, not exchange reactions. As far as I'm aware, the HMR_9736 reaction was a mistake, and was intended to be an exchange reaction (but I could be wrong). Regardless, it is clear that these sink/demand type reactions were not something that the designers of HMR felt were something that should be included in the model.

To elaborate on my previous comment, I think these sink/demand reactions should be removed, because they are not biologically relevant, and are unnecessary. I understand that such reactions can be useful in model curation purposes, e.g., to allow consumption/production of an intracellular metabolite without restriction. However, I think these reactions should only be added temporarily by the user (or they can instead modify the b-vector to allow for metabolite mass imbalances) for those particular applications, and should not be a permanent fixture in the model.

When I referred to potential problems with the sink/demand reactions, I was thinking of cases where a user may be constraining the exchange reactions for a particular condition they would like to simulate, but are unaware of the sink/demand reactions that allow for metabolites to be synthetically created or destroyed in their simulation (if the reactions are left "open"). Of course, they should be constraining the sink/demand reactions as well, but it leaves more opportunities for such problems to arise.

haowang-bioinfo commented 6 years ago

As far as I'm aware, the HMR_9736 reaction was a mistake, and was intended to be an exchange reaction (but I could be wrong).

It would worth to setup an independent issue regarding reaction HMR_9736 (cholesterol-ester pool[x] <=> cholesterol-ester pool[l]), which seems to be problematic.

JonathanRob commented 5 years ago

A script (curateExchangeRxns.m) was generated to address some of the points in this issue. Here is a description of the purpose/content of that script:

Script to curate exchange, sink, and demand reactions. The current model has a mix of exchange reactions that are formulated such that negative flux corresponds to export:

HMR_9730: albumin[x] <=> albumin[s]

whereas others have positive flux corresponding to export:

EX_adrn[e]: adrenic acid[s] <=> adrenic acid[x]

In order to standardize this, all exchange reactions will be converted to the second format (where positive flux = export).

However, before this can be incorporated, there are 7 exchange reactions in the model that are duplicated but are written in opposite directions:

HMR_9025: VLDL[x] => VLDL[s] HMR_9049: VLDL[s] => VLDL[x]

HMR_9026: HDL[x] => HDL[s] HMR_9050: HDL[s] => HDL[x]

HMR_9027: LDL[x] => LDL[s] HMR_9051: LDL[s] => LDL[x]

HMR_9028: chylomicron remnant[x] => chylomicron remnant[s] HMR_9052: chylomicron remnant[s] => chylomicron remnant[x]

HMR_9029: VLDL remnant[x] => VLDL remnant[s] HMR_9053: VLDL remnant[s] => VLDL remnant[x]

HMR_9030: HDL remnant[x] => HDL remnant[s] HMR_9054: HDL remnant[s] => HDL remnant[x]

HMR_9031: LDL remnant[x] => LDL remnant[s] HMR_9055: LDL remnant[s] => LDL remnant[x]

These reactions must first be merged, otherwise switching reaction directionality will result in these rxn pairs being completely identical. Therefore, for each pair, one reaction was deleted, and the second was made reversible. It was confirmed in advance that none of these reactions are associated with any genes.

The script also inactivates all sink and demand (DM) reactions by constraining their upper and lower bounds to zero. These reactions will be considered for full deletion in future model versions.

Finally, the upper and lower bounds of all exchange reactions are set to +/-1000, respectively. As a result, the model is by default completely "open", allowing free exchange of all metabolites.

haowang-bioinfo commented 5 years ago

Good work @JonathanRob. The re-formulation of exchange reactions (i.e. [s] <=> [x]) is a nice move, which improves the model compatibility with the COBRA toolbox.

Those Sink and DM reactions in current model don't make sense to me, and appear to be wrong. So it is safe to constrain them now, and check alternative routes (e.g. [c] <=> [x] vs [c] <=> [s] and [s] <=> [x]) in the model before removal.

JonathanRob commented 5 years ago

I would like to bring some attention back to this issue, specifically regarding the "sink" and "demand" reactions. These reactions are recognized by some functions (e.g., getExchangeRxns) as exchange reactions, which is incorrect, and can lead to problems/errors.

Currently, these sink and demand reactions are inactivated, meaning their upper and lower bounds are set to zero. This was done instead of a full removal, as we wanted to avoid accidental removal of reactions because it's difficult to undo this (see issue #44). However, this is a bit of a special case, because it is very often that exchange reactions have their bounds set to zero to simulate different nutrient environments. The sink and demand reactions, on the other hand, have their bounds set to zero because they are scheduled for eventual removal from the model.

This creates a confusing situation for functions, because they are unable to tell the difference between an exchange reaction that is simply turned off, and a sink/demand reaction that is inactivated.

Therefore, I recommend that we fully remove all sink_ and DM_ reactions from humanGEM in the next update (or sometime very soon).

haowang-bioinfo commented 5 years ago

The "sink" and "demand" reactions are recognized by some functions (e.g., getExchangeRxns) as exchange reactions, which is incorrect, and can lead to problems/errors.

It appears to me that these functions, such as getExchangeRxns, need to be adjusted to better recognize exchange reactions among others. Because "sink" and "demand" reactions do occur in GEMs.

I recommend that we fully remove all sink and DM reactions from humanGEM in the next update (or sometime very soon).

Since it's certain that sink and DM reactions don't fit in humanGEM, I agree to get rid of them soon.

haowang-bioinfo commented 5 years ago

Now we've reached a consensus regarding the long-lasting discussion upon Sink/DM reactions. As a generic model, humanGEM will NOT include them anymore. However, these reactions can be added by the users for their particular applications.