micom-dev / micom

Python package to study microbial communities using metabolic modeling.
https://micom-dev.github.io/micom
Apache License 2.0
82 stars 17 forks source link

optimize_all(fluxes=True) and optimize_single(fluxes=True) do not return fluxes #33

Closed michaelsilverstein closed 3 years ago

michaelsilverstein commented 3 years ago

Even with fluxes=True passed into community.optimize_all(), only maximal growth rates are returned. In fact, community.optimize_single() does not accept the fluxes argument.

Using version '0.22.6'

cdiener commented 3 years ago

Yeah, optimize_all shouldn't take the argument either, that is an error in the signature and docs. In most cases you will want to use community.optimize or community.cooperative_tradeoff. See https://micom-dev.github.io/micom/growth_fluxes.html. What were you trying to achieve?

michaelsilverstein commented 3 years ago

Thanks for the response, that makes sense!

I'm comparing carbon use efficiency (CUE) of single organisms and communities of those organisms. The calculation for CUE requires the rate of some exchange fluxes in addition the biomass flux.

My plan was to construct communities and use optimize_single() to get CUE for each community member independently and then compare with community.optimize(). I'm still able to get the monoculture CUE by loading the models into cobrapy separately, but that's obviously less efficient since the models need to be loaded multiple times.

michaelsilverstein commented 3 years ago

The main thing I'm trying to do is create subsets of a community without doing multiple model loads. A hacky solution I'm going to try and move forward with is to just update the abundance profile to reflect the subset I'm interested in.

michaelsilverstein commented 3 years ago

I'm having success with first setting community._rtol = 0 and then community.set_abundance()

cdiener commented 3 years ago

Reopening because I still need to fix the signature. Updating the abundance is good but may renormalize unless you use com.set_abundances(..., normalize=False). What also worak is to just knockout all reactions for some taxa, which will work with the context manager (this will require a fix from 0.22.7 which I'll release today).

from micom import Community
from micom.data import test_taxonomy

tax = test_taxonomy()
com = Community(tax)
knocked_taxa = com.taxa[0:2]
with com:
    [r.knock_out() for r in com.reactions.query(lambda ri: ri.community_id in knocked_taxa)]
    sol_knocked = com.cooperative_tradeoff()

# automatically reset on leaving the with block
sol_normal = com.cooperative_tradeoff()

print(sol_normal.members)
print(sol_knocked.members)

This is also a little bit more numerically stable and is the strategy used in knockout_taxa.

cdiener commented 3 years ago

Okay, fixed now. Feel free to open a discussion at https://github.com/micom-dev/micom/discussions for a strategy to implement partial models if you have additional questions.