q2-micom provides a Qiime 2 plugin for the MICOM
Python package and allows for functional analysis of microbial communities using metabolic modeling.
Analysis of function in microbial communities is often performed based solely on sequencing data and reference genomes. This is a good way to gauge the metabolic potential of a microbiome, but the mere presence of a particular gene does not guarantee a particular metabolic output. For instance, the same E. coli cell can grow either aerobically or anaerobically. However, it will consume and produce very different sets of metabolites in these two environmental settings.
Metabolic modeling tries to estimate the activity of biochemical networks by predicting metabolic fluxes (i.e. the rates of mass conversions within a cell, usually expressed in mmol per gram dry-weight per hour). It uses the genotype and environmental conditions to establish stoichiometry and flux limits for all biochemical reactions within a particular cell and then tries to find sets of fluxes that yield maximum growth under a set of assumptions. This optimization of biomass production under a given stoichiometry is called Flux Balance Analysis (FBA).
MICOM extends the concept of FBA beyond a single genotype to entire microbial communities. In particular, it tries to reduce the space of feasible fluxes by imposing a tradeoff between optimization of community-wide biomass production and the individual (“selfish”) maximization of biomass production for each taxon. To do so, MICOM makes the following assumptions:
Not all of these assumptions may be fulfilled by the system you are studying. The farther you are from fulfilling them, the less accurate predictions by MICOM are likely to be.
MICOM models all biochemical reactions in all taxa, which means that the optimization problem MICOM solves included hundreds of thousands of variables. There are only a few numerical solvers that can solve quadratic programming problems of that scale. Right, now we support CPLEX or Gurobi, which both have free academic licenses but will require you to sign up for them. We hope to change this in the future by making MICOM compatible with an open-source solver (stay tuned!). If you have one of the two solvers mentioned above, you can install
q2-micom using the following steps.
You will need a Qiime 2 environment with version
2021.2 or higher (how to install Qiime 2). For instance on Linux you would use:
wget https://data.qiime2.org/distro/core/qiime2-2021.2-py36-linux-conda.yml conda env create -n qiime2-2021.2 --file qiime2-2021.2-py36-linux-conda.yml # OPTIONAL CLEANUP rm qiime2-2021.2-py36-linux-conda.yml
Once installed, we also download the environment for q2-micom and install it into the QIIME 2 environment we just created.
wget https://raw.githubusercontent.com/micom-dev/q2-micom/main/q2-micom.yml conda env update -n qiime2-2021.2 -f q2-micom.yml # OPTIONAL CLEANUP rm q2-micom.yml
Finally, you activate your environment.
conda activate qiime2-2021.2
q2-micom comes with a QP solver. You may install alternative solvers if for an additional speed-up.
See the README for instructions.
If you installed
q2-micom in an already existing Qiime 2 environment, update the plugin cache:
conda activate qiime2-2021.2 # or whatever you called your environment qiime dev refresh-cache
You are now ready to run
Below is a graphical overview of a
To start, we will need taxon abundances and a model database. Let’s look at model databases first.
To build metabolic community models, you will first need taxon-specific metabolic models for as many taxa in your sample as possible. Unfortunately, there is no magic involved here and we will need to provide this. Building high quality metabolic models is in art in and of itself. It can be done from from a genome alone (for instance using ModelSEED or CarveME) but usually requires additional curation. Furthermore, MICOM requires these models to be collapsed to the specific taxonomic rank you want to simulate (e.g. species or genus level). To facilitate this step we provide pre-built databases based on the AGORA model collection, which currently comprises manually curated metabolic models for 818 bacterial strains. For this tutorial, we will use models summarized at the genus level, as this is usually the lowest rank with decent annotation accuracy in 16S data.
You may find many more model databases at https://doi.org/10.5281/zenodo.3755182.
The first artifacts you will want to produce are the metabolic community models. These are sample-specific metabolic models that contain all taxa in a sample (i.e. all taxa represented in AGORA) at the correct relative abundance (i.e. as inferred from the 16S or metagenomic sequencing data). To build these models, you will need the relative abundance and taxonomic annotation for each sequence variant in your sample. Luckily, one of the main function of Qiime 2 is to obtain these features from your sequencing data…phew.
For the purposes of this tutorial, we will provide pre-built artifacts for a 16S amplicon sequencing data set from 8 healthy and 8 colorectal cancer stool samples taken from https://doi.org/10.1158/1940-6207.CAPR-14-0129. These artifacts were generated using DADA2 and the taxonomy was inferred using the SILVA database (version 132).
Using our pre-baked model database, we can now build our community models with the
qiime micom build command. Note that most commands in
q2-micom take a
--p-threads parameter that specifies how many CPU cores to use. Using multiple threads will speed things up considerably. Additionally, using the
--verbose flag will usually show a progress bar. Finally, low abundance taxa are usually dropped from the models to improve computational efficiency. This abundance cutoff is controlled by the
--p-cutoff parameter, which is set to 0.01% by default. Okay, let’s build our community models:
qiime micom build --i-abundance crc_table.qza \ --i-taxonomy crc_taxa.qza \ --i-models agora103_genus.qza \ --p-cutoff 0.0001 \ --p-threads 8 \ --o-community-models models.qza \ --verbose
This will give you something like this:
Merged with the database using ranks: genus Each community model contains 17-36 taxa (average 27+-6). 0%| | 0/16 [00:00<?, ?sample(s)/s]
On average, we a mean of 27 taxa per sample (and a range of 17-36). Building the community models will take a while, but this usually only needs to be done once for each sample. You can then run many growth simulations for a given sample using the existing model.
Now that we have our community models, let’s simulate growth. For this we will need information about the environment, because growth rates will depend strongly on available nutrients. In contrast to classical growth media in the lab, growth conditions in metabolic modeling have to be provided in terms of fluxes and not concentrations. Fluxomics or time course metabolomics data can be used to parameterize the growth environment, but these data are almost never available. Alternatively, you might estimate a minimal medium using the ‘omics data itself (more about that further down below). Usually, however, we estimate these parameters based what is known a priori about the nutrient environment. For example, there are existing estimates for a number of common diets (Virtual Metabolic Human). These dietary parameterizations are based on European diet data, with a MICOM-specific manual depletion of metabolites that are commonly absorbed in the small intestine. In this tutorial, we will use the “average Western diet”.
Additionally, you will have to provide a parameter value that represents the tradeoff between maximal community growth and individual taxon growth. For now we will set this to 0.3 and have a better look at this parameter in the next section. Let’s dive in and simulate growth!
qiime micom grow --i-models models.qza \ --i-medium western_diet_gut.qza \ --p-tradeoff 0.3 \ --p-threads 8 \ --o-results growth.qza \ --verbose
You will again see a progress bar and will have the results when everything is done (for me, this took about 5m with 4 threads).
We can now start to look at growth rates and fluxes in our models, but we will first come back to our tradeoff parameter…
One feature specific to MICOM is the tradeoff between community growth rate and individual taxon growth rates. Here, MICOM applies “pressure” to the model to allow growth for as many taxa as possible. However, this dual optimization problem results in tug-of-war between high community-level biomass production vs. allowing more taxa to grow. The balance between the two is dictated by a tradeoff value which dictates what percentage of the maximal community growth rate must be maintained, which ranges from 0 (no growth enforced) to 1 (enforce maximum community growth).
The truth usually lies in the middle of these extremes, and in our study we found that a tradeoff of 0.5 seemed to correspond best to in vivo replication rates. However, the best tradeoff may also depend on your samples, environment, and taxonomic rank used to build the models. So it is always a good idea to run a validation step. The following command will run the growth simulation with many different tradeoff values and track some key metrics.
qiime micom tradeoff --i-models models.qza \ --i-medium western_diet_gut.qza \ --p-threads 8 \ --o-results tradeoff.qza \ --verbose
After this analysis is finished we can take a look at our first visualization.
qiime micom plot-tradeoff --i-results tradeoff.qza \ --o-visualization tradeoff.qzv
This will give you the following:
Here the distribution of growth rates is shown by the 2D histogram on the left and the fraction of growing taxa with its mean line is shown on the right. You can see that lowering the tradeoff gives you more and more taxa that grow. The elbow is around 0.5, but we might want to pick a value as low as 0.3 here because that is where we observe the largest jump. Due to the constant dilution rate in the gut system, we expect most of the taxa at appreciable relative abundances should have a positive growth rate. We use this assumption to tune the tradeoff parameter. Another way to validate model growth rates, if you are working with shotgun metagenomic data, is to compare them against independently-estimated replication rates. Alternatively, non-growth could also be due to a breakdown between MICOM and reality (e.g. incorrect dietary parameterization or gaps in the genome-scale metabolic models that do not accurately represent the activity of the strains that are in your sample).
Okay, now that we understand where the tradeoff value comes from we can come back to our growth simulations.
MicomResults artifact we generated earlier can be visualized in many different ways. The first thing we might want to have a look at are the growth rates themselves.
qiime micom plot-growth --i-results growth.qza \ --o-visualization growth_rates.qzv
We can see that growth rates are pretty heterogeneous across samples, but it’s still hard to see what is going on.
Now let’s have a look at consumption of metabolites. Depending on the community composition, we should see varying preferences for consumption/production of certain metabolites and it would be interesting to see how these fluxes line up with our phenotype. In general, even for a given growth rate, these fluxes are not unique. For instance, if E. coli needs 10 mmol of glucose per hour to grow it may still import 100mmol per hour and just not use the remaining 90, or import 80 and not use the remaining 70, etc. To make those imports more unique MICOM will report the minimal consumption rates that still provide the observed growth rates (i.e. 10 mmol/h for our E. coli example).
qiime micom exchanges-per-sample --i-results growth.qza \ --o-visualization exchanges.qzv
We do see that there is some separation between healthy and cancer samples in the consumption of certain metabolites. For instance, there is a set of amino acids that get consumed by healthy gut microbiota but not so much in the cancer-associated gut microbiota. One of these amino acids is glutamine. Many cancer cells show glutamine addiction, as a consequence of the Warburg effect. So, we might hypothesize that glutamine is somewhat depleted in cancer samples and the associated microbiota have adapted to grow without it.
We could also look at production of metabolites by passing the
--p-direction export parameter, but due to the enforced minimum import there is usually very little net production of metabolites. So this is only useful if you want to look at overshoot production.
The last visualization gave us a global view of metabolite consumption but we may also be interested in looking at consumption preferences for individual taxa. We will need to lower the complexity of these data, as it is impossible to visualize imports for each metaboliote within each taxon in each sample. One useful way to simplify this visualization is to first reduce dimensions on the metabolite axis with UMAP and represent metabolite imports as a single point in a 2D space. Then we can plot each taxon in each specific samples onto that 2D space. Taxa that are close to each other occupy a similar growth niche, whereas taxa that are far away from each other are dissimilar in thier growth requirements.
qiime micom exchanges-per-taxon --i-results growth.qza \ --o-visualization niche.qzv
Most taxa have a specific growth niche, which varies slightly from sample to sample. You can tune the UMAP reduction by using the parameters
--p-min-dist. You can look at metabolite production with
--p-direction export where you can observe weaker clustering, suggesting that there is considerably more overlap across taxa in metabolite production.
A common question is whether there is a conection between a host phenotype of interest and the fluxes estimated by MICOM. However, not all exchange fluxes will be informative for this question. The most predictive type of flux for a phenotype of interest is usually the overall production flux, which is the total rate of production of a metabolite without considering its consumption by taxa in the community. For instance if all bacteria produce 100 mmol/gDW/h of ethanol but other taxa also consume 90 mmol/gDW/h of ethanol, the overall production flux would still be 100 mmol/gDW/h. Why is that production flux important? So far, MICOM only models microbial communities, without including cells from other Domains of life (e.g. Eukarya) that may be present in the environment. But host-associated bacteria compete with eukaryotic cells for metabolites, like intestinal cells in the gut or plant root cells in soil. These other cell types will consume a fraction of all metabolites in the environment and thus will take up metabolites from this total production pool.
Overall production fluxes are the default set of fluxes used by q2-micom. MICOM also allows for the analysis of minimal import fluxes shown in
qiime micom exchanges-per-sample, but these are often not very informative.
To predict a phenotype from fluxes we will use the
fit-phenotype command. Here you will need a Qiime 2 Metadata file that specifies your phenotype of interest. Our metadata is pretty simple and just assigns a healthy or cancer status to each sample:
id status Healthy1 healthy Healthy2 healthy Healthy3 healthy Healthy4 healthy Healthy5 healthy Healthy6 healthy Healthy7 healthy Healthy8 healthy Cancer1 colorectal cancer Cancer2 colorectal cancer Cancer3 colorectal cancer Cancer4 colorectal cancer Cancer5 colorectal cancer Cancer6 colorectal cancer Cancer7 colorectal cancer Cancer8 colorectal cancer
q2-micom supports binary or continuous phenotype data (using L1 penalized logistic regression for binary data, and LASSO regression for continuous data).
qiime micom fit-phenotype --i-results growth.qza \ --m-metadata-file metadata.tsv \ --m-metadata-column status \ --o-visualization fit.qzv
We could define a continuous phenotype with
--p-variable-type continuous or use for minimal import fluxes with
The visualization now shows which production fluxes are predictive of the phenotype.
Disclaimer: Depending on the version of scikit-learn in your installation this may look different.
Basically, the metabolites on both extremes of the barplot are the most predictive. In this case, negative coefficients denote fluxes that are higher in cancer samples and positive coefficients denote fluxes that are higher in healthy samples. For instance, here we see that butyrate, a helath-romoting, short-chain fatty acid is produced in slightly higher levels in helathy individuals. We also see a lot of sulfate-related compounds associated with cancer which is consistent with the proliferation-promoting properties of hydrogen sulfite (https://doi.org/10.1042/CBI20090368).
Obviously, we would still have to validate these very speculative hypotheses, but at least we get some functional insight starting from 16S data alone (while making several assumptions as detailed above…).
You can create q2-micom artifacts containing only a subset of samples using the
filter-results actions. For instance, if we wanted
to create model and results artifacts that contain only the cancer samples we could do the following:
qiime micom filter-models --i-models models.qza \ --m-metadata-file metadata.tsv \ --p-query "status == 'colorectal cancer'" \ --o-filtered-models cancer_built_models.qza
and similarly for the results:
qiime micom filter-results --i-results growth.qza \ --m-metadata-file metadata.tsv \ --p-query "status == 'colorectal cancer'" \ --o-filtered-results cancer_results.qza
Note that you will need a metadata table either (1) only includes the samples of interest or (2) can be
queried with a pandas query
to select only the samples of interest. You can also add the parameter
--p-exclude which will
return an artifact covering all samples except the ones selected by the metadata and query.
Those can now be passed on to the visualization steps to generate visualizations for only the cancer samples.
Many of the metabolite IDs returned by MICOM may be hard to identify. As long as you use the AGORA models you can look them up at https://www.vmh.life/#microbes/metabolites/ after stripping away the final
(e) which denotes the intracellular compartment. For instance to identify the metabolite
nmn(e) use https://www.vmh.life/#microbes/metabolites/nmn.
Maybe your sample is not well represented by the AGORA models, or you may have better metabolic reconstructions for the taxa in your samples. In that case you may want to build your own database for
q2-micom. This is pretty simple – you will need 2 things:
With these two things you can use
qiime micom db to build your database, summarized to any given taxonomic rank. For instance the AGORA v1.03 artifact was built by downloading the SBML models from https://www.vmh.life/#downloadview and using an adapted metadata file (agora.tsv). The genus database was then built using the following command.
qiime micom db --m-meta-file agora.tsv \ --p-rank genus \ --p-threads 8 \ --o-metabolic-models agora_genus_103.qza
Obviously, you may want to perform different statistical tests, visualizations, or analyses on the growth rates or fluxes obtained from
q2-micom. For this we provide a convenience function in the
q2_micom module, which returns the data as pandas DataFrames.
In : from q2_micom import read_results In : res = read_results("growth.qza") # there may be some warnings here, they are not dangerous :) In : res.exchange_fluxes.head() Out: taxon sample_id reaction flux abundance metabolite direction 0 Actinomyces Healthy4 EX_12dgr180(e) -1.436615e-11 0.001053 12dgr180(e) import 4 Actinomyces Healthy4 EX_12ppd_S(e) 9.727212e+00 0.001053 12ppd_S(e) export 10 Actinomyces Healthy4 EX_26dap_M(e) -4.315640e+01 0.001053 26dap_M(e) import 14 Actinomyces Healthy4 EX_2dmmq8(e) -3.792135e-12 0.001053 2dmmq8(e) import 18 Actinomyces Healthy4 EX_2obut(e) -6.305684e+01 0.001053 2obut(e) import In : res.growth_rates.head() Out: compartments abundance growth_rate reactions metabolites taxon tradeoff sample_id 0 Actinomyces 0.001053 1.167306e+00 2445 1643 Actinomyces 0.3 Healthy4 1 Adlercreutzia 0.011055 1.392794e-10 942 904 Adlercreutzia 0.3 Healthy4 2 Akkermansia 0.031848 5.686213e-11 2274 1386 Akkermansia 0.3 Healthy4 3 Alistipes 0.000464 4.795338e-01 2650 1613 Alistipes 0.3 Healthy4 4 Anaerococcus 0.002369 7.572441e-01 1330 1097 Anaerococcus 0.3 Healthy4
In some cases you may have no information on the environment and thus lack a starting point for parameterizing a growth medium. In this case
q2-micom, lets you predict a minimal medium based on all importable metabolites in the community models by specifying a minimum growth rate all taxa have to achieve simultaneously.
Be aware that this does not necessarily correspond to any realalistic environment that your microbial community is exposed to and may lead to very unrealistic predictions!
For instance we can generate a minimal medium for our community models:
qiime micom minimal-medium --i-models models.qza \ --p-min-growth 0.01 \ --p-threads 8 \ --o-medium minimal_medium.qza \ --verbose
We can use the Qiime 2 artifact API to inspect the medium:
In : from qiime2 import Artifact In : import pandas as pd In : medium = Artifact.load("minimal_medium.qza").view(pd.DataFrame) In : medium.shape Out: (33, 3) In : medium.sort_values(by="flux", ascending=False).head() Out: reaction flux metabolite 5 EX_amp_m 0.011940 amp_m 3 EX_MGlcn9_m 0.009677 MGlcn9_m 24 EX_nmn_m 0.007725 nmn_m 4 EX_MGlcn9_rl_m 0.007656 MGlcn9_rl_m 2 EX_MGlcn103_m 0.006159 MGlcn103_m
So, the most ‘efficient’ growth requires at least 33 different metabolites and seems to be consuming mucin glycans (all of the MGln* metabolites), as well as nucleotide and NAD precursors, like AMP and nicotinamide ribotide.