Here, we demonstrate the standard workflow of the SIAMCAT package using as
an example the dataset from Nielsen et al. Nat Biotechnol
2014.
This dataset contains samples from patients with inflammatory bowel disease and
from controls.
More importantly, these samples have been collected in two different countries,
Spain and Denmark. Together with technical differences between these samples,
this introduces a potent confounding factor into the data. Here we are going
to explore how SIAMCAT would identify the confounding factor and what the
results would be if you account for the confounder or not.
First, we load the packages needed to perform the analyses.
library("tidyverse")
library("SIAMCAT")
library("ggpubr")There are two different ways to access the data for our example dataset. On
the one hand, it is available through the curatedMetagenomicData R package.
However, using it here would create many more dependencies for the SIAMCAT
package.
Therefore, we here use data available through the EMBL cluster.
In the SIAMCAT paper, we performed the presented analyses on the datasets
available through curatedMetagenomicData. If you want to reproduce the
analysis from the SIAMCAT paper, you can execute the code chunks in
the curatedMetageomicData section, otherwise execute the code in the
mOTUs2 section.
First, we load the package:
library("curatedMetagenomicData")The data are part of the combined_metadata
meta.nielsen.full <- combined_metadata %>%
filter(dataset_name=='NielsenHB_2014')One thing we have to keep in mind are repeated samples per subject (see also the Machine learning pitfalls vignette).
print(length(unique(meta.nielsen.full$subjectID)))
print(nrow(meta.nielsen.full))Some subjects (but not all) had been sampled multiple times. Therefore, we want to remove repeated samplings for the same subject, since the samples would otherwise not be indepdenent from another.
The visit number is encoded in the sampleID. Therefore, we can use this
information to extract when the samples have been taken and use only the
first visit for each subject.
meta.nielsen <- meta.nielsen.full %>%
select(sampleID, subjectID, study_condition, disease_subtype,
disease, age, country, number_reads, median_read_length, BMI) %>%
mutate(visit=str_extract(sampleID, '_[0-9]+$')) %>%
mutate(visit=str_remove(visit, '_')) %>%
mutate(visit=as.numeric(visit)) %>%
mutate(visit=case_when(is.na(visit)~0, TRUE~visit)) %>%
group_by(subjectID) %>%
filter(visit==min(visit)) %>%
ungroup() %>%
mutate(Sample_ID=sampleID) %>%
mutate(Group=case_when(disease=='healthy'~'CTR',
TRUE~disease_subtype))Now, we can restrict our metadata to samples with UC and healthy control
samples:
meta.nielsen <- meta.nielsen %>%
filter(Group %in% c('UC', 'CTR'))As a last step, we can adjust the column names for the metadata so that they
agree with the data available from the EMBL cluster. Also, we add rownames to
the dataframe since SIAMCAT needs rownames to match samples across metadata
and features.
meta.nielsen <- meta.nielsen %>%
mutate(Country=country)
meta.nielsen <- as.data.frame(meta.nielsen)
rownames(meta.nielsen) <- meta.nielsen$sampleIDWe can load the taxonomic profiles generated with MetaPhlAn2 via the curatedMetagenomicsData R package.
x <- 'NielsenHB_2014.metaphlan_bugs_list.stool'
feat <- curatedMetagenomicData(x=x, dryrun=FALSE)
feat <- feat[[x]]@assayData$exprsThe MetaPhlAn2 profiles contain information on different taxonomic levels. Therefore, we want to restrict them to species-level profiles. In a second step, we convert them into relative abundances (summing up to 1) instead of using the percentages (summing up to 100) that MetaPhlAn2 outputs.
feat <- feat[grep(x=rownames(feat), pattern='s__'),]
feat <- feat[grep(x=rownames(feat),pattern='t__', invert = TRUE),]
feat <- t(t(feat)/100)The feature names are very long and may be a bit un-wieldy for plotting later on, so we shorten them to only the species name:
rownames(feat) <- str_extract(rownames(feat), 's__.*$')Both metadata and features are available through the EMBL cluster:
# base url for data download
data.loc <- 'https://zenodo.org/api/files/d81e429c-870f-44e0-a44a-2a4aa541b6c1/'
## metadata
meta.nielsen <- read_tsv(paste0(data.loc, 'meta_Nielsen.tsv'))
## Rows: 396 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (5): Sample_ID, Individual_ID, Country, Gender, Group
## dbl (4): Sampling_day, Age, BMI, Library_Size
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# also here, we have to remove repeated samplings and CD samples
meta.nielsen <- meta.nielsen %>%
filter(Group %in% c('CTR', 'UC')) %>%
group_by(Individual_ID) %>%
filter(Sampling_day==min(Sampling_day)) %>%
ungroup() %>%
as.data.frame()
rownames(meta.nielsen) <- meta.nielsen$Sample_ID
## features
feat <- read.table(paste0(data.loc, 'metaHIT_motus.tsv'),
stringsAsFactors = FALSE, sep='\t',
check.names = FALSE, quote = '', comment.char = '')
feat <- feat[,colSums(feat) > 0]
feat <- prop.table(as.matrix(feat), 2)Now, we have everything ready to create a SIAMCAT object which stores
the feature matrix, the meta-variables, and the label. Here, the label is
created using the information in the metadata.
To demonstrate the normal SIAMCAT workflow, we remove the confounding factor
by only looking at samples from Spain. Below, we have a look what would have
happened if we had not removed them.
# remove Danish samples
meta.nielsen.esp <- meta.nielsen[meta.nielsen$Country == 'ESP',]
sc.obj <- siamcat(feat=feat, meta=meta.nielsen.esp, label='Group', case='UC')
## + starting create.label
## Label used as case:
## UC
## Label used as control:
## CTR
## + finished create.label.from.metadata in 0.004 s
## + starting validate.data
## +++ checking overlap between labels and features
## + Keeping labels of 128 sample(s).
## +++ checking sample number per class
## +++ checking overlap between samples and metadata
## + finished validate.data in 0.647 sNow, we can filter feature with low overall abundance and prevalence.
sc.obj <- filter.features(sc.obj, cutoff=1e-04,
filter.method = 'abundance')
## Features successfully filtered
sc.obj <- filter.features(sc.obj, cutoff=0.05,
filter.method='prevalence',
feature.type = 'filtered')
## Features successfully filteredThe check.assocation function calculates the significance of enrichment and
metrics of association (such as generalized fold change and single-feature
AUROC).
sc.obj <- check.associations(sc.obj, log.n0 = 1e-06, alpha=0.1)
association.plot(sc.obj, fn.plot = './association_plot_nielsen.pdf',
panels = c('fc'))