Using FLORAL for Microbiome Analysis
library(FLORAL)
library(dplyr)
library(patchwork)
Data
We will be using data from the curatedMetagenomicData
package. For easier installation, we saved a flat copy of the data, but
the steps below show how that was created.
if (! "BiocManager" %in% installed.packages()) install.packages("BiocManager")
if (! "curatedMetagenomicData" %in% installed.packages()) BiocManager::install("curatedMetagenomicData")
if (! "patchwork" %in% installed.packages()) install.packages("patchwork")
library(curatedMetagenomicData)
# Take a look at the summary of the studies available:
curatedMetagenomicData::sampleMetadata |> group_by(.data$study_name, .data$study_condition) |> count() |> arrange(.data$study_name)
#As an example, let us look at the `YachidaS_2019` study between healthy controls and colorectal cancer (CRC) patients.
curatedMetagenomicData::curatedMetagenomicData("YachidaS_2019")
# note -- if you are behind a firewall, see the solutions to 500 errors here:
# https://support.bioconductor.org/p/132864/
rawdata <- curatedMetagenomicData::curatedMetagenomicData("2021-03-31.YachidaS_2019.relative_abundance", dryrun = FALSE, counts = TRUE) |> mergeData()
x <- SummarizedExperiment::assays(rawdata)$relative_abundance %>% t()
y <- rawdata@colData$disease
save(list = c("x", "y"), file = file.path("inst", "extdata", "YachidaS_2019.Rdata"))
load(system.file("extdata", "YachidaS_2019.Rdata", package="FLORAL"))
Running FLORAL
We extracted the data from the TreeSummarizedExperiment
object to two objects: the taxa matrix x and the “outcomes”
vector y of whether a patient is healthy or has colorectal
cancer (CRC). Note that for binary outcomes, the input vector
y needs to be formatted with entries equal to either
0 or 1. In addition, we need to specify
family = "binomial" in FLORAL to fit the
logistic regression model. To print the progress bar as the algorithm
runs, please use progress = TRUE.
x <- x[y %in% c("CRC","healthy"),]
x <- x[,colSums(x >= 100) >= nrow(x)*0.2] # filter low abundance taxa
colnames(x) <- sapply(colnames(x), function(x) strsplit(x,split="[|]")[[1]][length(strsplit(x,split="[|]")[[1]])])
y <- as.numeric(as.factor(y[y %in% c("CRC","healthy")]))-1
fit <- FLORAL(x = x, y = y, family="binomial", ncv=10, progress=TRUE)
#> Using elastic net with a=1.Algorithm running for full dataset:
#> Algorithm running for cv dataset 1 out of 10:
#> Algorithm running for cv dataset 2 out of 10:
#> Algorithm running for cv dataset 3 out of 10:
#> Algorithm running for cv dataset 4 out of 10:
#> Algorithm running for cv dataset 5 out of 10:
#> Algorithm running for cv dataset 6 out of 10:
#> Algorithm running for cv dataset 7 out of 10:
#> Algorithm running for cv dataset 8 out of 10:
#> Algorithm running for cv dataset 9 out of 10:
#> Algorithm running for cv dataset 10 out of 10:
Interpreting the Model
FLORAL, like other methods that have an optimization step, has two
“best” solutions for \(\lambda\)
available: one minimizing the mean squared error (\(\lambda_\min\)), and one maximizing the
value of \(\lambda\) withing 1 standard
error of the minimum mean squared error (\(\lambda_{\text{1se}}\)). These are referred
to as the min and 1se solutions,
respectively.
We can see the mean squared error (MSE) and the coefficients vs
log(\(\lambda\)) as follows:
In both plots, the vertical dashed line and dotted line represent
\(\lambda_\min\) and \(\lambda_{\text{1se}}\), respectively. In
the MSE plot, the bands represent plus minus one standard error of the
MSE. In the coefficient plot, the colored lines represent individual
taxa, where taxa with non-zero values at \(\lambda_\min\) and \(\lambda_{\text{1se}}\) are selected as
predictive of the outcome.
To view specific names of the selected taxa, please see
fit$selected.feature$min or
fit$selected.feature$1se vectors. To view all coefficient
estimates, please see fit$best.beta$min or
fit$best.beta$1se. Without looking into ratios, one can
crudely interpret positive or negative association between a taxon and
the outcome by the positive or negative sign of the coefficient
estimates. However, we recommend referring to the two-step procedure
discussed below for a more rigorous interpretation based on ratios,
which is derived from the log-ratio model assumption.
head(fit$selected.feature$min)
#> [1] "s__Actinomyces_sp_HMSC035G02" "s__Actinomyces_sp_ICM47"
#> [3] "s__Anaerotignum_lactatifermentans" "s__Anaerotruncus_colihominis"
#> [5] "s__Asaccharobacter_celatus" "s__Bacteroides_caccae"
head(sort(fit$best.beta$min))
#> s__Parvimonas_micra s__Actinomyces_sp_HMSC035G02
#> -0.08252935 -0.05314047
#> s__Collinsella_aerofaciens s__Holdemania_filiformis
#> -0.04716041 -0.03963238
#> s__Anaerotruncus_colihominis s__Bacteroides_plebeius
#> -0.03450502 -0.02348108
The Two-step Procedure
In the previous section, we checked the lasso estimates without
identifying specific ratios that are predictive of the outcome (CRC in
this case). By default, FLORAL performs a two-step
selection procedure to use glmnet and step
regression to further identify taxa pairs which form predictive
log-ratios. To view those pairs, use fit$step2.ratios$min
or fit$step2.ratios$1se for names of ratios and
fit$step2.ratios$min.idx or
fit$step2.ratios$1se.idx for the pairs of indices in the
original input count matrix x. Note that one taxon can
occur in multiple ratios.
head(fit$step2.ratios$`1se`)
#> [1] "s__Bacteroides_plebeius/s__Bacteroides_dorei"
#> [2] "s__Eubacterium_eligens/s__Bacteroides_cellulosilyticus"
#> [3] "s__Collinsella_aerofaciens/s__Bifidobacterium_pseudocatenulatum"
#> [4] "s__Streptococcus_salivarius/s__Actinomyces_sp_HMSC035G02"
#> [5] "s__Erysipelatoclostridium_ramosum/s__Parvimonas_micra"
#> [6] "s__Holdemania_filiformis/s__Eubacterium_sp_CAG_38"
fit$step2.ratios$`1se.idx`
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
#> [1,] NA 1 4 12 22 26 27 27 44 73 81 114
#> [2,] NA 2 113 102 128 131 90 151 64 80 128 131
To further interpret the positive or negative associations between
the outcome, please refer to the output step regression
tables, where the effect sizes of the ratios can be found.
While the corresponding p-values are also available, we recommend
only using the p-values as a criterion to rank the strength of the
association. We do not recommend directly reporting the p-values for
inference, because these p-values were obtained after running the first
step lasso model without rigorous post-selective inference. However, it
is still valid to claim these selected log-ratios are predictive of the
outcome, as demonstrated by the improved 10-fold cross-validated
prediction errors.
fit$step2.tables$`1se`
#> Estimate
#> (Intercept) -1.10179159
#> s__Bacteroides_plebeius/s__Bacteroides_dorei -0.02659576
#> s__Eubacterium_eligens/s__Bacteroides_cellulosilyticus 0.04043200
#> s__Collinsella_aerofaciens/s__Bifidobacterium_pseudocatenulatum -0.04886115
#> s__Streptococcus_salivarius/s__Actinomyces_sp_HMSC035G02 0.09235497
#> s__Erysipelatoclostridium_ramosum/s__Parvimonas_micra 0.04784129
#> s__Holdemania_filiformis/s__Eubacterium_sp_CAG_38 -0.03798158
#> s__Holdemania_filiformis/s__Proteobacteria_bacterium_CAG_139 -0.03610014
#> s__Anaerotruncus_colihominis/s__Clostridium_aldenense -0.07966635
#> s__Lachnospira_pectinoschiza/s__Ruminococcus_torques 0.03247124
#> s__Roseburia_faecis/s__Actinomyces_sp_HMSC035G02 0.05287895
#> s__Coprobacter_fastidiosus/s__Parvimonas_micra 0.06757447
#> Std. Error
#> (Intercept) 0.26205522
#> s__Bacteroides_plebeius/s__Bacteroides_dorei 0.01099163
#> s__Eubacterium_eligens/s__Bacteroides_cellulosilyticus 0.01385914
#> s__Collinsella_aerofaciens/s__Bifidobacterium_pseudocatenulatum 0.01507927
#> s__Streptococcus_salivarius/s__Actinomyces_sp_HMSC035G02 0.02625928
#> s__Erysipelatoclostridium_ramosum/s__Parvimonas_micra 0.01965003
#> s__Holdemania_filiformis/s__Eubacterium_sp_CAG_38 0.01622550
#> s__Holdemania_filiformis/s__Proteobacteria_bacterium_CAG_139 0.01837978
#> s__Anaerotruncus_colihominis/s__Clostridium_aldenense 0.02173633
#> s__Lachnospira_pectinoschiza/s__Ruminococcus_torques 0.01300914
#> s__Roseburia_faecis/s__Actinomyces_sp_HMSC035G02 0.01743096
#> s__Coprobacter_fastidiosus/s__Parvimonas_micra 0.01818976
#> z value
#> (Intercept) -4.204425
#> s__Bacteroides_plebeius/s__Bacteroides_dorei -2.419638
#> s__Eubacterium_eligens/s__Bacteroides_cellulosilyticus 2.917353
#> s__Collinsella_aerofaciens/s__Bifidobacterium_pseudocatenulatum -3.240286
#> s__Streptococcus_salivarius/s__Actinomyces_sp_HMSC035G02 3.517042
#> s__Erysipelatoclostridium_ramosum/s__Parvimonas_micra 2.434668
#> s__Holdemania_filiformis/s__Eubacterium_sp_CAG_38 -2.340857
#> s__Holdemania_filiformis/s__Proteobacteria_bacterium_CAG_139 -1.964123
#> s__Anaerotruncus_colihominis/s__Clostridium_aldenense -3.665124
#> s__Lachnospira_pectinoschiza/s__Ruminococcus_torques 2.496033
#> s__Roseburia_faecis/s__Actinomyces_sp_HMSC035G02 3.033621
#> s__Coprobacter_fastidiosus/s__Parvimonas_micra 3.714973
#> Pr(>|z|)
#> (Intercept) 2.617464e-05
#> s__Bacteroides_plebeius/s__Bacteroides_dorei 1.553598e-02
#> s__Eubacterium_eligens/s__Bacteroides_cellulosilyticus 3.530165e-03
#> s__Collinsella_aerofaciens/s__Bifidobacterium_pseudocatenulatum 1.194097e-03
#> s__Streptococcus_salivarius/s__Actinomyces_sp_HMSC035G02 4.363849e-04
#> s__Erysipelatoclostridium_ramosum/s__Parvimonas_micra 1.490548e-02
#> s__Holdemania_filiformis/s__Eubacterium_sp_CAG_38 1.923954e-02
#> s__Holdemania_filiformis/s__Proteobacteria_bacterium_CAG_139 4.951583e-02
#> s__Anaerotruncus_colihominis/s__Clostridium_aldenense 2.472185e-04
#> s__Lachnospira_pectinoschiza/s__Ruminococcus_torques 1.255908e-02
#> s__Roseburia_faecis/s__Actinomyces_sp_HMSC035G02 2.416375e-03
#> s__Coprobacter_fastidiosus/s__Parvimonas_micra 2.032250e-04
Generating taxa selection probabilities
It is encouraged to run k-fold cross-validation for several times to
account for the random fold splits. FLORAL provides
mcv.FLORAL functions to repeat cross-validations for
mcv times and on ncore cores. The output
summarizes taxa selection probabilities, average coefficients based on
\(\lambda_\min\) and \(\lambda_{\text{1se}}\). Interpretable plots
can be created if plot = TRUE is specified.
mcv.fit <- mcv.FLORAL(mcv=2,
ncore=1,
x = x,
y = y,
family = "binomial",
ncv = 3,
progress=TRUE)
#> Warning in mcv.FLORAL(mcv = 2, ncore = 1, x = x, y = y, family = "binomial", :
#> Using 1 core for computation.
#> Random 3-fold cross-validation: 1
#> Random 3-fold cross-validation: 2
#Other options are also available
#mcv.fit$p_min_ratio
#mcv.fit$p_1se
#mcv.fit$p_1se_ratio
Elastic net
Beyond lasso model, FLORAL also supports elastic net
models by specifying the tuning parameter a between 0 and
1. Lasso penalty will be used when a=1 while ridge penalty
will be used when a=0.
The a.FLORAL function can help investigate the
prediction performance for different choices of a and
return a plot of the corresponding prediction metric trajectories
against the choice of \(\lambda\).
a.fit <- a.FLORAL(a = c(0.1,1),
ncore = 1,
x = x,
y = y,
family = "binomial",
ncv = 3,
progress=TRUE)
#> Warning in a.FLORAL(a = c(0.1, 1), ncore = 1, x = x, y = y, family =
#> "binomial", : Using 1 core for computation.
#> Running for a = 0.1
#> Running for a = 1
a.fit
