Practicing Corncob

The purpose of this lab is to introduce you all to a useful tool for analyzing and visualizing relative abundance. This is the analysis that came from the paper you all read. The tutorial below is written by the author of the paper. Use the following link to access the tutorial and the code:

https://cran.r-project.org/web/packages/corncob/vignettes/corncob-intro-no-phyloseq.html

As you go through the tutorial, copy and paste code below to perform the tasks. Additionally, I have a few questions you need to answer.

Install Packages and view data

library(corncob)
library(magrittr)
data(soil_phylo_sample)
data(soil_phylo_otu)
data(soil_phylo_taxa)

head(soil_phylo_sample)

##      Plants DayAmdmt Amdmt ID Day
## S009      1       01     1  D   0
## S204      1       21     1  D   2
## S112      0       11     1  B   1
## S247      0       22     2  F   2
## S026      0       00     0  A   0
## S023      1       00     0  C   0

soil_phylo_otu[1:5, 1:5]

##         S009 S204 S112 S247 S026
## OTU.43   350   74  300   70   43
## OTU.2   1796 4204 1752  695  945
## OTU.187  280  709  426  100  139
## OTU.150   33  151   18   13   28
## OTU.91     0    0  184    0    0

soil_phylo_taxa[1:3, ]

##         Kingdom    Phylum           Class                 Order             
## OTU.43  "Bacteria" "Nitrospirae"    "Nitrospira"          "Nitrospirales"   
## OTU.2   "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"     
## OTU.187 "Bacteria" "Acidobacteria"  "Acidobacteriia"      "Acidobacteriales"
##         Family              Genus            Species
## OTU.43  "Nitrospiraceae"    "Nitrospira"     ""     
## OTU.2   "Bradyrhizobiaceae" "Bradyrhizobium" ""     
## OTU.187 "Koribacteraceae"   ""               ""

Fit a Model

# Subset data and group samples by Phylum
data(soil_phylum_small_sample)
sample_data <- soil_phylum_small_sample
data(soil_phylum_small_otu)
data <- soil_phylum_small_otu

#Make dataframe to include sample data with Proteobacteria counts
pro_data <- cbind(sample_data, 
                  W = unlist(data["Proteobacteria", ]),
                  M = colSums(data))

corncob <- bbdml(formula = cbind(W, M - W) ~ 1,
             phi.formula = ~ 1,
             data = pro_data)

Interpret Model

plot(corncob, B = 50)

#Plot with 95% prediction intervals 
plot(corncob, total = TRUE, B = 50)

#With color
plot(corncob, total = TRUE, color = "DayAmdmt", B = 50)

plot(corncob, color = "DayAmdmt", B = 50)

Add Covariates

#Model data with expected relative abundance and variability counts as a covariate by modifying formula
corncob_da <- bbdml(formula = cbind(W, M - W) ~ DayAmdmt,
             phi.formula = ~ DayAmdmt,
             data = pro_data)
plot(corncob_da, color = "DayAmdmt", total = TRUE, B = 50)

plot(corncob_da, color = "DayAmdmt", B = 50)

Model Selection

#Liklihood ratio test to select final model and test the null: covariates model = noncovariates model
lrtest(mod_null = corncob, mod = corncob_da)

## [1] 4.550571e-05

#low p-value

Parameter Interpretation

summary(corncob_da)

## 
## Call:
## bbdml(formula = cbind(W, M - W) ~ DayAmdmt, phi.formula = ~DayAmdmt, 
##     data = pro_data)
## 
## 
## Coefficients associated with abundance:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.44595    0.03604 -12.375 7.18e-13 ***
## DayAmdmt21  -0.16791    0.04067  -4.129 0.000297 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Coefficients associated with dispersion:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.3077     0.3537 -15.008 6.44e-15 ***
## DayAmdmt21   -1.3518     0.5029  -2.688    0.012 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Log-likelihood: -286.53

Analysis for Multiple Taxa

#To test all the taxa in our data to look for differential abundance or variability. 

#View a list of differentially abundant taxa
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)
da_analysis

