Potential New Automated Tests

[Devlin-Moyer]

I've been working on a Python package for automatically detecting several kinds of errors in GEMs as a formalization/generalization of much of the error correcting I've done in Human-GEM. It is now ready to use; I've successfully run it on every version of Human-GEM and thousands of other GEMs for different organisms, and made it pip-installable. When I first mentioned this a while ago, @haowang-bioinfo suggested that it could be turned into a GitHub action bot, and that I should start a discussion once I'd finished developing the tool, so here is that discussion.

[edkerk]

Nice, and being a Python package would make this relatively easy to implement. One issue I foresee though is what to do with the output. It currently generates a Pandas.Dataframe, but how do you imagine this could be represented once the action is run? Some options:

• Write a datafile that will be committed (possible in a dedicated branch), similar as is done for memote reports in yeast-GEM *1
• Prepare a commit text that is added to PR's, similar as is done for RAVEN. *2
• Just have the test results available via GitHub actions *2

*1 = Pandas.Dataframe should be exported to some form of flat-text or CSV format?
*2 = Instead of constructing a Pandas.Dataframe, the results should be output to command-line

I'm unaware of what the macaw output looks like, but I would imagine that *2 could be quite long, so perhaps *1 is more feasible? Can it provide summary output and a detailed output?

[Devlin-Moyer]

(it's not super relevant to GitHub actions for Human-GEM, but it's maybe interesting to note that the runtime for the dilution and loop tests strongly depends on the number of reactions in the model; it takes under an hour to run all tests on iML1515 (2,712 reactions) and 10-15 minutes on many of the AGORA2 models of various human gut microbes (generally 1,000-1,500 reactions))

[feiranl]

In that case, I would say we can put the dead-end and duplicate tests in the test, but the loop and dilution in each release, but I am not sure about the diphospahte, it seems not related to each PR and maybe can also be done in each release.

[Devlin-Moyer]

I'm not sure how it would be most helpful to set them up as tests that run for every PR, since, as you can see above, the current version of Human-GEM has many reactions flagged by every test except the diphosphate test (because I addressed all of those in #557 and #641). Also adding/creating new dead-ends might not necessarily be a problem, i.e. I'm not sure if it would be helpful to have a test "fail" a new PR if it increases the number of dead-ends. Similarly, the duplicate test doesn't just look for reactions that are completely identical to each other, but also reactions that are very similar that could represent an error but could be accurate; e.g. it will flag MAR09985 and MAR11307 as duplicates because they both move 5-aminolevulinate between [c] and [e] along with either 1 proton (MAR11307) or 2 protons (MAR09985), but SLC15A1 is known to be capable of doing both (but I only figured that out because I saw that the duplicate test flagged them and decided to investigate). It does seem like it would be nice to have GitHub warn you that a particular PR increases the number of dead-ends or duplicates so that you make sure that that was an anticipated consequence of those changes, but I don't know enough about how GitHub Actions work to figure out what the best way to set that up is

[edkerk]

I propose that test results do not throw a failure as GitHub Action, precisely for what you argue. Meanwhile, there would be some value to report such test results, but running all of it would probably fail on GitHub Actions (there is a 6-hour maximum runtime, and I doubt that it runs on 8+ threads with 256 GB RAM). As @feiranl suggested, a complete test could be nice to implement for each release, but this would be performed locally. For this purpose it might be useful to prepare a docker that can easily be deployed?

For the GitHub actions, we can just test how this looks like by implementing it in a feature branch:

• What would be the command to just run the fast tests, show a summary output in the terminal, and write the detailed results to an output file?

[Devlin-Moyer]

import cobra
from macaw.main import dead_end_test, duplicate_test

model = cobra.io.read_sbml_model('model.xml')
(dead_end_results, dead_end_edges) = dead_end_test(model)
(duplicate_results, duplicate_edges) = duplicate_test(model)
output = dead_end_results.merge(duplicate_results)
output.to_csv('macaw_dead_end_and_duplicate_test_results.csv', index = False)

Running the above as a Python script would print the following to the terminal if model.xml was version 1.18 of Human-GEM:

Starting dead-end test...
 - Found 1370 dead-end metabolites.
 - Found 1123 reactions incapable of sustaining steady-state fluxes in either direction due to these dead-ends.
 - Found 1594 reversible reactions that can only carry steady-state fluxes in a single direction due to dead-ends.
Starting duplicate test...
 - Skipping redox duplicates because no redox_pairs and/or proton_ids were provided.
 - Found 706 reactions that were some type of duplicate:
   - 0 were completely identical to at least one other reaction.
   - 13 involve the same metabolites but go in the opposite direction or have the opposite reversibility as at least one other reaction.
   - 445 involve the same metabolites but with different coefficients as at least one other reaction.

