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

exchanges not identified correctly in CarveME models #9

Closed
arianccbasile opened this issue Sep 5, 2019 · 14 comments
Closed

exchanges not identified correctly in CarveME models #9

arianccbasile opened this issue Sep 5, 2019 · 14 comments

Comments

@arianccbasile
Copy link

Goodmorning.
I have recently used your program to analyze metabolic exchanges in a community. I have only a question: in the output, I have for each metabolite two columns one called, for example, EX_but_e and the other EX_but_m. The EX_but_m appears only in one row which is the medium one. I supposed that summing the EX_but_e for all the microbes in my community I should have got the value in the medium row at EX_but_m column. But it is not the case. Why?

Sincerely,
Arianna Basile

@cdiener
Copy link
Collaborator

cdiener commented Sep 5, 2019

Hi,

fluxes for individual taxa and exchange fluxes for the external medium have different units in micom. For an individual taxon the unit is mmol/[gDW taxon * h] but for the external medium the unit is mmol/[gDW total biomass * h]. This is done so that fluxes for individual taxa are comparable. So in order to get that to line up you have to scale all individual taxon fluxes by their abundance. So in your case: EX_but_m = Σ a_i * EX_but_e_i where i denotes individual taxa and a_i their relative abundances.

@arianccbasile
Copy link
Author

Hello,
thank you for the kind answer. Now everything is clearer.

Sincerely,
Arianna Basile

@cdiener
Copy link
Collaborator

cdiener commented Sep 6, 2019

Cool, happy that helped. Will close the issue for now. Feel free to reopen if there are more issues.

@cdiener cdiener closed this as completed Sep 6, 2019
@arianccbasile
Copy link
Author

Good morning,
I have another question: the fluxes with negative sign should refer to metabolites consumed while those with positive sign should refer to metabolites produced. Is it right? I don't understand why for some metabolites I have more absorption than production on the whole of the community. How is it possible?

Sincerely,
Arianna Basile

@cdiener
Copy link
Collaborator

cdiener commented Sep 10, 2019

For the AGORA models exchange fluxes are usually directed towards import, ergo ∅ --> M_c, so positive fluxes mean imports and negative fluxes mean exports. Could you expand on what you mean that you have more adsorption then production for a single metabolite? For a single metabolite there should only be net adsorption or net production. Usually you can't have both at the same time. Or if you have an example of specific output and want to mark what seems wrong that would help too. Thanks!

@cdiener cdiener reopened this Sep 10, 2019
@arianccbasile
Copy link
Author

arianccbasile commented Sep 12, 2019

Thank you very much for previous support!

I am using models created with carveme and I am thinking it works the other way around because in my simulation archaea do produce methane (flux >0).

For the second point I fear I wasn't clear enough. This is a column of my output file.
L-Histidine his__L 0 -32.5477 0 0 -15.626 -17.6141 0 0 -7.59349 -0.4545 0 0 -0.0696 0 -0.59577 -0.09332 0 0 0 0 0 0 0 -0.22369 0 -0.07805 -0.11424 0 0 0 -0.07728 0 -0.06839 -0.15902 -0.17386 -0.0789 -0.30852 0 0 -0.03839 -0.06967 -0.05813 0 -0.0834 -0.00579 0 0 -0.05972 0 -0.12363 0 -0.06973 0 0 -0.05215 -0.03324 -0.01289 0 0 0 -0.01275 -0.00994 0 -0.02471 0 -0.04649 -0.03175 -0.02737 0
I have some import and some exports because I have a complex community but the production is different from the consumed and sometimes, like in this example, the total metabolite consumption is different than 0 but the production is equal to 0. So where that product come from?

Sincerely,
Arianna

@cdiener
Copy link
Collaborator

cdiener commented Sep 12, 2019

Hi Anna,

you are absolutely right, I confused the direction. In AGORA as well as CarveMe exhanges are directed towards export so negative flux means import and positive export. I was confusing it with the output of minimal_medium which will turn that around.

