Nothing
Introduction to corncob, no phyloseq
In corncob: Count Regression for Correlated Observations with the Beta-Binomial
\addtolength{\headheight}{-.025\textheight}
\thispagestyle{fancyplain}
\rhead{\includegraphics[height=.1\textheight]{../docs/logo.png}}
\renewcommand{\headrulewidth}{0pt}
Vignette Information
This is a version of the corncob-intro.Rmd vignette that does not rely on the package phyloseq. To see how to use corncob with phyloseq objects, or to see an additional analysis with an IBD microbiome dataset, check out the original corncob-intro.Rmd vignette.
We thank Dr. Thea Whitman for kindly providing us with the example data set we use for this vignette. You can read more about this data in Whitman, Thea, et al. "Dynamics of microbial community composition and soil organic carbon mineralization in soil following addition of pyrogenic and fresh organic matter." The ISME Journal 10.12 (2016): 2918.
Introduction
Effectively modeling microbial relative abundance poses a number of statistical challenges, including:
\begin{itemize}
\item different sequencing depth,
\item excessive zeros from unobserved taxa,
\item high variability of empirical relative abundances (overdispersion),
\item within-taxon correlation,
\item hypothesis testing with categorical and continuous covariates.
\end{itemize}
knitr::opts_chunk$set(fig.width=8, fig.height=4)
Here, we introduce corncob, an individual taxon regression model that uses abundance tables and sample data. corncob is able to model differential abundance and differential variability, and addresses each of the challenges presented above.
Install corncob using:
remotes::install_github("statdivlab/corncob")
To begin, we load our example data set we load our example data set as three different data frames
library(corncob)
library(magrittr)
data(soil_phylo_sample)
data(soil_phylo_otu)
data(soil_phylo_taxa)
We can see that we have 5 sample variables. Let's look at the first few observations.
head(soil_phylo_sample)
Our covariates are as follows:
\begin{itemize}
\item Plants: Indicator of whether plants are in the soil for this sample.
\item Amdmt: Categorical variable representing one of three soil additives; none, biochar, and freshbiomass, respectively.
\item ID: Categorical variable representing different plots of soil.
\item Day: Categorical variable representing one of three days of measurement; day 1, day 21, and day 81, respectively.
\item DayAmdmt: Categorical variable combining the Day and Amdmt variables into a single variable.
\end{itemize}
Additionally, we have a table that gives us counts of the phyla by sample. Let's take a look at the first 5 phyla and the first 5 samples.
soil_phylo_otu[1:5, 1:5]
Finally, we have a taxonomy table with 7 taxonomic ranks.
soil_phylo_taxa[1:3, ]
Fitting a Model
First, let's subset our samples to only include those with the DayAmdmt covariate equal to 11 or 21 and then collapse the samples to the phylum level. We have already done this and saved the resulting sample data frame and otu data frame.
data(soil_phylum_small_sample)
sample_data <- soil_phylum_small_sample
data(soil_phylum_small_otu)
data <- soil_phylum_small_otu
Note that collapsing the samples is not necessary, and this model can work at any taxonomic rank. However, we will later be fitting a model to every taxa. We can see that by agglomerating taxa to the phylum level, we have gone from from 7770 to 40 taxa. Thus we collapse in order to increase the speed for the purposes of this tutorial.
Now we fit our model. We will demonstrate with Proteobacteria. We will need to make a data frame that includes our sample data as well as the Proteobacteria counts and sequencing depths for each sample.
pro_data <- cbind(sample_data,
W = unlist(data["Proteobacteria", ]),
M = colSums(data))
For now, we will not include any covariates, so we use ~ 1 as our model formula responses.
corncob <- bbdml(formula = cbind(W, M - W) ~ 1,
phi.formula = ~ 1,
data = pro_data)
\newpage
Interpreting a Model
First, let's plot the data with our model fit on the relative abundance scale. To do this, we simply type:
plot(corncob, B = 50)
The points represent the relative abundances. The bars represent the 95\% prediction intervals for the observed relative abundance by sample.
The parameter B determines the number of bootstrap simulations used to approximate the prediction intervals. For purposes of this tutorial, we use a small value B = 50 for computational purposes, but recommend a higher setting for more accurate intervals, such as the default B = 1000.
Now let's look at the same plot, but on the counts scale with 95\% prediction intervals (since counts is not a parameter). To do this, we add the option total = TRUE to our plotting code.
plot(corncob, total = TRUE, B = 50)
Finally, let's color the plot by the DayAmdmt covariate. To do this, we add the option color = "DayAmdmt" to our plotting code.
plot(corncob, total = TRUE, color = "DayAmdmt", B = 50)
plot(corncob, color = "DayAmdmt", B = 50)
Notice that this plot also reorders our samples so that groups appear together so that they are easier to compare.
We can observe on this plot that it might be of interest to distinguish between the two groups with covariates. The average empirical relative abundance for the samples with DayAmdmt = 21 tends to be lower and less variable than the samples with DayAmdmt = 11.
Adding covariates
Let's try modeling the expected relative abundance and the variability of the counts with DayAmdmt as a covariate. We do this by modifying formula and phi.formula as:
corncob_da <- bbdml(formula = cbind(W, M - W) ~ DayAmdmt,
phi.formula = ~ DayAmdmt,
data = pro_data)
Let's also plot this data on both the total count and relative abundance scales.
plot(corncob_da, color = "DayAmdmt", total = TRUE, B = 50)
plot(corncob_da, color = "DayAmdmt", B = 50)
Visually, the model with covariates seems to provide a much better fit to the data, but how can we compare the two models statistically?
Model Selection
Let's use a likelihood ratio test to select our final model for this taxon. We want to test the null hypothesis that the likelihood of the model with covariates is equal to the likelihood of the model without covariates. To do this test, we use:
lrtest(mod_null = corncob, mod = corncob_da)
We obtain a p-value much smaller than a cut-off of 0.05. Therefore we conclude that there is a statistically significant difference in the likelihood of the two models. Thus, we probably want to use the model with covariates for this taxon.
Parameter Interpretation
Now that we have chosen our model, let's interpret our model output. To see a summary of the model, type:
summary(corncob_da)
This output will look familiar if you have done regression analysis in R in the past. Covariates associated with the expected relative abundance are presented separately from covariates associated with the variance of the counts are preceded by.
From this model summary, we can see that the DayAmdmt21 abundance coefficient is negative and statistically significant. This suggests that this taxon is differentially-abundant across DayAmdmt, and that samples with DayAmdmt = 21 are expected to have a lower relative abundance. This matches what we saw from the observed abundances.
We can also see that the DayAmdmt21 dispersion coefficient is negative and statistically significant. This suggests that this taxon is differentially-variable across DayAmdmt, and that samples with DayAmdmt = 21 are expected to have a lower variability. This matches what we saw from the observed abundances.
Analysis for Multiple Taxa
What if we want to test all the taxa in our data to see if they are differentially-abundant or differentially-variable? We use the differentialTest function. It will perform the above tests on all taxa, and it will control the false discovery rate to account for multiple comparisons.
Next, we use the differentialTest command. We specify the covariates of our model using formula and phi.formula as before, except we no longer include the response term because we are testing multiple taxa. We also specify which covariates we want to test for by removing them in the formula_null and phi.formula_null arguments.
The difference between the formulas and the null version of the formulas are the variables that we test. We will go through several examples, starting with a test for differential abundance across the DayAmdmt coefficient.
We set fdr_cutoff to be our controlled false discovery rate.
Here, we will pass in our entire count table for all phyla as the data argument (setting taxa_are_rows to TRUE), and the sample data frame as the sample_data argument.
set.seed(1)
da_analysis <- differentialTest(formula = ~ DayAmdmt,
phi.formula = ~ DayAmdmt,
formula_null = ~ 1,
phi.formula_null = ~ DayAmdmt,
test = "Wald", boot = FALSE,
data = data,
sample_data = sample_data,
taxa_are_rows = TRUE,
fdr_cutoff = 0.05)
We can see the output of the function by calling it:
da_analysis
We can see a list of differentially-abundant taxa using:
da_analysis$significant_taxa
In this case, we identified r length(da_analysis$significant_taxa) taxa that are differentially-abundant across DayAmdmt (out of the 39 taxa tested).
We can see a list of differentially-variable taxa using:
set.seed(1)
dv_analysis <- differentialTest(formula = ~ DayAmdmt,
phi.formula = ~ DayAmdmt,
formula_null = ~ DayAmdmt,
phi.formula_null = ~ 1,
test = "LRT", boot = FALSE,
data = data,
sample_data = sample_data,
taxa_are_rows = TRUE,
fdr_cutoff = 0.05)
dv_analysis$significant_taxa
In this case, we identified r length(dv_analysis$significant_taxa) taxa that are differentially-variable across DayAmdmt (out of the 40 taxa tested).
We can examine a subset of the p-values of our tests using:
da_analysis$p[1:5]
We can examine a subset of the p-values after controlling for the false discovery rate using:
da_analysis$p_fdr[1:5]
where the values are now adjusted to control the false discovery rate at 0.05.
We can also plot the model coefficients of our results:
plot(da_analysis)
Here, we can see that for Bacteria_Armatimonadetes, the effect of DayAmdmt21 is positive compared to the baseline (in this case, DayAmdmt11).
Finally, we can see a list of any taxa for which we were not able to fit a model using:
which(is.na(da_analysis$p)) %>% names
In this case, we weren't able to fit GN04 automatically. It's worthwhile to investigate the OTU individually if this is the case.
It may be that the model is overparameterized because there aren't enough observations, or it may just be that the initializations were invalid for that taxa and it needs to be re-evaluated with new initializations.
Let's first try examining the data.
data["GN04", ]
We see that the observed counts of OTU is zero in all samples except for S102, where we observed a single count. We should be skeptical of any statistical model fit on a single observed count!
corncob is stable, but if you notice any issues, please log them on Github to help us help you!
Examples of Answering Scientific Questions
We will now walk through several scientific questions of interest and show how they can be answered using hypothesis testing with corncob. Note that Day and Amdmt are both factor covariates with levels 0, 1, and 2.
Note that some of these are rather strange tests, and shown for demonstration of the flexibility of the model only. Normally, when testing for differential variability across a covariate, we recommend always controlling for the effect of that covariate on the abundance. We first demonstrate examples with the small version of the soil dataset (note, this is a slightly different subset of the data from that used in corncob-intro.Rmd).
Testing for differential abundance across Day, without controlling for anything else:
ex1 <- differentialTest(formula = ~ Day,
phi.formula = ~ 1,
formula_null = ~ 1,
phi.formula_null = ~ 1,
data = data,
taxa_are_rows = TRUE,
sample_data = sample_data,
test = "Wald", boot = FALSE,
fdr_cutoff = 0.05)
plot(ex1)
Testing for differential abundance across Day, controlling for the effect of Day on dispersion:
ex2 <- differentialTest(formula = ~ Day,
phi.formula = ~ Day,
formula_null = ~ 1,
phi.formula_null = ~ Day,
data = data,
taxa_are_rows = TRUE,
sample_data = sample_data,
test = "Wald", boot = FALSE,
fdr_cutoff = 0.05)
plot(ex2)
Jointly testing for differential abundance and differential variability across Day:
ex3 <- differentialTest(formula = ~ Day,
phi.formula = ~ Day,
formula_null = ~ 1,
phi.formula_null = ~ 1,
data = data,
taxa_are_rows = TRUE,
sample_data = sample_data,
test = "Wald", boot = FALSE,
fdr_cutoff = 0.05)
plot(ex3)
More examples of answering scientific questions with a larger subset of this soil data set and with an IBD data set can be seen in corncob-intro.Rmd.
Try the corncob package in your browser
Any scripts or data that you put into this service are public.
corncob documentation built on May 29, 2024, 7:15 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
\addtolength{\headheight}{-.025\textheight} \thispagestyle{fancyplain} \rhead{\includegraphics[height=.1\textheight]{../docs/logo.png}} \renewcommand{\headrulewidth}{0pt}
Vignette Information
This is a version of the corncob-intro.Rmd vignette that does not rely on the package phyloseq. To see how to use corncob with phyloseq objects, or to see an additional analysis with an IBD microbiome dataset, check out the original corncob-intro.Rmd vignette.
We thank Dr. Thea Whitman for kindly providing us with the example data set we use for this vignette. You can read more about this data in Whitman, Thea, et al. "Dynamics of microbial community composition and soil organic carbon mineralization in soil following addition of pyrogenic and fresh organic matter." The ISME Journal 10.12 (2016): 2918.
Introduction
Effectively modeling microbial relative abundance poses a number of statistical challenges, including: \begin{itemize} \item different sequencing depth, \item excessive zeros from unobserved taxa, \item high variability of empirical relative abundances (overdispersion), \item within-taxon correlation, \item hypothesis testing with categorical and continuous covariates. \end{itemize}
knitr::opts_chunk$set(fig.width=8, fig.height=4)
Here, we introduce corncob, an individual taxon regression model that uses abundance tables and sample data. corncob is able to model differential abundance and differential variability, and addresses each of the challenges presented above.
Install corncob using:
remotes::install_github("statdivlab/corncob")
To begin, we load our example data set we load our example data set as three different data frames
library(corncob) library(magrittr) data(soil_phylo_sample) data(soil_phylo_otu) data(soil_phylo_taxa)
We can see that we have 5 sample variables. Let's look at the first few observations.
head(soil_phylo_sample)
Our covariates are as follows:
\begin{itemize}
\item Plants: Indicator of whether plants are in the soil for this sample.
\item Amdmt: Categorical variable representing one of three soil additives; none, biochar, and freshbiomass, respectively.
\item ID: Categorical variable representing different plots of soil.
\item Day: Categorical variable representing one of three days of measurement; day 1, day 21, and day 81, respectively.
\item DayAmdmt: Categorical variable combining the Day and Amdmt variables into a single variable.
\end{itemize}
Additionally, we have a table that gives us counts of the phyla by sample. Let's take a look at the first 5 phyla and the first 5 samples.
soil_phylo_otu[1:5, 1:5]
Finally, we have a taxonomy table with 7 taxonomic ranks.
soil_phylo_taxa[1:3, ]
Fitting a Model
First, let's subset our samples to only include those with the DayAmdmt covariate equal to 11 or 21 and then collapse the samples to the phylum level. We have already done this and saved the resulting sample data frame and otu data frame.
data(soil_phylum_small_sample) sample_data <- soil_phylum_small_sample data(soil_phylum_small_otu) data <- soil_phylum_small_otu
Note that collapsing the samples is not necessary, and this model can work at any taxonomic rank. However, we will later be fitting a model to every taxa. We can see that by agglomerating taxa to the phylum level, we have gone from from 7770 to 40 taxa. Thus we collapse in order to increase the speed for the purposes of this tutorial.
Now we fit our model. We will demonstrate with Proteobacteria. We will need to make a data frame that includes our sample data as well as the Proteobacteria counts and sequencing depths for each sample.
pro_data <- cbind(sample_data, W = unlist(data["Proteobacteria", ]), M = colSums(data))
For now, we will not include any covariates, so we use ~ 1 as our model formula responses.
corncob <- bbdml(formula = cbind(W, M - W) ~ 1, phi.formula = ~ 1, data = pro_data)
\newpage
Interpreting a Model
First, let's plot the data with our model fit on the relative abundance scale. To do this, we simply type:
plot(corncob, B = 50)
The points represent the relative abundances. The bars represent the 95\% prediction intervals for the observed relative abundance by sample.
The parameter B determines the number of bootstrap simulations used to approximate the prediction intervals. For purposes of this tutorial, we use a small value B = 50 for computational purposes, but recommend a higher setting for more accurate intervals, such as the default B = 1000.
Now let's look at the same plot, but on the counts scale with 95\% prediction intervals (since counts is not a parameter). To do this, we add the option total = TRUE to our plotting code.
plot(corncob, total = TRUE, B = 50)
Finally, let's color the plot by the DayAmdmt covariate. To do this, we add the option color = "DayAmdmt" to our plotting code.
plot(corncob, total = TRUE, color = "DayAmdmt", B = 50)
plot(corncob, color = "DayAmdmt", B = 50)
Notice that this plot also reorders our samples so that groups appear together so that they are easier to compare.
We can observe on this plot that it might be of interest to distinguish between the two groups with covariates. The average empirical relative abundance for the samples with DayAmdmt = 21 tends to be lower and less variable than the samples with DayAmdmt = 11.
Adding covariates
Let's try modeling the expected relative abundance and the variability of the counts with DayAmdmt as a covariate. We do this by modifying formula and phi.formula as:
corncob_da <- bbdml(formula = cbind(W, M - W) ~ DayAmdmt, phi.formula = ~ DayAmdmt, data = pro_data)
Let's also plot this data on both the total count and relative abundance scales.
plot(corncob_da, color = "DayAmdmt", total = TRUE, B = 50)
plot(corncob_da, color = "DayAmdmt", B = 50)
Visually, the model with covariates seems to provide a much better fit to the data, but how can we compare the two models statistically?
Model Selection
Let's use a likelihood ratio test to select our final model for this taxon. We want to test the null hypothesis that the likelihood of the model with covariates is equal to the likelihood of the model without covariates. To do this test, we use:
lrtest(mod_null = corncob, mod = corncob_da)
We obtain a p-value much smaller than a cut-off of 0.05. Therefore we conclude that there is a statistically significant difference in the likelihood of the two models. Thus, we probably want to use the model with covariates for this taxon.
Parameter Interpretation
Now that we have chosen our model, let's interpret our model output. To see a summary of the model, type:
summary(corncob_da)
This output will look familiar if you have done regression analysis in R in the past. Covariates associated with the expected relative abundance are presented separately from covariates associated with the variance of the counts are preceded by.
From this model summary, we can see that the DayAmdmt21 abundance coefficient is negative and statistically significant. This suggests that this taxon is differentially-abundant across DayAmdmt, and that samples with DayAmdmt = 21 are expected to have a lower relative abundance. This matches what we saw from the observed abundances.
We can also see that the DayAmdmt21 dispersion coefficient is negative and statistically significant. This suggests that this taxon is differentially-variable across DayAmdmt, and that samples with DayAmdmt = 21 are expected to have a lower variability. This matches what we saw from the observed abundances.
Analysis for Multiple Taxa
What if we want to test all the taxa in our data to see if they are differentially-abundant or differentially-variable? We use the differentialTest function. It will perform the above tests on all taxa, and it will control the false discovery rate to account for multiple comparisons.
Next, we use the differentialTest command. We specify the covariates of our model using formula and phi.formula as before, except we no longer include the response term because we are testing multiple taxa. We also specify which covariates we want to test for by removing them in the formula_null and phi.formula_null arguments.
The difference between the formulas and the null version of the formulas are the variables that we test. We will go through several examples, starting with a test for differential abundance across the DayAmdmt coefficient.
We set fdr_cutoff to be our controlled false discovery rate.
Here, we will pass in our entire count table for all phyla as the data argument (setting taxa_are_rows to TRUE), and the sample data frame as the sample_data argument.
set.seed(1) da_analysis <- differentialTest(formula = ~ DayAmdmt, phi.formula = ~ DayAmdmt, formula_null = ~ 1, phi.formula_null = ~ DayAmdmt, test = "Wald", boot = FALSE, data = data, sample_data = sample_data, taxa_are_rows = TRUE, fdr_cutoff = 0.05)
We can see the output of the function by calling it:
da_analysis
We can see a list of differentially-abundant taxa using:
da_analysis$significant_taxa
In this case, we identified r length(da_analysis$significant_taxa) taxa that are differentially-abundant across DayAmdmt (out of the 39 taxa tested).
We can see a list of differentially-variable taxa using:
set.seed(1) dv_analysis <- differentialTest(formula = ~ DayAmdmt, phi.formula = ~ DayAmdmt, formula_null = ~ DayAmdmt, phi.formula_null = ~ 1, test = "LRT", boot = FALSE, data = data, sample_data = sample_data, taxa_are_rows = TRUE, fdr_cutoff = 0.05) dv_analysis$significant_taxa
In this case, we identified r length(dv_analysis$significant_taxa) taxa that are differentially-variable across DayAmdmt (out of the 40 taxa tested).
We can examine a subset of the p-values of our tests using:
da_analysis$p[1:5]
We can examine a subset of the p-values after controlling for the false discovery rate using:
da_analysis$p_fdr[1:5]
where the values are now adjusted to control the false discovery rate at 0.05.
We can also plot the model coefficients of our results:
plot(da_analysis)
Here, we can see that for Bacteria_Armatimonadetes, the effect of DayAmdmt21 is positive compared to the baseline (in this case, DayAmdmt11).
Finally, we can see a list of any taxa for which we were not able to fit a model using:
which(is.na(da_analysis$p)) %>% names
In this case, we weren't able to fit GN04 automatically. It's worthwhile to investigate the OTU individually if this is the case.
It may be that the model is overparameterized because there aren't enough observations, or it may just be that the initializations were invalid for that taxa and it needs to be re-evaluated with new initializations.
Let's first try examining the data.
data["GN04", ]
We see that the observed counts of OTU is zero in all samples except for S102, where we observed a single count. We should be skeptical of any statistical model fit on a single observed count!
corncob is stable, but if you notice any issues, please log them on Github to help us help you!
Examples of Answering Scientific Questions
We will now walk through several scientific questions of interest and show how they can be answered using hypothesis testing with corncob. Note that Day and Amdmt are both factor covariates with levels 0, 1, and 2.
Note that some of these are rather strange tests, and shown for demonstration of the flexibility of the model only. Normally, when testing for differential variability across a covariate, we recommend always controlling for the effect of that covariate on the abundance. We first demonstrate examples with the small version of the soil dataset (note, this is a slightly different subset of the data from that used in corncob-intro.Rmd).
Testing for differential abundance across Day, without controlling for anything else:
ex1 <- differentialTest(formula = ~ Day, phi.formula = ~ 1, formula_null = ~ 1, phi.formula_null = ~ 1, data = data, taxa_are_rows = TRUE, sample_data = sample_data, test = "Wald", boot = FALSE, fdr_cutoff = 0.05) plot(ex1)
Testing for differential abundance across Day, controlling for the effect of Day on dispersion:
ex2 <- differentialTest(formula = ~ Day, phi.formula = ~ Day, formula_null = ~ 1, phi.formula_null = ~ Day, data = data, taxa_are_rows = TRUE, sample_data = sample_data, test = "Wald", boot = FALSE, fdr_cutoff = 0.05) plot(ex2)
Jointly testing for differential abundance and differential variability across Day:
ex3 <- differentialTest(formula = ~ Day, phi.formula = ~ Day, formula_null = ~ 1, phi.formula_null = ~ 1, data = data, taxa_are_rows = TRUE, sample_data = sample_data, test = "Wald", boot = FALSE, fdr_cutoff = 0.05) plot(ex3)
More examples of answering scientific questions with a larger subset of this soil data set and with an IBD data set can be seen in corncob-intro.Rmd.
Try the corncob package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.