## Object of class differentialTest 
## 
## $p: p-values 
## $p_fdr: FDR-adjusted p-values 
## $significant_taxa: taxa names of the statistically significant taxa 
## $significant_models: model summaries of the statistically significant taxa 
## $all_models: all model summaries 
## $restrictions_DA: covariates tested for differential abundance 
## $restrictions_DV: covariates tested for differential variability 
## $discriminant_taxa_DA: taxa for which at least one covariate associated with the abundance was perfectly discriminant 
## $discriminant_taxa_DV: taxa for which at least one covariate associated with the dispersion was perfectly discriminant 
## 
## plot( ) to see a plot of tested coefficients from significant taxa

da_analysis$significant_taxa

##  [1] "Proteobacteria"   "Gemmatimonadetes" "Bacteroidetes"    "Cyanobacteria"   
##  [5] "Firmicutes"       "Planctomycetes"   "Armatimonadetes"  "Spirochaetes"    
##  [9] "Elusimicrobia"    "BRC1"             "OP3"              "FBP"             
## [13] "Chlorobi"         "TM6"

#View a list of differentially variable taxa
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

## [1] "Acidobacteria" "Cyanobacteria" "Spirochaetes"  "Elusimicrobia"
## [5] "FBP"

#examine p-values
da_analysis$p[1:5]

##    Acidobacteria   Proteobacteria Gemmatimonadetes   Actinobacteria 
##     6.509418e-01     3.642734e-05     3.270448e-13     3.703095e-01 
##         [Thermi] 
##     1.031419e-01

#subset p-values
da_analysis$p_fdr[1:5]

##    Acidobacteria   Proteobacteria Gemmatimonadetes   Actinobacteria 
##     7.811302e-01     1.457094e-04     3.924537e-12     5.172744e-01 
##         [Thermi] 
##     2.062838e-01

#view model coefficients 
plot(da_analysis)

## `height` was translated to `width`.

df <- plot(da_analysis, data_only = TRUE)

#Make custom plots 
# we can easily remove special characters used in our formatting steps
df <- df %>%
  dplyr::mutate(variable = gsub("\nDifferential Abundance", "",
                                variable, fixed = TRUE))

head(df)

##            x        xmin        xmax             taxa   variable
## 1 -0.1679129 -0.24761731 -0.08820844   Proteobacteria DayAmdmt21
## 2  0.3119623  0.22800588  0.39591870 Gemmatimonadetes DayAmdmt21
## 3  0.2306765  0.09148981  0.36986312    Bacteroidetes DayAmdmt21
## 4  1.7090454  1.33732670  2.08076404    Cyanobacteria DayAmdmt21
## 5 -0.3087530 -0.49665466 -0.12085125       Firmicutes DayAmdmt21
## 6  0.2013158  0.11895313  0.28367838   Planctomycetes DayAmdmt21

#View unfit taxa not in the model
which(is.na(da_analysis$p)) %>% names

## [1] "GN04"   "GN02"   "MVP-21"

#re-examine data to see wht it did not fit the model
data["GN04", ]

##      S204 S112 S134 S207 S202 S139 S122 S212 S117 S104 S214 S109 S217 S229 S132
## GN04    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
##      S209 S227 S107 S237 S224 S127 S137 S114 S124 S119 S219 S232 S129 S102 S234
## GN04    0    0    0    0    0    0    0    0    0    0    0    0    0    1    0
##      S222 S239
## GN04    0    0

Answering scientific questions

#Test for differential abundance across Day without controlling for anything else and plot
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)

## `height` was translated to `width`.

#Control for the effect of Day on dispresion
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)

## `height` was translated to `width`.

#Joint test 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)

## `height` was translated to `width`.

QUESTIONS: 1. How do you interpret the first relative abundance plot? (Black and White one) - The relative abundance of Proteobacteria varies across samples. Each point represents an observed relative abundance, each line represetns that 95% confidence interval.

Why do we add covariates and what does it do to the relative abundance plot (the red and blue one at the end of the section, compare to previous red and blue plot)

Adding covariates allows the model to account for variables that may explain differences in relative abundance that are not explicitly being tested for. After adding covariates, the model can account for differences between the group and produces a better fit with lower expected relative abundance and variability.

In the final section of the tutorial “Examples of Answering Scientific Questions”, why is it useful/helpful to make plots like these? How would you interpret these plots?

It is helpful to have a visual to summarize the results of hypothesis tests across multiple taxa. The plots also make it easier to identify differences between taxa. For these plots the coefficients indicate higher (positive), lower (negative), or no (near zero) change in estimated response relative to the reference group. This helps to identify patterns in different microbial taxa and help researchers decide which to further investigate.