For the second part I am still not entirely sure if I get what you mean, so let me try to understand better. For total metabolite consumption do you mean the overall import flux (sum of all import fluxes for the specific taxon)? In that case it can definitely happen that the total import flux is larger than the total export flux since individual taxa produce biomass which consumes metabolites. There may also be sink reactions consuming internal metabolites to fulfill a maintenance requirement. For instance sometimes one knows the total ATP need of a bacterium and may insert an additional reaction that consumes/destroys ATP to compensate for missing pathways in the model. The opposite is generally not possible (total export > total import) unless you have source reactions in your model (reactions that create mass from nothing, usually not the case).

@arianccbasile
Copy link
Author

Thank you for clarification for the first part of the question.
For what it regards the second part I fear I was not clear enough. I mean for each compound in the community I have some microbes as takers and some of them as givers but sometimes I have that the takers as a whole take more than what the givers produce. How is it possible?

Sincerely,
Arianna

@cdiener
Copy link
Collaborator

cdiener commented Sep 19, 2019

Again this would have to be scaled by the abundances first since those would denote exchange fluxes. Also the community can take up metabolites from the environment. This is dictated by the media exchanges with names like EX_*_m. So your metabolite may also be taken up from the outside additionally to what is produced by the taxa. Also see https://micom-dev.github.io/micom/media.html#Applying-growth-media on how to modify those global exchanges.

@arianccbasile
Copy link
Author

arianccbasile commented Jan 20, 2020


I am trying to apply my own growth medium as you suggested but I am having some issues. I tried many media files but nothing seems to work properly.
I have a tsv file with only one row with
EX_glc__D_m\t10.00
I am reading the file with:

acetate_medium=pandas.Series.from_csv("medium.tsv",sep="\t")

to import it as a pandas.core.series.Series and then I am trying to use it as my medium with

com.medium = acetate_medium

But it seems not to work properly since I get the following error


KeyError Traceback (most recent call last)
in
8 taxonomy=pandas.read_csv(in_file)
9 com = Community(taxonomy, solver="gurobi")
---> 10 com.medium = acetate_medium
11 cherrypy.log("Build a community with a total of {} reactions.".format(len(com.reactions)))
12 sol=com.optimize(fluxes=True,pfba=True)

/mnt/data_SSD/anaconda3/envs/flux_balance/lib/python3.6/site-packages/cobra/core/model.py in medium(self, medium)
283 exchange_rxns = frozenset(self.exchanges)
284 for rxn_id, bound in iteritems(medium):
--> 285 rxn = self.reactions.get_by_id(rxn_id)
286 if rxn not in exchange_rxns:
287 logger.warn("%s does not seem to be an"

/mnt/data_SSD/anaconda3/envs/flux_balance/lib/python3.6/site-packages/cobra/core/dictlist.py in get_by_id(self, id)
56 def get_by_id(self, id):
57 """return the element with a matching id"""
---> 58 return list.getitem(self, self._dict[id])
59
60 def list_attr(self, attribute):

KeyError: 'EX_glc__D_m'

Thank you in advance for your help.
Sincerely,
Arianna Basile

@cdiener
Copy link
Collaborator

cdiener commented Jan 20, 2020

Hi Arianna,

the names of the exchange reactions depend a lot on the original models you put in the community construction and need to match the global import names of your model. You can see a list of all global imports by using Community.exchanges. For instance:

In [1]: import micom

In [2]: tax = micom.data.test_taxonomy()

In [3]: com = micom.Community(tax)
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:02<00:00,  2.43models/s]

