detectDeadEnds has buggy behavior when removeExternalMets=1 about cobratoolbox HOT 14 CLOSED

@akaviaLab , may we close this issue?

from cobratoolbox.

Well, what should happen with bvite in this example?
The code hasn't changed since 2015. This is still a bug.

from cobratoolbox.

Hi,
The documented behaviour for removeExternalMets is that it should "remove metabolites that participate in reactions only with themselves". To me there are two interpretations to this statement:
A: Those metabolites should not be returned (which is exactly what the function currently does)
B: Those metabolites should be removed, before checking for dead ends.

Since A is the implemented behaviour, I would assume, that this is what is desired.

I personally think, that it is intended behaviour that it removes everything that has an exchanger (yes, this also includes explicit demand reactions in the cytosol, which might not be intentional), as that is what the documentation says it should do.

So I don't really see the "bug" here.

About your example:
Both strch1[e] and bvite[c] and bvite[e] should be removed from the list if the option is set (which currently happens), since they both have "reactions with only themselves".

One might ask the question, whether it makes sense to have this option, and/or whether an option to ignore Exchange Reactions during determination of dead ends would make sense, but from the function behaviour I don't see a bug.

from cobratoolbox.

from cobratoolbox.

Okay - I think we might have been using different definitions of "only with itself." I agree with you about bvite now, but strch1 participates in the reaction AMY1e
which has the formula h2o[e] + strch1[e] => 8 glc_D[e] + strch2[e]

In that case, I would argue that strch1 participates in reactions with other metabolites as well, and it shouldn't be removed. Is that correct? If not, can you clarify your definition of "only with itself"?

I agree with Ronan's comment about doing it algorithmically. That sounds a potentially good idea.

from cobratoolbox.

Ahh ok, Thats where you are coming from.
I didn't read this exclusively i.e. my interpretation was that metabolites that "participate in reactions only with themselves" can also participate in reactions with other metabolites. Thats what the code is looking at. "Does the metabolite exist in an exchange reaction? Yes: Remove it. No: Keep it" (with the definition of an Exchange reaction being that it only has a single metabolite).

If you read it as metabolites that only participate in reactions where they are the exclusive metabolite, I would say that both bvite and strch are definitely not part of that set, so they should be returned (for bvite[e]/bvite[c] they both are in the transport reaction from e to c, which is a reaction containing another metabolite).

from cobratoolbox.

from cobratoolbox.

I have been thinking about it a bit, and bvite should never be in the list anyways, as it has both reactions consuming and producing it (both e and c), so from the function description (check if values are all -1 or all +1 and also check if lb is zero or not) bvite[c] should never be added to the list, as it has coefficients of -1 and 1 and bvite[e] takes part in a reaction with a lb < 0 and thus can also be produced or consumed.

Importantly, the function itself does not find all "unproducible/consumable" metabolites but only those which are clear dead ends (i.e. don't have any reaction that produce them or don't have any reaction that consumes them).

I also noticed, that the function itself is extremely inefficient, and updated the function along with its documentation and submitted a pull request.

For the sake of simplicity (at least at the moment) I have only added Ronans suggestions in a commented part, as I am not quite sure, how an external metabolite would have to be defined in this instance? From my understanding, it is a metabolite which is either inconsistent according to the findStoichConsistentSubset function ( @rmtfleming , is there a clear definition that you use for those metabolites marked in InConsistentMetBool ?) or detected as exchanged metabolite by findSExRxnInd, but I currently lack a clear mathematical defintion of those. Is there a paper available?

from cobratoolbox.

from cobratoolbox.

Ok, do I get this right than that:
SInConsistentMetBool is 1 if the metabolite is present in any Reaction that is inconsistent?
How is stoichiometric consistency defined in this scenario? It seems not to be capability to carry flux (at least if I compare the SConsistentRxnBool with min/max fluxes of an FVA, it contains rxns flagged as consistent which can't carry flux).
I assume this is referring to the stoichiometric consistency defined in Gevorgyan et al:
https://academic.oup.com/bioinformatics/article/24/19/2245/248822/Detection-of-stoichiometric-inconsistencies-in

Is that correct?
In this instance I would adapt the description of the removeExternalMets option as follows:
Don't return metabolites that either participate in exchange reactions ("A <=> " or " <=> A") or are restricted to having 0 mass according to the definition of stoichiometric consistency in Gevorgyan et al, Bioinformatics 2008

Does that make sense?

from cobratoolbox.

from cobratoolbox.

from cobratoolbox.

Ok,
I updated the PR to reflect that and changed the effect of the optional argument accordingly.

from cobratoolbox.

PR #391 merged.

from cobratoolbox.