Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor: transform production envelope to wide format with multiple carbon sources #587

Merged
merged 7 commits into from
Sep 13, 2017

Conversation

Midnighter
Copy link
Member

@Midnighter Midnighter commented Aug 29, 2017

fix #582

So far this PR re-organizes the code for production_envelope a little bit and catches the potential ZeroDivisionError in mass_yield. However, I suggest the following further changes which is why I marked this PR as WIP.

  1. Allow specification of more than one carbon source.
  2. Transform the resulting data frame into wide format (the way it's done in cameo) meaning that maximum and minimum are separate columns.
  3. Use the configured solver threshold to convert below threshold values to zero.

@Midnighter Midnighter added the WIP work in progress label Aug 29, 2017
@Midnighter Midnighter self-assigned this Aug 29, 2017
@hredestig
Copy link
Contributor

On 2, being used to ggplot's interface, I find long format much more convenient as the column names remain fixed and then easier to write general plotting code. But, I guess that may be somewhat niche case and that wide format feels more familiar.

@Midnighter
Copy link
Member Author

Midnighter commented Aug 31, 2017

Actually, even for ggplot2 I find the wide variant more convenient because I opted to use geom_ribbon.

ggplot2::ggplot(total_df, ggplot2::aes(x = biomass,
                                       ymin = flux_minimum,
                                       ymax = flux_maximum,
                                       color = label,
                                       fill = label)) +
  ggplot2::geom_ribbon(alpha=0.3) +
  ggplot2::scale_color_brewer(name = "Model", palette = "Set1") +
  ggplot2::scale_fill_brewer(name = "Model", palette = "Set1") +
  ggplot2::labs(x = expression(paste("Growth [", g[DW], h^{-1}, "]")),
                y = expression(paste("Flux [", mmol, g[DW]^{-1}, h^{-1}, "]")))

Which can produce something like the following plot, for example.

production envelopes

@hredestig
Copy link
Contributor

Sure that works as well, I find it easiest to know the right format when having tried to implement many different plots which I didn't do yet (e.g. to recycle code with plotting e.g. fva, however that should be done the best way.). Something to think of when implementing overall better plot functions for cobra. In the mean time, switch to wide format here totally fine with me.

@Midnighter
Copy link
Member Author

Midnighter commented Sep 6, 2017

  1. Done
  2. Done
  3. TBD

Any other ideas that should be included here? I think testing could be expanded for sure.

@Midnighter Midnighter changed the title refactor: catch zero division in mass yield [WIP] refactor: transform production envelope to wide format with multiple carbon sources Sep 6, 2017
@codecov-io
Copy link

codecov-io commented Sep 6, 2017

Codecov Report

Merging #587 into devel will increase coverage by 0.01%.
The diff coverage is 82.05%.

Impacted file tree graph

@@            Coverage Diff             @@
##            devel     #587      +/-   ##
==========================================
+ Coverage   71.43%   71.45%   +0.01%     
==========================================
  Files          65       65              
  Lines        8778     8784       +6     
  Branches     1483     1491       +8     
==========================================
+ Hits         6271     6277       +6     
  Misses       2239     2239              
  Partials      268      268
Impacted Files Coverage Δ
cobra/test/test_flux_analysis.py 86.33% <100%> (+0.08%) ⬆️
cobra/core/model.py 90.45% <66.66%> (+0.02%) ⬆️
cobra/flux_analysis/phenotype_phase_plane.py 68.12% <80.88%> (+0.25%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ea866f3...69efe6a. Read the comment docs.

results['carbon_source'] = carbon_io[0].id
return results
with model:
model.solver.objective.direction = sense
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't use context manager, perhaps better formulate as model.objective = .. ? Otherwise I guess will exit with objective changed to max if starting with min..

for rxn in reactions:
point = grid.at[i, rxn.id]
rxn.bounds = point, point
model.slim_optimize()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

perhaps obj_value = model.slim_optimize() and then use obj_value below

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed.

model.slim_optimize()
if model.solver.status != OPTIMAL:
continue
if abs(model.solver.objective.value) < \
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

how about making the tolerance an optional argument and have thereby have this clipping documented?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea.

try:
return carbon_output_flux / carbon_input_flux
except ZeroDivisionError:
return nan


def mass_yield(c_input_output):
"""Gram product divided by gram of carbon input source
def reaction_components(reaction, atom='C'):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

element not atom.. And since you add this, perhaps do it all the way and allow user to calculate any element yield by having e.g. yield_element='N' as argument to production_envelope? Or do later and remove parametrisation here..

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removing parametrization here for now. It's a little annoying with column naming so unless there is a need for it I'd skip it.

@Midnighter
Copy link
Member Author

This PR is now based on #588 and the new Model.objective_direction property.

@Midnighter Midnighter force-pushed the fix/production_envelope branch 3 times, most recently from 01329b2 to f292c28 Compare September 12, 2017 14:34
grid.at[i, 'carbon_yield_{}'.format(direction)] = \
total_yield([rxn.flux for rxn in c_input],
input_components,
model.solver.objective.value,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

obj_val ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Absolutely!

@Midnighter Midnighter merged commit 996ba7b into devel Sep 13, 2017
@Midnighter Midnighter deleted the fix/production_envelope branch September 13, 2017 09:28
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
WIP work in progress
Projects
None yet
Development

Successfully merging this pull request may close these issues.

ZeroDivisionError in production_envelope
3 participants