and macaw_dead_end_and_duplicate_test_results.csv would look like:

reaction_id	reaction_equation	dead_end_test	duplicate_test_exact	duplicate_test_directions	duplicate_test_coefficients	duplicate_test_redox
MAR03905	MAM01796c + MAM02552c --> MAM01249c + MAM02039c + MAM02553c	ok	ok	ok	ok	NA
MAR03907	MAM01796c + MAM02554c --> MAM01249c + MAM02039c + MAM02555c	ok	ok	ok	ok	NA
(one row for each reaction in the model)	⋮	⋮	⋮	⋮	⋮	⋮

Also testing for what I have dubbed "redox duplicates" would require providing a list of lists of pairs of metabolite IDs corresponding to the oxidized and reduced forms of things like NAD(H) and NADP(H), like so:

import cobra
from macaw.main import dead_end_test, duplicate_test

model = cobra.io.read_sbml_model('model.xml')
redox_pairs = [
    # NAD+ and NADH
    ('MAM02552c', 'MAM02553c'), ('MAM02552e', 'MAM02553e'),
    ('MAM02552m', 'MAM02553m'), ('MAM02552n', 'MAM02553n'),
    ('MAM02552r', 'MAM02553r'), ('MAM02552x', 'MAM02553x'),
    # NADP+ and NADPH
    ('MAM02554c', 'MAM02555c'), ('MAM02554l', 'MAM02555l'),
    ('MAM02554m', 'MAM02555m'), ('MAM02554n', 'MAM02555n'),
    ('MAM02554r', 'MAM02555r'), ('MAM02554x', 'MAM02555x'),
    # FAD and FADH2
    ('MAM01802c', 'MAM01803c'), ('MAM01802m', 'MAM01803m'),
    ('MAM01802r', 'MAM01803r'), ('MAM01802x', 'MAM01803x'),
    # FMN and FMNH2
    ('MAM01828r', 'MAM20019r'), ('MAM01828i', 'MAM20019i'),
    # ubiquinone and ubiquinol
    ('MAM03103c', 'MAM03102c'), ('MAM03103e', 'MAM03102e'),
    ('MAM03103m', 'MAM03102m'),
    # ferricytochrome c and ferrocytochrome c
    ('MAM01824c', 'MAM01826c'), ('MAM01824m', 'MAM01826m'),
    # ferricytochrome b5 and ferrocytochrome b5
    ('MAM01823c', 'MAM01825c'),
    # thioredoxin
    ('MAM02990c', 'MAM02666c'), ('MAM02990m', 'MAM02666m'),
    # mitothioredoxin
    ('MAM02487m', 'MAM02486m'),
    # oxygen and hydrogen peroxide
    ('MAM02630c', 'MAM02041c'), ('MAM02630e', 'MAM02041e'),
    ('MAM02630l', 'MAM02041l'), ('MAM02630m', 'MAM02041m'),
    ('MAM02630n', 'MAM02041n'), ('MAM02630r', 'MAM02041r'),
    ('MAM02630x', 'MAM02041x')
]
proton_ids = [m.id for m in model.metabolites if m.id.startswith('MAM02039')]

(dead_end_results, dead_end_edges) = dead_end_test(model)
(duplicate_results, duplicate_edges) = duplicate_test(model, redox_pairs, proton_ids)
output = dead_end_results.merge(duplicate_results)
output.to_csv('macaw_dead_end_and_duplicate_test_results.csv', index = False)

which would print this to the terminal:

Starting dead-end test...
 - Found 1370 dead-end metabolites.
 - Found 1123 reactions incapable of sustaining steady-state fluxes in either direction due to these dead-ends.
 - Found 1594 reversible reactions that can only carry steady-state fluxes in a single direction due to dead-ends.
Starting duplicate test...
 - Found 706 reactions that were some type of duplicate:
   - 0 were completely identical to at least one other reaction.
   - 13 involve the same metabolites but go in the opposite direction or have the opposite reversibility as at least one other reaction.
   - 445 involve the same metabolites but with different coefficients as at least one other reaction.
   - 261 redox reactions oxidize and/or reduce the same metabolites as another redox reaction but use a different pair of the given redox carriers.

and make macaw_dead_end_and_duplicate_test_results.csv look like:

reaction_id	reaction_equation	dead_end_test	duplicate_test_exact	duplicate_test_directions	duplicate_test_coefficients	duplicate_test_redox
MAR03905	MAM01796c + MAM02552c --> MAM01249c + MAM02039c + MAM02553c	ok	ok	ok	ok	MAR03907
MAR03907	MAM01796c + MAM02554c --> MAM01249c + MAM02039c + MAM02555c	ok	ok	ok	ok	MAR03905
⋮	⋮	⋮	⋮	⋮	⋮	⋮