Growth media VMH high fat low carb diet

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.

In [1]:
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
Out[1]:
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.

Identifying human adsorption

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.

In [2]:
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()
Out[2]:
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.

In [3]:
medium["dilution"] = 1.0
medium.loc[medium.reaction.isin(exchanges), "dilution"] = 0.1
medium.dilution.value_counts()
Out[3]:
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.

In [4]:
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
Out[4]:
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

Checking the growth medium against the DB

But can the bacteria in our model database actually grow on this medium? Let's check and start by downbloading the AGORA model database.

In [5]:
# !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.

In [6]:
from micom.workflows.db_media import check_db_medium

check = check_db_medium("../data/agora103_genus.qza", medium, threads=20)

check now includes the entire manifest plus two new columns: the growth rate and whether the models can grow.

In [7]:
check.can_grow.value_counts()
Out[7]:
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.

Supplementing a growth medium from a skeleton

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.

In [8]:
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")
In [9]:
manifest.can_grow.value_counts()
Out[9]:
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.

In [10]:
manifest.added.describe()
Out[10]:
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.

In [11]:
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
Out[11]:
(122, 4)

Let's also export the medium as Qiime 2 artifact which can be read with q2-micom or the normal micom package.

In [12]:
from qiime2 import Artifact

arti = Artifact.import_data("MicomMedium[Global]", completed)
arti.save("../media/vmh_high_fat_low_carb_agora.qza")
Out[12]:
'../media/vmh_high_fat_low_carb_agora.qza'

Validation

As a last step we validate the created medium.

In [13]:
check = check_db_medium("../data/agora103_genus.qza", completed, threads=20)
check.can_grow.value_counts()
Out[13]:
True    227
Name: can_grow, dtype: int64
In [14]:
check.growth_rate.describe()
Out[14]:
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