# SIAMCAT: Statistical Inference of Associations between Microbial Communities And host phenoTypes

Konrad Zych, Jakob Wirbel, and Georg Zeller1*

1EMBL Heidelberg

*georg.zeller@embl.de

#### Date last modified: 2018-09-24

# Contents

# 1 About This Vignette

This vignette aims to be a short tutorial for the main functionalities of `SIAMCAT`. Examples of additional workflows or more detailed tutorials can be found in other vignettes (see the BioConductor page).

`SIAMCAT` is part of the suite of computational microbiome analysis tools hosted at EMBL by the groups of Peer Bork and Georg Zeller. Find out more at EMBL-microbiome tools.

# 2 Introduction

Associations between microbiome and host phenotypes are ideally described by quantitative models able to predict host status from microbiome composition. `SIAMCAT` can do so for data from hundreds of thousands of microbial taxa, gene families, or metabolic pathways over hundreds of samples. `SIAMCAT` produces graphical output for convenient assessment of the quality of the input data and statistical associations, for model diagnostics and inference revealing the most predictive microbial biomarkers.

# 3 Quick Start

For this vignette, we use an example dataset included in the `SIAMCAT` package. As example dataset we use the data from the publication of Zeller et al, which demonstrated the potential of microbial species in fecal samples to distinguish patients with colorectal cancer (CRC) from healthy controls.

    library(SIAMCAT)

    data("feat_crc_zeller", package="SIAMCAT")
    data("meta_crc_zeller", package="SIAMCAT")

First, `SIAMCAT` needs a feature matrix (can be either a `matrix`, a `data.frame`, or a `phyloseq-otu_table`), which contains values of different features (in rows) for different samples (in columns). For example, the feature matrix included here contains relative abundances for bacterial species calculated with the mOTU profiler for 141 samples:

    feat.crc.zeller[1:3, 1:3]

    ##                                  CCIS27304052ST-3-0 CCIS15794887ST-4-0
    ## UNMAPPED                                   0.589839          0.7142157
    ## Methanoculleus marisnigri [h:1]            0.000000          0.0000000
    ## Methanococcoides burtonii [h:10]           0.000000          0.0000000
    ##                                  CCIS74726977ST-3-0
    ## UNMAPPED                                  0.7818674
    ## Methanoculleus marisnigri [h:1]           0.0000000
    ## Methanococcoides burtonii [h:10]          0.0000000

    dim(feat.crc.zeller)

    ## [1] 1754  141

> Please note that `SIAMCAT` is supposed to work with relative abundances. Other types of data (e.g. counts) will also work, but not all functions of the package will result in meaningful outputs.

Secondly, we also have metadata about the samples in another `data.frame`:

    head(meta.crc.zeller)

    ##                    Age BMI Gender AJCC_stage     FOBT Group
    ## CCIS27304052ST-3-0  52  20      F         -1 Negative   CTR
    ## CCIS15794887ST-4-0  37  18      F         -1 Negative   CTR
    ## CCIS74726977ST-3-0  66  24      M         -1 Negative   CTR
    ## CCIS16561622ST-4-0  54  26      M         -1 Negative   CTR
    ## CCIS79210440ST-3-0  65  30      M         -1 Positive   CTR
    ## CCIS82507866ST-3-0  57  24      M         -1 Negative   CTR

In order to tell `SIAMCAT`, which samples are cancer cases and which are healthy controls, we can construct a label object from the `Group` column in the metadata.

    label.crc.zeller <- create.label(meta=meta.crc.zeller,
        label='Group', case='CRC')

    ## Label used as case:
    ##    CRC
    ## Label used as control:
    ##    CTR

    ## + finished create.label.from.metadata in    0 s

Now we have all the ingredients to create a `SIAMCAT` object. Please have a look at the vignette about input formats for more information about supported formats and other ways to create a `SIAMCAT` object.

    siamcat <- siamcat(feat=feat.crc.zeller,
        label=label.crc.zeller,
        meta=meta.crc.zeller)

    ## + starting validate.data

    ## +++ checking overlap between labels and features

    ## + Keeping labels of 141 sample(s).

    ## +++ checking sample number per class

    ## +++ checking overlap between samples and metadata

    ## + finished validate.data in 0.17 s

