Skip to content

Latest commit

 

History

History
124 lines (92 loc) · 7.56 KB

exercise_2.md

File metadata and controls

124 lines (92 loc) · 7.56 KB

🕸️ 2. Simulate community interactions between bacteria and yeast

We will now use flux balance analysis based methods to enumerate feasible interactions between metabolic networks representing community members. Refer to the SMETANA documentation for more details regarding usage.

$ smetana -h
Click to see helpfile

usage: smetana [-h] [-c COMMUNITIES.TSV] [-o OUTPUT] [--flavor FLAVOR]
               [-m MEDIA] [--mediadb MEDIADB]
               [-g | -d | -a ABIOTIC | -b BIOTIC] [-p P] [-n N] [-v] [-z]
               [--solver SOLVER] [--molweight] [--exclude EXCLUDE]
               [--no-coupling]
               MODELS [MODELS ...]

Calculate SMETANA scores for one or multiple microbial communities.

positional arguments:
  MODELS                
                        Multiple single-species models (one or more files).
                        
                        You can use wild-cards, for example: models/*.xml, and optionally protect with quotes to avoid automatic bash
                        expansion (this will be faster for long lists): "models/*.xml". 

optional arguments:
  -h, --help            show this help message and exit
  -c COMMUNITIES.TSV, --communities COMMUNITIES.TSV
                        
                        Run SMETANA for multiple (sub)communities.
                        The communities must be specified in a two-column tab-separated file with community and organism identifiers.
                        The organism identifiers should match the file names in the SBML files (without extension).
                        
                        Example:
                            community1	organism1
                            community1	organism2
                            community2	organism1
                            community2	organism3
                        
  -o OUTPUT, --output OUTPUT
                        Prefix for output file(s).
  --flavor FLAVOR       Expected SBML flavor of the input files (cobra or fbc2).
  -m MEDIA, --media MEDIA
                        Run SMETANA for given media (comma-separated).
  --mediadb MEDIADB     Media database file
  -g, --global          Run global analysis with MIP/MRO (faster).
  -d, --detailed        Run detailed SMETANA analysis (slower).
  -a ABIOTIC, --abiotic ABIOTIC
                        Test abiotic perturbations with given list of compounds.
  -b BIOTIC, --biotic BIOTIC
                        Test biotic perturbations with given list of species.
  -p P                  Number of components to perturb simultaneously (default: 1).
  -n N                  
                        Number of random perturbation experiments per community (default: 1).
                        Selecting n = 0 will test all single species/compound perturbations exactly once.
  -v, --verbose         Switch to verbose mode
  -z, --zeros           Include entries with zero score.
  --solver SOLVER       Change default solver (current options: 'gurobi', 'cplex').
  --molweight           Use molecular weight minimization (recomended).
  --exclude EXCLUDE     List of compounds to exclude from calculations (e.g.: inorganic compounds).
  --no-coupling         Don't compute species coupling scores.

🤝 SMETANA detailed algorithm enumerates all feasible interactions between metabolic networks

Refer to the methods sections of the SMETANA paper for details regarding the implementation of MILP probelms that are solved in order to compute the following:

  1. The species coupling score (SCS) measures the dependence of growth of species A on species B (SCSA,B)
    • calculated by enumerating all possible community member subsets where species A can grow, SCSA,B is the fraction of subsets where both species A and B can grow.
  2. The metabolite uptake score (MUS) measures the dependence of growth of species A on metabolite m (MUSA,m)
    • calculated by enumerating all possible metabolite requirement subsets where species A can grow, MUSA,m is the fraction of subsets where both species A grows and metabolite m is taken up.
  3. The metabolite production score (MPS) is a binary score indicating whether a given species B can produce metabolite m (MPS = 1) or not (MPS = 0) in the community of N members (MPSB,m)
  4. The SMETANA score ranges from 0 to 1
    • measures how strongly a receiver species relies on a donor species for a particular metabolite
    • SMETANAA,B,m = SCSA,B * MUSA,m * MPSB,m

🥼 Run 3 computational experiments

Run each simulation on your own machine as shown in the following code chunks.

🐄 CDM35 with lactose as main carbon source

$ smetana -d -v --flavor cobra --mediadb $ROOT/media.tsv -m CDM35_lcts -o $ROOT/CDM35_lcts $ROOT/models/*.xml && paste $ROOT/CDM35_lcts_detailed.tsv 

🍬 CDM35 with glucose as main carbon source

$ smetana -d -v --flavor cobra --mediadb $ROOT/media.tsv -m CDM35_glc -o $ROOT/CDM35_glc $ROOT/models/*.xml && paste $ROOT/CDM35_glc_detailed.tsv 

🌌 CDM35 with galactose as main carbon source

$ smetana -d -v --flavor cobra --mediadb $ROOT/media.tsv -m CDM35_gal -o $ROOT/CDM35_gal $ROOT/models/*.xml && paste $ROOT/CDM35_gal_detailed.tsv 

📺 Visualizing interactions

Alluvial diagrams, AKA sankey or flow diagrams, are a quick and easy way to visualize interaction predictions generated by SMETANA. We provide the following Rscript to be used for plotting your generated results.

Click to see Rscript

# install ggalluvial if missing
install.packages("ggalluvial")
# load libraries
library(tidyverse)
library(ggalluvial)
# load data, replace filename with your results to generate new plot
read.delim("CDM35_lcts_detailed.tsv") -> yeastlactis_smetana # SMETANA results
read.delim("bigg_classes.tsv") -> bigg_mets # BiGG metabolite ids information
# plot, play around with the parameters in this function to modify plot
ggplot(yeastlactis_smetana %>%
mutate(compound = gsub("M_","",compound),compound = gsub("_e","",compound)) %>%
left_join(.,bigg_mets) %>% filter(smetana >= 0.2,name!="DL-Glutamate"),
aes(axis1 = donor, axis2 = name, axis3 = receiver,
y = smetana)) +
scale_x_discrete(limits = c("Donor", "Metabolite", "Reciever")) +
xlab("Interaction") +
geom_alluvium(aes(fill = name)) +
geom_stratum(width=0.3) +
theme_minimal() + geom_text(stat = "stratum", aes(label = after_stat(stratum)),min.y=0.2)+theme_bw() +
theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),axis.line.y = element_blank(),axis.ticks.y = element_blank(),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.line.x = element_blank(),axis.ticks.x = element_blank(),legend.position = "none")

To start plotting your own data, first open Rstudio by clicking on the bottom left icon and searching for Rstudio. Then click to navigate to your $ROOT folder in the bottom right file explorer. Click on the little gear button labeled More and set your current folder as the working directory. Then click on the plot_interactions.R script, which should bring up the code shown above. Now you may run the script with the default parameters to generate an alluvial diagram.

The script is hardcoded to read the CDM35 lactose SMETANA results file CDM35_lcts_detailed.tsv and plot to a file called CDM35_lcts.pdf. After generating prediction for each media condition, modify these parameters to read in and generate plots for each of the above experimental conditions. You may also have a look at some precomputed results for comparison, or if you run into any trouble generating results of your own.

💎 Discussion questions

  • How does switching carbon source affect the metabolic interactions between yeast and bacteria?
  • To what extent do these predicted interactions reflect those experimentally verified in the case study publication?
  • Come up with some variation of CDM35 media by removing glucose/lactose/galactose and adding a new carbon source, e.g. trehalose, maltose, etc. Then enumerate interactions with SMETANA in that media, plot, and compare the interactions to the simulations above. How do the interactions look? Was there something unexpected?

Move on to exercise 3