Minimal gut bioreactor examples

Minimal gut microbiome

Dense samples of the minimal gut microbiome. In the initial hours, MDb-MM was grown under batch condition and 24 h onwards, continuous feeding of media with pulse feeding cycles. This information is stored in the colData.

library(miaTime)
#> Loading required package: mia
#> Loading required package: MultiAssayExperiment
#> Loading required package: SummarizedExperiment
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
#>     table, tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians
#> Loading required package: SingleCellExperiment
#> Loading required package: TreeSummarizedExperiment
#> Loading required package: Biostrings
#> Loading required package: XVector
#> 
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#> 
#>     strsplit
#> This is mia version 1.15.3
#> - Online documentation and vignettes: https://microbiome.github.io/mia/
#> - Online book 'Orchestrating Microbiome Analysis (OMA)': https://microbiome.github.io/OMA/docs/devel/

data(minimalgut)
tse <- minimalgut

# Quick check of number of samples 
table(tse[["StudyIdentifier"]], tse[["condition_1"]])
#>               
#>                batch_carbs DoS pulse Overnight
#>   Bioreactor A           4        38        19
#>   Bioreactor B           4        38        19
#>   Bioreactor C           4        38        19

Visualize samples available for each of the bioreactors. This allows to identify if there are any missing samples for specific times.

library(ggplot2)

colData(tse) |>
    ggplot() +
    geom_tile(
        aes(x = as.factor(Time.hr), y = StudyIdentifier, fill = condition_1))

Community dynamics

The minimalgut dataset, mucus-diet based minimal microbiome (MDbMM-16), consists of 16 species assembled in three bioreactors. We can investigate the succession of mdbMM16 from the start of experiment here hour zero until the end of the experiment.

# Transform data to relativeS
tse <- transformAssay(tse, method = "relabundance")
# Divergence from baseline i.e from hour zero
tse <- addBaselineDivergence(
    tse,
    assay.type = "relabundance",
    method = "bray",
    group = "StudyIdentifier",
    time.col = "Time.hr",
    )

Let’s then visualize the divergence.

library(scater)
#> Loading required package: scuttle

# Create a time series plot for divergence
p <- plotColData(
    tse, x = "Time.hr", y = "divergence", colour_by = "StudyIdentifier") +
    # Add line between points
    geom_line(aes(group = .data[["colour_by"]], colour =  .data[["colour_by"]]))
p

Visualizing selected taxa

Now visualize abundance of Blautia hydrogenotrophica using the miaViz::plotSeries() function.

library(miaViz)
#> Loading required package: ggraph
#> 
#> Attaching package: 'miaViz'
#> The following object is masked from 'package:mia':
#> 
#>     plotNMDS

# Plot certain feature by time
p <- plotSeries(
    tse,
    x = "Time.hr", y = "Blautia_hydrogenotrophica", colour_by = "Species",
    assay.type = "relabundance")
p

Visualize the rate (slope) of divergence

Sample dissimilarity between consecutive time steps(step size n >= 1) within a group(subject, age, reaction chamber, etc.) can be calculated by addStepwiseDivergence.

# Divergence between consecutive time points
tse <- addStepwiseDivergence(
    tse,
    assay.type = "relabundance",
    method = "bray",
    group = "StudyIdentifier",
    time.interval = 1,
    time.col = "Time.hr",
    name = "divergence_from_previous_step",
    name.time = "time_from_previous_step"
    )

The results are again stored in colData. We calculate the speed of divergence change by dividing each divergence change by the corresponding change in time. Then we use similar plotting methods as previously.

# Calculate slope for the change
tse[["divergence_change"]] <- tse[["divergence_from_previous_step"]] /
    tse[["time_from_previous_step"]]

# Create a time series plot for divergence
p <- plotColData(
    tse,
    x = "Time.hr",
    y = "divergence_change",
    colour_by = "StudyIdentifier"
    ) +
    # Add line between points
    geom_line(aes(group = .data[["colour_by"]], colour =  .data[["colour_by"]]))
p
#> Warning: Removed 3 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Warning: Removed 3 rows containing missing values or values outside the scale range
#> (`geom_line()`).

Moving average of the slope

This shows how to calculate and plot moving average for the variable of interest (here: slope).

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:mia':
#> 
#>     full_join, inner_join, left_join, right_join
#> The following objects are masked from 'package:Biostrings':
#> 
#>     collapse, intersect, setdiff, setequal, union
#> The following object is masked from 'package:XVector':
#> 
#>     slice
#> The following object is masked from 'package:Biobase':
#> 
#>     combine
#> The following objects are masked from 'package:GenomicRanges':
#> 
#>     intersect, setdiff, union
#> The following object is masked from 'package:GenomeInfoDb':
#> 
#>     intersect
#> The following objects are masked from 'package:IRanges':
#> 
#>     collapse, desc, intersect, setdiff, slice, union
#> The following objects are masked from 'package:S4Vectors':
#> 
#>     first, intersect, rename, setdiff, setequal, union
#> The following objects are masked from 'package:BiocGenerics':
#> 
#>     combine, intersect, setdiff, union
#> The following object is masked from 'package:matrixStats':
#> 
#>     count
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

# Calculate moving average with time window of 3 time points
tse[["sliding_divergence"]] <- colData(tse) |>
    as.data.frame() |>
    # Group based on reactor
    group_by(StudyIdentifier) |>
    # Calculate moving average
    mutate(sliding_avg = (
        # We get the previous 2 samples
        lag(divergence_change, 2) +
        lag(divergence_change, 1) +
        # And the current sample
        divergence_change
        # And take average
        ) / 3
    ) |>
    # Get only the values as vector
    ungroup() |>
    pull(sliding_avg)

After calculating the moving average of divergences, we can visualize the result in a similar way to our previous approach.

# Create a time series plot for divergence
p <- plotColData(
    tse,
    x = "Time.hr",
    y = "sliding_divergence",
    colour_by = "StudyIdentifier"
    ) +
    # Add line between points
    geom_line(aes(group = .data[["colour_by"]], colour =  .data[["colour_by"]]))
p
#> Warning: Removed 9 rows containing missing values or values outside the scale range
#> (`geom_point()`).
#> Warning: Removed 3 rows containing missing values or values outside the scale range
#> (`geom_line()`).

Session info