A few information about the `siamcat` object can be accessed with the `show` function from `phyloseq` (`SIAMCAT` builds on the `phyloseq` data structure):

    show(siamcat)

    ## siamcat-class object
    ## label()                Label object:         88 CTR and 53 CRC samples
    ##
    ## contains phyloseq-class experiment-level object @phyloseq:
    ## phyloseq@otu_table()   OTU Table:            [ 1754 taxa and 141 samples ]
    ## phyloseq@sam_data()    Sample Data:          [ 141 samples by 6 sample variables ]

Since we have quite a lot of microbial species in the dataset at the moment, we can perform unsupervised feature selection using the function `filter.features`.

    siamcat <- filter.features(siamcat,
        filter.method = 'abundance',
        cutoff = 0.001)

    ## Features successfully filtered

# 4 Association Testing

Associations between microbial species and the label can be tested with the `check.associations` function. The function computes for each species the significance using a non-parametric Wilcoxon test and different effect sizes for the association (e.g. AUC or fold change).

    siamcat <- check.associations(
        siamcat,
        sort.by = 'fc',
        alpha = 0.05,
        mult.corr = "fdr",
        detect.lim = 10 ^-6,
        plot.type = "quantile.box",
        panels = c("fc", "prevalence", "auroc"))

The function produces a pdf file as output, since the plot is optimized for a landscape DIN-A4 layout, but can also used to plot on an active graphic device, e.g. in RStudio. The resulting plot then looks like that:

# 5 Model Building

One strength of `SIAMCAT` is the versatile but easy-to-use interface for the construction of machine learning models on the basis of microbial species. `SIAMCAT` contains functions for data normalization, splitting the data into cross-validation folds, training the model, and making predictions based on cross-validation instances and the trained models.

## 5.1 Data Normalization

Data normalization is performed with the `normalize.features` function. Here, we use the `log.unit` method, but several other methods and customization options are available (please check the documentation).

    siamcat <- normalize.features(
        siamcat,
        norm.method = "log.unit",
        norm.param = list(
            log.n0 = 1e-06,
            n.p = 2,
            norm.margin = 1
        )
    )

    ## Features normalized successfully.

## 5.2 Prepare Cross-Validation

Preparation of the cross-validation fold is a crucial step in machine learning. `SIAMCAT` greatly simplifies the set-up of cross-validation schemes, including stratification of samples or keeping samples inseperable based on metadata. For this small example, we choose a twice-repeated 5-fold cross-validation scheme. The data-split will be saved in the `data_split` slot of the `siamcat` object.

    siamcat <-  create.data.split(
        siamcat,
        num.folds = 5,
        num.resample = 2
    )

    ## Features splitted for cross-validation successfully.

## 5.3 Model Training

The actual model training is performed using the function `train.model`. Again, multiple options for customization are available, ranging from the machine learning method to the measure for model selection or customizable parameter set for hyperparameter tuning.

    siamcat <- train.model(
        siamcat,
        method = "lasso"
    )

The models are saved in the `model_list` slot of the `siamcat` object. The model building is performed using the `mlr` R package. All models can easily be accessed.

    # get information about the model type
    model_type(siamcat)

    ## [1] "lasso"

    # access the models
    models <- models(siamcat)
    models[[1]]

    ## Model for learner.id=classif.cvglmnet; learner.class=classif.cvglmnet
    ## Trained on: task.id = data; obs = 112; features = 207
    ## Hyperparameters: nlambda=100,alpha=1

## 5.4 Make Predictions

Using the data-split and the models trained in previous step, we can use the function `make.predictions` in order to apply the models on the test instances in the data-split. The predictions will be saved in the `pred_matrix` slot of the `siamcat` object.

    siamcat <- make.predictions(siamcat)
    pred_matrix <- pred_matrix(siamcat)

    head(pred_matrix)

    ##                      CV_rep1    CV_rep2
    ## CCIS27304052ST-3-0 0.1355932 0.04910265
    ## CCIS15794887ST-4-0 0.1052626 0.06721616
    ## CCIS74726977ST-3-0 0.6010033 0.45839135
    ## CCIS16561622ST-4-0 0.2552595 0.14979968
    ## CCIS79210440ST-3-0 0.2767179 0.57783899
    ## CCIS82507866ST-3-0 0.1456438 0.05918578

# 6 Model Evaluation and Interpretation