In [4]: com.exchanges
Out[4]:
[<Reaction EX_ac_m at 0x7f9da9d5be48>,
 <Reaction EX_acald_m at 0x7f9da9d5b710>,
 <Reaction EX_akg_m at 0x7f9da9d6c390>,
 <Reaction EX_co2_m at 0x7f9da9d6c588>,
 <Reaction EX_etoh_m at 0x7f9da9d6c748>,
 <Reaction EX_for_m at 0x7f9da9d6c908>,
 <Reaction EX_fru_m at 0x7f9da9d6cac8>,
 <Reaction EX_fum_m at 0x7f9da9d6cc88>,
 <Reaction EX_glc__D_m at 0x7f9da9d6cdd8>,
 <Reaction EX_gln__L_m at 0x7f9da9d6cf60>,
 <Reaction EX_glu__L_m at 0x7f9da9d80128>,
 <Reaction EX_h_m at 0x7f9da9d80320>,
 <Reaction EX_h2o_m at 0x7f9da9d80518>,
 <Reaction EX_lac__D_m at 0x7f9da9d80668>,
 <Reaction EX_mal__L_m at 0x7f9da9d807f0>,
 <Reaction EX_nh4_m at 0x7f9da9d809b0>,
 <Reaction EX_o2_m at 0x7f9da9d80b70>,
 <Reaction EX_pi_m at 0x7f9da9d80d68>,
 <Reaction EX_pyr_m at 0x7f9da9d80f98>,
 <Reaction EX_succ_m at 0x7f9da9d87198>]

The reaction names in your medium have to correspond exactly to those IDs you see there. One common pitfall is that the model may no have import for some of the metabolites. So I usually do something like that in my workflows (here diet is the growth medium).

ex_ids = [r.id for r in com.exchanges]
logger.info(
    "%d/%d import reactions found in model.",
    diet.index.isin(ex_ids).sum(),
    len(diet)
)
com.medium = diet[diet.index.isin(ex_ids)]

If the log tells you that no imports were found something is going awry.

I haven't really used CarveME reconstructions so there might be some issues there with reaction IDs. But I would be interested to support those more and would be happy to help if you send me one of the community reconstructions to investigate. You can do so in private by saving the community model with com.to_pickle("model.pickle") and sharing this file with with Google Drive or something similar via E-mail to cdiener [at] isbscience.org (or attaching directly if it is small enough).

@arianccbasile
Copy link
Author

arianccbasile commented Jan 21, 2020

Thank you for your fast and accurate answer. I tried what you proposed and

com.exchanges

gives me an empty dictionary so yes I think there is a problem. I gave a glance also to

com.reactions

and I spotted some reactions with the ID starting with "EX_" flag so I am trying changing the bounds to those reactions.
I am also sending you the pickle file of the community model.

Sincerely,
Arianna Basile

@cdiener cdiener changed the title Issues on output exchanges not identified correctly in CarveME models Jan 21, 2020
@cdiener
Copy link
Collaborator

cdiener commented Jan 21, 2020

I played around with some CarveME models and MICOM does indeed not identify exchange reactions correctly there. I added a fix and it now seems to work:

In [1]: import micom as mm

In [2]: path = "/Users/cdiener/Downloads/Bacteroides_vulgatus_ATCC_8482.xml.gz"

In [3]: import pandas as pd

In [4]: tax = pd.DataFrame({"id": ["bac1", "bac2"], "file": [path, path]})

In [5]: com = mm.Community(tax)
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2/2 [00:07<00:00,  3.90s/models]

