Similar to the western-style diet we will again start by loading the diet and depleting components absorbed by the host. In this case we have no manual annotation for which components should be diluted so we will use a generic human metabolic model to find those. The growth medium supllied here was created the following way:
Let's start by reading the diet which was downloaded from https://www.vmh.life/#diet/High%20fat,%20low%20carb. Flux is in mmol/human/day. This has to be adjusted to 1 hour. Also the VMH site has a bug where it will clip fluxes after 4 digits, so we will set values like 0.0000 to 0.0001.
import pandas as pd
medium = pd.read_csv("../data/vmh_high_fat_low_carb.tsv", index_col=False, sep="\t")
medium.columns = ["reaction", "flux"]
medium.reaction = medium.reaction.str.replace("(\[e\]$)|(\(e\)$)", "", regex=True)
medium.loc[medium.flux < 1e-4, "flux"] = 1e-4
medium.flux = medium.flux / 24
medium
reaction | flux | |
---|---|---|
0 | EX_etoh | 0.000004 |
1 | EX_h2o | 6608.413339 |
2 | EX_caro | 0.000068 |
3 | EX_retinol | 0.129458 |
4 | EX_thm | 0.075371 |
... | ... | ... |
86 | EX_lgnc | 0.008570 |
87 | EX_fol | 0.000005 |
88 | EX_strch1 | 0.002862 |
89 | EX_i | 0.000320 |
90 | EX_starch1200 | 0.000026 |
91 rows × 2 columns
Now we will try to identify components that can be taken up by human cells.
To achieve this we will load the Recon3 human model. AGORA and Recon IDs are very similar so we should be able to match them. We just have to adjust the Recon3 ones a bit. We start by identifying all available exchanges in Recon3 and adjusting the IDs.
from cobra.io import read_sbml_model
import pandas as pd
recon3 = read_sbml_model("../data/Recon3D.xml.gz")
exchanges = pd.Series([r.id for r in recon3.exchanges])
exchanges = exchanges.str.replace("__", "_").str.replace("_e$", "", regex=True)
exchanges.head()
0 EX_5adtststerone 1 EX_5adtststerones 2 EX_5fthf 3 EX_5htrp 4 EX_5mthf dtype: object
Now we will check which ones we can find in our set and add in the dilution factors (again going with 1:10.
medium["dilution"] = 1.0
medium.loc[medium.reaction.isin(exchanges), "dilution"] = 0.1
medium.dilution.value_counts()
0.1 79 1.0 12 Name: dilution, dtype: int64
Okay, so 79/91 components can be adsorbed by humans. We end by filling in the additional info.
medium["metabolite"] = medium.reaction.str.replace("^EX_", "", regex=True) + "_m"
medium["global_id"] = medium.reaction + "(e)"
medium["reaction"] = medium.reaction + "_m"
medium.loc[medium.flux < 1e-4, "flux"] = 1e-4
medium
reaction | flux | dilution | metabolite | global_id | |
---|---|---|---|---|---|
0 | EX_etoh_m | 0.000100 | 0.1 | etoh_m | EX_etoh(e) |
1 | EX_h2o_m | 6608.413339 | 0.1 | h2o_m | EX_h2o(e) |
2 | EX_caro_m | 0.000100 | 0.1 | caro_m | EX_caro(e) |
3 | EX_retinol_m | 0.129458 | 0.1 | retinol_m | EX_retinol(e) |
4 | EX_thm_m | 0.075371 | 0.1 | thm_m | EX_thm(e) |
... | ... | ... | ... | ... | ... |
86 | EX_lgnc_m | 0.008570 | 0.1 | lgnc_m | EX_lgnc(e) |
87 | EX_fol_m | 0.000100 | 0.1 | fol_m | EX_fol(e) |
88 | EX_strch1_m | 0.002862 | 0.1 | strch1_m | EX_strch1(e) |
89 | EX_i_m | 0.000320 | 0.1 | i_m | EX_i(e) |
90 | EX_starch1200_m | 0.000100 | 1.0 | starch1200_m | EX_starch1200(e) |
91 rows × 5 columns
But can the bacteria in our model database actually grow on this medium? Let's check and start by downbloading the AGORA model database.
# !wget https://zenodo.org/record/3755182/files/agora103_genus.qza?download=1 -O data/agora103_genus.qza
No we we will check for growth by running the growth medium against any single model.
from micom.workflows.db_media import check_db_medium
check = check_db_medium("../data/agora103_genus.qza", medium, threads=20)
Output()
check
now includes the entire manifest plus two new columns: the growth rate and whether the models can grow.
check.can_grow.value_counts()
False 227 Name: can_grow, dtype: int64
Okay nothing can grow. We probably miss some important cofactor such as manganese or copper.
Let's complete the medium so that all taxa in Refseq can grow at a rate of at least 1e-4.
Sometimes you may start from a few componenents and will want to complete this skeleton medium to reach a certain minimum growth rate across all models in the database. This can be done with complete_db_medium
. We can minimize either the added total flux, mass or presence of any atom. Since, we want to build a low carb diet here we will minimize the presence of added carbon.
from micom.workflows.db_media import complete_db_medium
manifest, imports = complete_db_medium("../data/agora103_genus.qza", medium, growth=0.001, threads=20, max_added_import=10, weights="C")
Output()
manifest.can_grow.value_counts()
True 227 Name: can_grow, dtype: int64
manifest
is the amended manifest as before and imports
contains the used import fluxes for each model. A new column in the manifest also tells us how many import were added.
manifest.added.describe()
count 227.000000 mean 6.678414 std 3.959711 min 1.000000 25% 3.000000 50% 7.000000 75% 9.000000 max 22.000000 Name: added, dtype: float64
So we added 7 metabolites on average (1-22).
From this we build up our new medium.
fluxes = imports.max()
fluxes = fluxes[(fluxes > 1e-6) | fluxes.index.isin(medium.reaction)]
completed = pd.DataFrame({
"reaction": fluxes.index,
"metabolite": fluxes.index.str.replace("^EX_", "", regex=True),
"global_id": fluxes.index.str.replace("_m$", "(e)", regex=True),
"flux": fluxes
})
completed.shape
(122, 4)
Let's also export the medium as Qiime 2 artifact which can be read with q2-micom
or the normal micom package.
from qiime2 import Artifact
arti = Artifact.import_data("MicomMedium[Global]", completed)
arti.save("../media/vmh_high_fat_low_carb_agora.qza")
'../media/vmh_high_fat_low_carb_agora.qza'
As a last step we validate the created medium.
check = check_db_medium("../data/agora103_genus.qza", completed, threads=20)
check.can_grow.value_counts()
Output()
True 227 Name: can_grow, dtype: int64
check.growth_rate.describe()
count 227.000000 mean 0.002235 std 0.001484 min 0.001000 25% 0.001000 50% 0.002030 75% 0.002564 max 0.006467 Name: growth_rate, dtype: float64