In the final part, we want to find out how well the model performed and which microbial species had been selected in the model. In order to do so, we first calculate how well the predictions fit the real data using the function `evaluate.predictions`. This function calculates the Area Under the Receiver Operating Characteristic (ROC) Curve (AU-ROC) and the Precision Recall (PR) Curve for each resampled cross-validation run.

    siamcat <-  evaluate.predictions(siamcat)

    ## Evaluated predictions successfully.

## 6.1 Evaluation Plot

To plot the results of the evaluation, we can use the function `model.evaluation.plot`, which produces a pdf-file showing the ROC and PR Curves for the different resamples runs as well as the mean ROC and PR Curve.

    model.evaluation.plot(siamcat)

## 6.2 Interpretation Plot

The final plot produced by `SIAMCAT` is the model interpretation plot, created by the `model.interpretation.plot` function. The plot shows for the top selected features the

- model weights (and how robust they are) as a barplot,

- a heatmap with the z-scores or fold changes for the top selected features, and

- a boxplot showing the proportions of weight per model which is captured by the top selected features.

Additionally, the distribution of metadata is shown in a heatmap below.

The function again produces a pdf-file optimized for a landscape DIN-A4 plotting region.

    model.interpretation.plot(
        siamcat,
        fn.plot = 'interpretation.pdf',
        consens.thres = 0.5,
        norm.models = TRUE,
        limits = c(-3, 3),
        heatmap.type = 'zscore',
    )

The resulting plot looks like this:

# 7 Session Info

    sessionInfo()

    ## R version 3.5.2 (2018-12-20)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] SIAMCAT_1.2.1     phyloseq_1.26.1   mlr_2.13          ParamHelpers_1.11
    ## [5] BiocStyle_2.10.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] Biobase_2.42.0      jsonlite_1.6        splines_3.5.2
    ##  [4] foreach_1.4.4       LiblineaR_2.10-8    PRROC_1.3.1
    ##  [7] assertthat_0.2.0    BiocManager_1.30.4  stats4_3.5.2
    ## [10] yaml_2.2.0          corrplot_0.84       pillar_1.3.1
    ## [13] backports_1.1.3     lattice_0.20-38     glue_1.3.0
    ## [16] pROC_1.13.0         digest_0.6.18       RColorBrewer_1.1-2
    ## [19] XVector_0.22.0      checkmate_1.8.5     colorspace_1.3-2
    ## [22] htmltools_0.3.6     Matrix_1.2-15       plyr_1.8.4
    ## [25] XML_3.98-1.16       pkgconfig_2.0.2     bookdown_0.9
    ## [28] zlibbioc_1.28.0     purrr_0.2.5         scales_1.0.0
    ## [31] parallelMap_1.3     tibble_2.0.0        beanplot_1.2
    ## [34] mgcv_1.8-26         IRanges_2.16.0      ggplot2_3.1.0
    ## [37] infotheo_1.2.0      BiocGenerics_0.28.0 lazyeval_0.2.1
    ## [40] survival_2.43-3     magrittr_1.5        crayon_1.3.4
    ## [43] evaluate_0.12       nlme_3.1-137        MASS_7.3-51.1
    ## [46] vegan_2.5-3         tools_3.5.2         data.table_1.11.8
    ## [49] matrixStats_0.54.0  gridBase_0.4-7      BBmisc_1.11
    ## [52] stringr_1.3.1       Rhdf5lib_1.4.2      S4Vectors_0.20.1
    ## [55] munsell_0.5.0       glmnet_2.0-16       cluster_2.0.7-1
    ## [58] bindrcpp_0.2.2      Biostrings_2.50.2   ade4_1.7-13
    ## [61] compiler_3.5.2      rlang_0.3.0.1       rhdf5_2.26.2
    ## [64] grid_3.5.2          iterators_1.0.10    biomformat_1.10.1
    ## [67] igraph_1.2.2        rmarkdown_1.11      gtable_0.2.0
    ## [70] codetools_0.2-16    multtest_2.38.0     reshape2_1.4.3
    ## [73] R6_2.3.0            gridExtra_2.3       knitr_1.21
    ## [76] dplyr_0.7.8         bindr_0.1.1         fastmatch_1.1-0
    ## [79] permute_0.9-4       ape_5.2             stringi_1.2.4
    ## [82] parallel_3.5.2      Rcpp_1.0.0          tidyselect_0.2.5
    ## [85] xfun_0.4