In [6]: com.exchanges
Out[6]:
[<Reaction EX_glc__D_e_m at 0x7fdc6226e390>,
 <Reaction EX_h2o_e_m at 0x7fdc702cb358>,
 <Reaction EX_23camp_e_m at 0x7fdc6226e6d8>,
 <Reaction EX_23ccmp_e_m at 0x7fdc6226e080>,
 <Reaction EX_23cgmp_e_m at 0x7fdc6226e550>,
 <Reaction EX_23cump_e_m at 0x7fdc6226e208>,
 <Reaction EX_h_e_m at 0x7fdc6226e240>,
 <Reaction EX_2pglyc_e_m at 0x7fdc6226e6a0>,
 <Reaction EX_leu__L_e_m at 0x7fdc62263be0>,
 <Reaction EX_3hpp_e_m at 0x7fdc62263940>,
 <Reaction EX_ala__L_e_m at 0x7fdc62263f28>,
 <Reaction EX_cl_e_m at 0x7fdc62263a20>,
 <Reaction EX_arab__L_e_m at 0x7fdc62263c88>,
 <Reaction EX_4abut_e_m at 0x7fdc622630f0>,
 <Reaction EX_acald_e_m at 0x7fdc622632b0>,
 <Reaction EX_ac_e_m at 0x7fdc62263160>,
 <Reaction EX_chol_e_m at 0x7fdc62256c88>,
 <Reaction EX_pi_e_m at 0x7fdc62256828>,
 <Reaction EX_adn_e_m at 0x7fdc62256f28>,
 <Reaction EX_ins_e_m at 0x7fdc62256470>,
 <Reaction EX_nh4_e_m at 0x7fdc62256198>,
 <Reaction EX_ade_e_m at 0x7fdc62256208>,
 <Reaction EX_mal__L_e_m at 0x7fdc62256128>,
 <Reaction EX_alaala_e_m at 0x7fdc6224bf98>,
 <Reaction EX_asn__L_e_m at 0x7fdc62256358>,
 <Reaction EX_cys__L_e_m at 0x7fdc6224bfd0>,
 <Reaction EX_ala__D_e_m at 0x7fdc6224b8d0>,
 <Reaction EX_gln__L_e_m at 0x7fdc6224bac8>,
 <Reaction EX_arg__L_e_m at 0x7fdc6224b208>,
 <Reaction EX_cgly_e_m at 0x7fdc6224b9e8>,
 <Reaction EX_malt_e_m at 0x7fdc6224b978>,
 <Reaction EX_malttr_e_m at 0x7fdc6224b2e8>,
 <Reaction EX_starch1200_e_m at 0x7fdc6224b470>,
 <Reaction EX_amylose300_e_m at 0x7fdc6224b240>,
 <Reaction EX_arab__D_e_m at 0x7fdc622406d8>,
 <Reaction EX_fe3_e_m at 0x7fdc62240d68>,
 <Reaction EX_asp__L_e_m at 0x7fdc62240e80>,
 <Reaction EX_aso3_e_m at 0x7fdc622407f0>,
 <Reaction EX_aso4_e_m at 0x7fdc62240048>,
 <Reaction EX_k_e_m at 0x7fdc622408d0>,
 <Reaction EX_pro__L_e_m at 0x7fdc62240710>,
 <Reaction EX_meoh_e_m at 0x7fdc622407b8>,
 <Reaction EX_buts_e_m at 0x7fdc62240358>,
 <Reaction EX_but_e_m at 0x7fdc62240278>,
 <Reaction EX_ca2_e_m at 0x7fdc62235be0>,
 <Reaction EX_cell6_e_m at 0x7fdc622356a0>,
 <Reaction EX_chols_e_m at 0x7fdc62235160>,
 <Reaction EX_cit_e_m at 0x7fdc62235ac8>,
 <Reaction EX_mg2_e_m at 0x7fdc622355c0>,
 <Reaction EX_mn2_e_m at 0x7fdc6222add8>,
 <Reaction EX_cobalt2_e_m at 0x7fdc6222a7f0>,
 <Reaction EX_zn2_e_m at 0x7fdc6222a668>,
 <Reaction EX_succ_e_m at 0x7fdc6222af28>,
 <Reaction EX_for_e_m at 0x7fdc6222a550>,
 <Reaction EX_co2_e_m at 0x7fdc6222ac18>,
 <Reaction EX_co_e_m at 0x7fdc6222a2e8>,
 <Reaction EX_csn_e_m at 0x7fdc6222a048>,
 <Reaction EX_cu2_e_m at 0x7fdc62225eb8>,
 <Reaction EX_o2_e_m at 0x7fdc62225978>,
 <Reaction EX_glu__L_e_m at 0x7fdc622254a8>,
 <Reaction EX_cytd_e_m at 0x7fdc62225550>,
 <Reaction EX_dad_2_e_m at 0x7fdc622258d0>,
 <Reaction EX_din_e_m at 0x7fdc6221e2b0>,
 <Reaction EX_dca_e_m at 0x7fdc6221eeb8>,
 <Reaction EX_dcyt_e_m at 0x7fdc6221e5c0>,
 <Reaction EX_dgsn_e_m at 0x7fdc6221a550>,
 <Reaction EX_dha_e_m at 0x7fdc6221a160>,
 <Reaction EX_lac__D_e_m at 0x7fdc6221acf8>,
 <Reaction EX_drib_e_m at 0x7fdc6221a1d0>,
 <Reaction EX_duri_e_m at 0x7fdc6221a198>,
 <Reaction EX_enlipa_e_m at 0x7fdc6221a978>,
 <Reaction EX_galur_e_m at 0x7fdc62213908>,
 <Reaction EX_ttdca_e_m at 0x7fdc62213fd0>,
 <Reaction EX_etha_e_m at 0x7fdc62213198>,
 <Reaction EX_etoh_e_m at 0x7fdc6220ee10>,
 <Reaction EX_h2_e_m at 0x7fdc6220e550>,
 <Reaction EX_fe2_e_m at 0x7fdc6220eb38>,
 <Reaction EX_enter_e_m at 0x7fdc62213b38>,
 <Reaction EX_feenter_e_m at 0x7fdc62208a58>,
 <Reaction EX_fol_e_m at 0x7fdc62208828>,
 <Reaction EX_fuc__L_e_m at 0x7fdc62208400>,
 <Reaction EX_fum_e_m at 0x7fdc62201e10>,
 <Reaction EX_g3pe_e_m at 0x7fdc622019b0>,
 <Reaction EX_gal_e_m at 0x7fdc62201438>,
 <Reaction EX_gal_bD_e_m at 0x7fdc62201128>,
 <Reaction EX_galt_e_m at 0x7fdc621facc0>,
 <Reaction EX_gthrd_e_m at 0x7fdc621fa390>,
 <Reaction EX_glcur_e_m at 0x7fdc621fac18>,
 <Reaction EX_glucan6_e_m at 0x7fdc621facf8>,
 <Reaction EX_glucan4_e_m at 0x7fdc621fa208>,
 <Reaction EX_glyclt_e_m at 0x7fdc621faa90>,
 <Reaction EX_glycogen1500_e_m at 0x7fdc621f4470>,
 <Reaction EX_glyc_e_m at 0x7fdc621f4f98>,
 <Reaction EX_gsn_e_m at 0x7fdc621f4080>,
 <Reaction EX_gthox_e_m at 0x7fdc621fa2b0>,
 <Reaction EX_h2o2_e_m at 0x7fdc621ef9e8>,
 <Reaction EX_gua_e_m at 0x7fdc621ef438>,
 <Reaction EX_h2s_e_m at 0x7fdc621ef630>,
 <Reaction EX_3hcinnm_e_m at 0x7fdc621efef0>,
 <Reaction EX_4hphac_e_m at 0x7fdc621e8198>,
 <Reaction EX_indole_e_m at 0x7fdc621e8ba8>,
 <Reaction EX_inost_e_m at 0x7fdc621e81d0>,
 <Reaction EX_kdo2lipid4_e_m at 0x7fdc621e8ac8>,
 <Reaction EX_Larab_e_m at 0x7fdc621e2c50>,
 <Reaction EX_lcts_e_m at 0x7fdc621e2e10>,
 <Reaction EX_LalaLglu_e_m at 0x7fdc621e2470>,
 <Reaction EX_lac__L_e_m at 0x7fdc621e8780>,
 <Reaction EX_lmn2_e_m at 0x7fdc621ddb00>,
 <Reaction EX_lmn30_e_m at 0x7fdc621dd8d0>,
 <Reaction EX_lyx__L_e_m at 0x7fdc621dd6a0>,
 <Reaction EX_manttr_e_m at 0x7fdc621e2898>,
 <Reaction EX_melib_e_m at 0x7fdc621d7e10>,
 <Reaction EX_n2o_e_m at 0x7fdc621d7f98>,
 <Reaction EX_no2_e_m at 0x7fdc621d7550>,
 <Reaction EX_no3_e_m at 0x7fdc621d2f98>,
 <Reaction EX_no_e_m at 0x7fdc621d2630>,
 <Reaction EX_uri_e_m at 0x7fdc621d20f0>,
 <Reaction EX_udcpo5_e_m at 0x7fdc621d29b0>,
 <Reaction EX_octa_e_m at 0x7fdc621cda20>,
 <Reaction EX_pac_e_m at 0x7fdc621cd630>,
 <Reaction EX_pdima_e_m at 0x7fdc621cd588>,
 <Reaction EX_pnto__R_e_m at 0x7fdc621cdda0>,
 <Reaction EX_ppa_e_m at 0x7fdc621c6b00>,
 <Reaction EX_progly_e_m at 0x7fdc621c6eb8>,
 <Reaction EX_pullulan1200_e_m at 0x7fdc621c64e0>,
 <Reaction EX_pyr_e_m at 0x7fdc621c6198>,
 <Reaction EX_so4_e_m at 0x7fdc621c6908>,
 <Reaction EX_val__L_e_m at 0x7fdc621c07b8>,
 <Reaction EX_raffin_e_m at 0x7fdc621c0240>,
 <Reaction EX_rbt_e_m at 0x7fdc621c04a8>,
 <Reaction EX_spmd_e_m at 0x7fdc621c0ac8>,
 <Reaction EX_rib__D_e_m at 0x7fdc621b8128>,
 <Reaction EX_rmn_e_m at 0x7fdc621b8fd0>,
 <Reaction EX_sbt__D_e_m at 0x7fdc621b8518>,
 <Reaction EX_s_e_m at 0x7fdc621b1ef0>,
 <Reaction EX_tartr__D_e_m at 0x7fdc621b8908>,
 <Reaction EX_thymd_e_m at 0x7fdc621b1780>,
 <Reaction EX_thym_e_m at 0x7fdc621b1cc0>,
 <Reaction EX_xan_e_m at 0x7fdc621b11d0>,
 <Reaction EX_xtsn_e_m at 0x7fdc621b1a90>,
 <Reaction EX_xylan12_e_m at 0x7fdc621aba90>,
 <Reaction EX_xylan4_e_m at 0x7fdc621abf98>,
 <Reaction EX_xylan8_e_m at 0x7fdc621ab6a0>,
 <Reaction EX_xyl3_e_m at 0x7fdc621abdd8>,
 <Reaction EX_xyl__D_e_m at 0x7fdc621a69b0>,
 <Reaction EX_xylb_e_m at 0x7fdc621a6240>]

(before com.exchanges would return []). This was due to CarveME using a different naming scheme for compartments so for instance the compartment for EX_glc__D_e is named C_e in CarveMe and not e as the name would imply. MICOM now uses the new cobrapy heuristics for identifying boundary types so that should be more robust. The naming is a bit weird due to the _e_m part but I would not know how to make this better if compartment indicators are not consistent for reactions in CarveMe.

This fix will be released in version 0.10.1 of MICOM. So an update via pip should be enough.

@arianccbasile
Copy link
Author

The weird name doesn't seem a problem to me. Thank you very much, Christian. :)

Sincerely,
Arianna Basile

@cdiener cdiener closed this as completed Feb 5, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants