Packages Used

library(tidyverse)
library(dplyr)
library(phyloseq)
library(corncob)
library(lme4)
library(decontam)
library(ggplot2)
library(vegan)
library(egg)

Tissue and Swab Analysis

Data organization

ts.meta <- read.csv("R_files/ts_meta_samples.csv", row.names = 1)
ts.counts <- read.csv("R_files/ts_counts_decontam.csv", row.names = 1)
ts.tax <- read.csv("R_files/ts_tax_decontam.csv", row.names = 1)

ts.counts <- as.data.frame(t(ts.counts))
str(ts.meta)
## 'data.frame':    46 obs. of  7 variables:
##  $ exp.code  : chr  "GH" "GH" "GH" "GH" ...
##  $ form      : chr  "field.swab" "field.swab" "field.swab" "field.swab" ...
##  $ type      : chr  "cloacal.swab" "cloacal.swab" "cloacal.swab" "cloacal.swab" ...
##  $ toe.clip  : int  101 2 4 13 17 18 23 24 4 23 ...
##  $ note      : chr  "" "" "" "" ...
##  $ location  : chr  "" "" "North Fork" "" ...
##  $ is.control: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
ts.meta$toe.clip <- as.factor(ts.meta$toe.clip)
ts.meta$form <- as.factor(ts.meta$form)
ts.meta$type <- as.factor(ts.meta$type)

ts.meta <- ts.meta[order(rownames(ts.meta)),]
ts.counts <- ts.counts[order(rownames(ts.counts)),]

str(ts.meta)
## 'data.frame':    46 obs. of  7 variables:
##  $ exp.code  : chr  "GH" "GH" "GH" "GH" ...
##  $ form      : Factor w/ 3 levels "field.swab","lab.swab",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ type      : Factor w/ 5 levels "cloaca","cloacal.swab",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ toe.clip  : Factor w/ 9 levels "2","4","13","17",..: 9 1 2 3 4 5 7 8 2 7 ...
##  $ note      : chr  "" "" "" "" ...
##  $ location  : chr  "" "" "North Fork" "" ...
##  $ is.control: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...

Remove unused field swabs

unused.field <- c("GH04","GH13","GH17","GH18","GH23","GH24")
ts.meta <- ts.meta[!(rownames(ts.meta) %in% (unused.field)),]
ts.counts <- ts.counts[rownames(ts.counts) %in% rownames(ts.meta),]
ts.counts <- ts.counts[,colSums(ts.counts) > 27]
ts.tax <- ts.tax[rownames(ts.tax) %in% colnames(ts.counts),]

ts.counts.ps <- otu_table(ts.counts, taxa_are_rows = FALSE)
colnames(ts.counts.ps) <- colnames(ts.counts)

ts.tax.ps <- tax_table(ts.tax)
taxa_names(ts.tax.ps) <- taxa_names(ts.counts.ps)
colnames(ts.tax.ps) <- c("Kingdom","Phylum","Class","Order","Family","Genus")
ts.meta.ps <- sample_data(ts.meta)


ts.ps <- phyloseq(ts.counts.ps, ts.tax.ps, ts.meta.ps)

ts.fam.ps <- tax_glom(ts.ps, "Family", NArm = FALSE)
ts.ps <- tax_glom(ts.ps, "Family", NArm = FALSE)
richness <- estimate_richness(ts.ps, measures = "Observed")
ts.ps@otu_table <- transform_sample_counts(ts.ps@otu_table, 
                                                  function(x) log10(x +1))
shannon<- estimate_richness(ts.ps, measures = "Shannon")
ts.meta <- cbind(ts.meta, shannon,richness)
ts.meta.ps <- sample_data(ts.meta)
ts.ps@sam_data <- ts.meta.ps

Remove low read tissue sample, and single sample from toe clip 21

ts.meta <- ts.meta[rownames(ts.meta) != "TS20",]
ts.counts <- ts.counts[rownames(ts.counts) %in% rownames(ts.meta),]
ts.ps <- subset_samples(ts.ps, sample_names(ts.ps) != "TS20")
ts.fam.ps <- subset_samples(ts.fam.ps, sample_names(ts.fam.ps) != "TS20")

ts.meta <- ts.meta[rownames(ts.meta) != "TS06",]
ts.counts <- ts.counts[rownames(ts.counts) %in% rownames(ts.meta),]
ts.ps <- subset_samples(ts.ps, sample_names(ts.ps) != "TS06")
ts.fam.ps <- subset_samples(ts.fam.ps, sample_names(ts.fam.ps) != "TS06")

Alpha Diversity

Shannon Diversity Index

Check assumptions

ts.aov <- aov(Shannon ~ type, data = ts.meta)

plot(density(ts.aov$residuals))

qqnorm(ts.aov$residuals)
qqline(ts.aov$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))

plot(ts.aov$residuals~ts.aov$fitted.values)
lines(lowess(ts.aov$fitted.values,ts.aov$residuals), col="blue")

ok

Anova with random effect

ts.aov <- aov(Shannon ~ type + Error(toe.clip), data = ts.meta)

summary(ts.aov)
## 
## Error: toe.clip
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       1  0.717  0.7173   1.326  0.293
## Residuals  6  3.246  0.5410               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       4  3.524  0.8811   1.352  0.278
## Residuals 26 16.948  0.6518

Richness Check assumptions

ts.aov.rich <- aov(Observed ~ type, data = ts.meta)

plot(density(ts.aov.rich$residuals))

qqnorm(ts.aov.rich$residuals)
qqline(ts.aov.rich$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))

plot(ts.aov.rich$residuals~ts.aov.rich$fitted.values)
lines(lowess(ts.aov.rich$fitted.values,ts.aov.rich$residuals), col="blue")

ok anova with random effects

ts.aov.rich <- aov(Observed ~ type + Error(toe.clip), data = ts.meta)

summary(ts.aov.rich)
## 
## Error: toe.clip
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       1  61.66   61.66   2.691  0.152
## Residuals  6 137.49   22.91               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       4  183.8   45.96   1.306  0.294
## Residuals 26  915.2   35.20

Beta Diversity

All sample types

set.seed(1)
ts.mds <- ordinate(ts.ps, method = "NMDS", distance = "bray", k = 4, trymax =  100, maxit = 100)
## Run 0 stress 0.07661727 
## Run 1 stress 0.07661648 
## ... New best solution
## ... Procrustes: rmse 0.001068949  max resid 0.005313108 
## ... Similar to previous best
## Run 2 stress 0.07662647 
## ... Procrustes: rmse 0.001280951  max resid 0.006205534 
## ... Similar to previous best
## Run 3 stress 0.07662087 
## ... Procrustes: rmse 0.001341989  max resid 0.006382977 
## ... Similar to previous best
## Run 4 stress 0.07661761 
## ... Procrustes: rmse 0.001140824  max resid 0.005682364 
## ... Similar to previous best
## Run 5 stress 0.0766154 
## ... New best solution
## ... Procrustes: rmse 0.00063105  max resid 0.002917568 
## ... Similar to previous best
## Run 6 stress 0.07722737 
## Run 7 stress 0.07708785 
## ... Procrustes: rmse 0.007175722  max resid 0.03068457 
## Run 8 stress 0.07661663 
## ... Procrustes: rmse 0.0006593384  max resid 0.003079713 
## ... Similar to previous best
## Run 9 stress 0.07702336 
## ... Procrustes: rmse 0.00969492  max resid 0.04832688 
## Run 10 stress 0.07927013 
## Run 11 stress 0.08391244 
## Run 12 stress 0.07661722 
## ... Procrustes: rmse 0.0008870524  max resid 0.004527009 
## ... Similar to previous best
## Run 13 stress 0.2433405 
## Run 14 stress 0.07661555 
## ... Procrustes: rmse 0.0004104975  max resid 0.002085157 
## ... Similar to previous best
## Run 15 stress 0.07663305 
## ... Procrustes: rmse 0.001785753  max resid 0.006247106 
## ... Similar to previous best
## Run 16 stress 0.0766199 
## ... Procrustes: rmse 0.001395348  max resid 0.007308458 
## ... Similar to previous best
## Run 17 stress 0.07776051 
## Run 18 stress 0.07661548 
## ... Procrustes: rmse 0.0003039463  max resid 0.001399574 
## ... Similar to previous best
## Run 19 stress 0.07662036 
## ... Procrustes: rmse 0.00144499  max resid 0.007517465 
## ... Similar to previous best
## Run 20 stress 0.07663836 
## ... Procrustes: rmse 0.001826092  max resid 0.00881322 
## ... Similar to previous best
## *** Solution reached

Dispersion

ts.dist <- vegdist(ts.ps@otu_table, method = "bray")

anova(betadisper(ts.dist, ts.ps@sam_data[["type"]]))
## Analysis of Variance Table
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq F value  Pr(>F)  
## Groups     4 0.085574 0.0213936  3.5424 0.01641 *
## Residuals 33 0.199299 0.0060394                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Composition

ts.perma <- adonis(ts.dist ~ type, data = ts.meta)
ts.perma
## 
## Call:
## adonis(formula = ts.dist ~ type, data = ts.meta) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## type       4    1.8634 0.46586   2.392 0.22477  0.003 **
## Residuals 33    6.4271 0.19476         0.77523          
## Total     37    8.2905                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Split “low” and “high” groups

ps.low <- subset_samples(ts.ps, type != "up.intestine")
ps.low <- subset_samples(ps.low, type != "oviduct")

ps.high <- subset_samples(ts.ps, type != "low.intestine")
ps.high <- subset_samples(ps.high, type != "cloaca")
ps.high <- subset_samples(ps.high, type != "cloacal.swab")

meta.low <- ts.meta[ts.meta$type != "up.intestine" & ts.meta$type != "oviduct",]
meta.high <- ts.meta[ts.meta$type == "up.intestine" | ts.meta$type == "oviduct",]
set.seed(1)
low.mds <- ordinate(ps.low, method = "NMDS", distance = "bray", k = 3, trymax = 100, maxit = 100 )
## Run 0 stress 0.0541883 
## Run 1 stress 0.05654577 
## Run 2 stress 0.0529863 
## ... New best solution
## ... Procrustes: rmse 0.0340552  max resid 0.1203461 
## Run 3 stress 0.05417035 
## Run 4 stress 0.05354653 
## Run 5 stress 0.05430399 
## Run 6 stress 0.05426918 
## Run 7 stress 0.05724595 
## Run 8 stress 0.05655692 
## Run 9 stress 0.05733475 
## Run 10 stress 0.05709366 
## Run 11 stress 0.05329056 
## ... Procrustes: rmse 0.01467522  max resid 0.04421703 
## Run 12 stress 0.05460982 
## Run 13 stress 0.05323083 
## ... Procrustes: rmse 0.06451203  max resid 0.1834151 
## Run 14 stress 0.05417019 
## Run 15 stress 0.05320061 
## ... Procrustes: rmse 0.009223069  max resid 0.0277863 
## Run 16 stress 0.05342451 
## ... Procrustes: rmse 0.01993474  max resid 0.08239169 
## Run 17 stress 0.05650636 
## Run 18 stress 0.05737242 
## Run 19 stress 0.05664218 
## Run 20 stress 0.05635574 
## Run 21 stress 0.0574207 
## Run 22 stress 0.05531926 
## Run 23 stress 0.05318169 
## ... Procrustes: rmse 0.008747289  max resid 0.02724001 
## Run 24 stress 0.0532311 
## ... Procrustes: rmse 0.06450666  max resid 0.1834468 
## Run 25 stress 0.05418718 
## Run 26 stress 0.05329126 
## ... Procrustes: rmse 0.0144107  max resid 0.04257044 
## Run 27 stress 0.05418626 
## Run 28 stress 0.0532358 
## ... Procrustes: rmse 0.0643848  max resid 0.1835616 
## Run 29 stress 0.05417051 
## Run 30 stress 0.05417478 
## Run 31 stress 0.05417292 
## Run 32 stress 0.05319398 
## ... Procrustes: rmse 0.009036223  max resid 0.02860076 
## Run 33 stress 0.05298823 
## ... Procrustes: rmse 0.0004777605  max resid 0.001126366 
## ... Similar to previous best
## *** Solution reached

Lower region dispersion

low.dist <- vegdist(ps.low@otu_table, method = "bray")

anova(betadisper(low.dist, ps.low@sam_data[["type"]]))
## Analysis of Variance Table
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq F value  Pr(>F)  
## Groups     2 0.048931 0.0244656  3.4755 0.05065 .
## Residuals 20 0.140788 0.0070394                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Lower region composition

set.seed(1)
low.perma <- adonis(low.dist ~ type, data = meta.low)
low.perma
## 
## Call:
## adonis(formula = low.dist ~ type, data = meta.low) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)
## type       2    0.4870 0.24352  1.4129 0.1238  0.223
## Residuals 20    3.4471 0.17235         0.8762       
## Total     22    3.9341                 1.0000

Upper region dispersion

high.dist <- vegdist(ps.high@otu_table, method = "bray")

anova(betadisper(high.dist, ps.high@sam_data[["type"]]))
## Analysis of Variance Table
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq F value Pr(>F)
## Groups     1 0.000099 0.0000995  0.0236 0.8803
## Residuals 13 0.054842 0.0042186

Upper region composition

set.seed(1)
high.perma <- adonis(high.dist ~ type, data = meta.high)
high.perma
## 
## Call:
## adonis(formula = high.dist ~ type, data = meta.high) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## type       1    0.3889 0.38888  1.6965 0.11543  0.058 .
## Residuals 13    2.9800 0.22923         0.88457         
## Total     14    3.3689                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Differential Abundance

ts.meta$cc <- 0

for (i in 1:38){
  if(ts.meta$type[i] == "cloacal.swab") {
    ts.meta$cc[i] <- "A_clo.swab" 
  }else if(ts.meta$type[i] == "cloaca") {
    ts.meta$cc[i] <- "B_clo" 
  }else if(ts.meta$type[i] == "low.intestine") {
    ts.meta$cc[i] <- "C_LI" 
  }else if(ts.meta$type[i] == "up.intestine") {
    ts.meta$cc[i] <- "D_UI" 
  }else if(ts.meta$type[i] == "oviduct") {
    ts.meta$cc[i] <- "E_ovi"
  } 
}
ts.meta.ps <- sample_data(ts.meta)
ts.fam.ps@sam_data <- ts.meta.ps

Explore whole community

set.seed(1)
cc.ts <- differentialTest(formula = ~ type,
                                 phi.formula = ~ type,
                                 formula_null = ~ 1,
                                 phi.formula_null = ~ type,
                                 test = "LRT", boot = FALSE,
                                 data = ts.fam.ps,
                                 fdr_cutoff = 0.1)
cc.ts$significant_taxa
## [1] "ASV_1"  "ASV_11" "ASV_14"
otu_to_taxonomy(cc.ts$significant_taxa, ts.fam.ps, level = "Family")
##                ASV_1               ASV_11               ASV_14 
## "Enterobacteriaceae"    "Ruminococcaceae"     "Bacteroidaceae"

check models of significant taxa

enter.cc <- bbdml(formula = ASV_1 ~ cc, 
                  phi.formula = ~ cc,
                  data = ts.fam.ps)
enter.cc
## 
## Call:
## bbdml(formula = ASV_1 ~ cc, phi.formula = ~cc, data = ts.fam.ps)
## 
## 
## Coefficients associated with abundance:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   1.4569     0.4785   3.045  0.00503 **
## ccB_clo       0.2147     0.7984   0.269  0.78998   
## ccC_LI       -0.5141     0.6945  -0.740  0.46534   
## ccD_UI       -2.0470     0.5977  -3.425  0.00191 **
## ccE_ovi      -1.2724     0.5863  -2.170  0.03864 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Coefficients associated with dispersion:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.8432     0.5113  -1.649    0.110
## ccB_clo       0.5980     0.8059   0.742    0.464
## ccC_LI        0.7841     0.6938   1.130    0.268
## ccD_UI       -0.2747     0.6799  -0.404    0.689
## ccE_ovi      -0.5087     0.7020  -0.725    0.475
## 
## 
## Log-likelihood: -319.6
rumin.cc <- bbdml(formula = ASV_11 ~ type, 
                  phi.formula = ~ type,
                  data = ts.fam.ps)
rumin.cc
## 
## Call:
## bbdml(formula = ASV_11 ~ type, phi.formula = ~type, data = ts.fam.ps)
## 
## 
## Coefficients associated with abundance:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -10.273      1.029  -9.983  1.0e-10 ***
## typecloacal.swab     5.572      1.410   3.952 0.000478 ***
## typelow.intestine    5.958      1.372   4.341 0.000167 ***
## typeoviduct          4.260      1.711   2.490 0.018966 *  
## typeup.intestine     9.237      1.164   7.939  1.2e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Coefficients associated with dispersion:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -8.560      1.485  -5.765 3.45e-06 ***
## typecloacal.swab     5.917      1.849   3.200  0.00341 ** 
## typelow.intestine    6.198      1.798   3.446  0.00181 ** 
## typeoviduct          5.135      2.168   2.368  0.02501 *  
## typeup.intestine     8.795      1.565   5.619 5.13e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Log-likelihood: -152.71
bacter.cc <- bbdml(formula = ASV_14 ~ type, 
                  phi.formula = ~ type,
                  data = ts.fam.ps)
bacter.cc
## 
## Call:
## bbdml(formula = ASV_14 ~ type, phi.formula = ~type, data = ts.fam.ps)
## 
## 
## Coefficients associated with abundance:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -7.5984     0.6028 -12.605 4.63e-13 ***
## typecloacal.swab    4.5184     1.0210   4.425 0.000133 ***
## typelow.intestine   4.7356     0.9647   4.909 3.56e-05 ***
## typeoviduct         4.3698     1.2160   3.594 0.001234 ** 
## typeup.intestine    5.1195     1.0102   5.068 2.30e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Coefficients associated with dispersion:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -6.7178     0.8413  -7.985 1.07e-08 ***
## typecloacal.swab    5.5436     1.2293   4.510 0.000106 ***
## typelow.intestine   5.5681     1.1681   4.767 5.25e-05 ***
## typeoviduct         5.8806     1.4161   4.153 0.000279 ***
## typeup.intestine    6.3898     1.1922   5.360 1.04e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Log-likelihood: -191.4

Feces and Cloacal swab analysis

Data organization

pfs.meta <- read.csv("pfs/R_files/pfs_meta_samples.csv", row.names = 1)
pfs.counts <- read.csv("pfs/R_files/pfs_counts_decontam.csv", row.names = 1)
pfs.tax <- read.csv("pfs/R_files/pfs_tax_decontam.csv", row.names = 1)

pfs.counts <- as.data.frame(t(pfs.counts))
str(pfs.meta)
## 'data.frame':    60 obs. of  6 variables:
##  $ date      : chr  "" "" "31-May" "31-May" ...
##  $ form      : chr  "insect" "insect" "cloacal.swab" "cloacal.swab" ...
##  $ type      : chr  "insect" "insect" "baseline" "baseline" ...
##  $ sex       : chr  "" "" "F" "M" ...
##  $ toe.clip  : int  NA NA 2404 2411 2402 4012 2953 2953 2420 2420 ...
##  $ is.control: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
pfs.meta$toe.clip <- as.factor(pfs.meta$toe.clip)
pfs.meta$form <- as.factor(pfs.meta$form)
pfs.meta$type <- as.factor(pfs.meta$type)

pfs.meta <- pfs.meta[order(rownames(pfs.meta)),]
pfs.counts <- pfs.counts[order(rownames(pfs.counts)),]
pfs.meta$type <- factor(pfs.meta$type, levels = c("baseline","T0","6.hour","24.hour","insect","feces"))
str(pfs.meta)
## 'data.frame':    60 obs. of  6 variables:
##  $ date      : chr  "" "" "31-May" "31-May" ...
##  $ form      : Factor w/ 3 levels "cloacal.swab",..: 3 3 1 1 1 1 1 2 1 2 ...
##  $ type      : Factor w/ 6 levels "baseline","T0",..: 5 5 1 1 1 1 2 6 2 6 ...
##  $ sex       : chr  "" "" "F" "M" ...
##  $ toe.clip  : Factor w/ 12 levels "2402","2403",..: NA NA 3 6 1 12 10 10 9 9 ...
##  $ is.control: logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
pfs.counts <- pfs.counts[,colSums(pfs.counts) > 27]
pfs.tax <- pfs.tax[rownames(pfs.tax) %in% colnames(pfs.counts),]

pfs.counts.ps <- otu_table(pfs.counts, taxa_are_rows = FALSE)
colnames(pfs.counts.ps) <- colnames(pfs.counts)

pfs.tax.ps <- tax_table(pfs.tax)
taxa_names(pfs.tax.ps) <- taxa_names(pfs.counts.ps)
colnames(pfs.tax.ps) <- c("Kingdom","Phylum","Class","Order","Family","Genus")
pfs.meta.ps <- sample_data(pfs.meta)


pfs.ps <- phyloseq(pfs.counts.ps, pfs.tax.ps, pfs.meta.ps)

pfs.fam.ps <- tax_glom(pfs.ps, "Family", NArm = FALSE)
pfs.ps <- tax_glom(pfs.ps, "Family", NArm = FALSE)
richness <- estimate_richness(pfs.ps, measures = "Observed")
pfs.ps@otu_table <- transform_sample_counts(pfs.ps@otu_table, 
                                                  function(x) log10(x +1))
shannon<- estimate_richness(pfs.ps, measures = "Shannon")
pfs.meta <- cbind(pfs.meta, shannon,richness)
pfs.meta.ps <- sample_data(pfs.meta)
pfs.ps@sam_data <- pfs.meta.ps
total.reads <- rowSums(pfs.counts)
total.reads
## Cricket1 Cricket2     PFS1    PFS10    PFS11    PFS12    PFS13    PFS14 
##    13614    19032     5443     5999     3486     6616     9135    27256 
##    PFS15    PFS16    PFS17    PFS18    PFS19     PFS2    PFS20    PFS21 
##     4492    89754     2370    94822     1750     7241    14721     9011 
##    PFS22    PFS23    PFS24    PFS25    PFS26    PFS27    PFS28    PFS29 
##    66590    36361    31919    10420     3713    10671    11278    15611 
##     PFS3    PFS30    PFS31    PFS32    PFS34    PFS35    PFS37    PFS39 
##     4933    13343    16878    24037    11193    20672    13504    17077 
##     PFS4    PFS41    PFS43    PFS45    PFS46    PFS47    PFS48    PFS49 
##     5952    28980    14183    24097    24727    10539   103934    12868 
##     PFS5    PFS50    PFS51    PFS52    PFS53    PFS55    PFS57    PFS58 
##      730     5336     6558    23374    12621    10478    13576    96615 
##     PFS6    PFS65    PFS66    PFS67    PFS68     PFS7    PFS70    PFS71 
##      980    24223     3989    12488    12842     1038    10037    21660 
##    PFS72    PFS73     PFS8     PFS9 
##    18244    23135     1697    16636
pfs.meta <- cbind(pfs.meta, total.reads)

Alpha Diversity

Shannon Diversity

Feces compared to baseline swabs

liz.meta <- pfs.meta[pfs.meta$form != "insect",]
liz.form.meta <- liz.meta[liz.meta$type == "baseline" | liz.meta$type == "feces",]

Assumption Check

liz.form.aov <- aov(Shannon ~ form * sex, data = liz.form.meta)
plot(density(liz.form.aov$residuals))

qqnorm(liz.form.aov$residuals)
qqline(liz.form.aov$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))

plot(liz.form.aov$residuals~liz.form.aov$fitted.values)
lines(lowess(liz.form.aov$fitted.values,liz.form.aov$residuals), col="blue")

Model

liz.form.aov <- aov(Shannon ~ form * sex + Error(toe.clip), data = liz.form.meta)
summary(liz.form.aov)
## 
## Error: toe.clip
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## form       1 2.3240  2.3240  11.987 0.00854 **
## sex        1 0.1221  0.1221   0.630 0.45032   
## form:sex   1 0.0972  0.0972   0.501 0.49906   
## Residuals  8 1.5510  0.1939                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## form       1  9.521   9.521  93.183 7.08e-05 ***
## form:sex   1  0.040   0.040   0.388    0.556    
## Residuals  6  0.613   0.102                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Difference in time point swabs

swab.meta <- liz.meta[liz.meta$form != "feces",]
swab.aov <- aov(Shannon ~ type *sex, data = swab.meta)

plot(density(swab.aov$residuals))

qqnorm(swab.aov$residuals)
qqline(swab.aov$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))

plot(swab.aov$residuals~swab.aov$fitted.values)
lines(lowess(swab.aov$fitted.values,swab.aov$residuals), col="blue")

Model

swab.aov <- aov(Shannon ~ type * sex + Error(toe.clip), data = swab.meta)
summary(swab.aov)
## 
## Error: toe.clip
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       3  2.898  0.9661   0.694  0.589
## sex        1  0.112  0.1116   0.080  0.787
## type:sex   1  0.154  0.1542   0.111  0.751
## Residuals  6  8.353  1.3921               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       3  0.882  0.2939   0.517  0.674
## type:sex   3  3.334  1.1115   1.954  0.141
## Residuals 32 18.198  0.5687

Richness

Feces v. Swabs assumption check

liz.form.aov.rich <- aov(Observed ~ form * sex, data = liz.form.meta)
plot(density(liz.form.aov.rich$residuals))

qqnorm(liz.form.aov.rich$residuals)
qqline(liz.form.aov.rich$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))

plot(liz.form.aov.rich$residuals~liz.form.aov.rich$fitted.values)
lines(lowess(liz.form.aov.rich$fitted.values,liz.form.aov.rich$residuals), col="blue")

model

liz.form.aov.rich <- aov(Observed ~ form * sex + Error(toe.clip), data = liz.form.meta)
summary(liz.form.aov.rich)
## 
## Error: toe.clip
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## form       1  418.6   418.6  27.705 0.000761 ***
## sex        1   12.8    12.8   0.847 0.384272    
## form:sex   1   12.0    12.0   0.795 0.398605    
## Residuals  8  120.9    15.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## form       1 1958.1  1958.1 194.590 8.46e-06 ***
## form:sex   1    5.1     5.1   0.503    0.505    
## Residuals  6   60.4    10.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Swab time point assumption check

swab.aov.rich <- aov(Observed ~ type * sex, data = swab.meta)

plot(density(swab.aov.rich$residuals))

qqnorm(swab.aov.rich$residuals)
qqline(swab.aov.rich$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))

plot(swab.aov.rich$residuals~swab.aov.rich$fitted.values)
lines(lowess(swab.aov.rich$fitted.values,swab.aov.rich$residuals), col="blue")

Model

swab.aov.rich <- aov(Observed ~ type * sex + Error(toe.clip), data = swab.meta)
summary(swab.aov.rich)
## 
## Error: toe.clip
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       3  348.4  116.12   1.015  0.449
## sex        1    0.3    0.26   0.002  0.964
## type:sex   1    9.6    9.61   0.084  0.782
## Residuals  6  686.3  114.39               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## type       3  216.4   72.15   1.486  0.237
## type:sex   3  286.4   95.46   1.966  0.139
## Residuals 32 1553.4   48.54

Beta Diversity

Subset data

liz.ps <- subset_samples(pfs.ps, type != "insect")
swab.ps <- subset_samples(liz.ps, type != "feces")
form.ps <- subset_samples(liz.ps, type != "24.hour")
form.ps <- subset_samples(form.ps, type != "6.hour")
form.ps <- subset_samples(form.ps, type != "T0")

Feces and baseline swabs

Dispersion

form.dist <- vegdist(form.ps@otu_table, method = "bray")

anova(betadisper(form.dist, form.ps@sam_data[["form"]]))
## Analysis of Variance Table
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq F value   Pr(>F)   
## Groups     1 0.036876 0.036876  9.9354 0.005514 **
## Residuals 18 0.066808 0.003712                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(betadisper(form.dist, form.ps@sam_data[["sex"]]))
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq F value Pr(>F)
## Groups     1 0.00011 0.0001098  0.0138 0.9078
## Residuals 18 0.14330 0.0079611

Composition

set.seed(1)
form.perma <- adonis(form.dist ~ form * sex, data = liz.form.meta)
form.perma
## 
## Call:
## adonis(formula = form.dist ~ form * sex, data = liz.form.meta) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## form       1    2.3533 2.35334 29.7237 0.62601  0.001 ***
## sex        1    0.0425 0.04251  0.5369 0.01131  0.656    
## form:sex   1    0.0966 0.09662  1.2203 0.02570  0.255    
## Residuals 16    1.2668 0.07917         0.33698           
## Total     19    3.7592                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Defecation time point

Dispersion

swab.dist <- vegdist(swab.ps@otu_table, method = "bray")

anova(betadisper(swab.dist, swab.ps@sam_data[["type"]]))
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value   Pr(>F)   
## Groups     3 0.10083 0.033611  4.5118 0.007418 **
## Residuals 46 0.34268 0.007450                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(betadisper(swab.dist, swab.ps@sam_data[["sex"]]))
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq F value Pr(>F)
## Groups     1 0.00648 0.0064795  0.8744 0.3544
## Residuals 48 0.35568 0.0074100

Composition

set.seed(1)
swab.perma <- adonis(swab.dist ~ type, data = swab.meta)
swab.perma
## 
## Call:
## adonis(formula = swab.dist ~ type, data = swab.meta) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## type       3    0.5729 0.19096  1.2234 0.07389  0.278
## Residuals 46    7.1800 0.15609         0.92611       
## Total     49    7.7529                 1.00000

Differential Abundance

Feces compared to baseline cloacal swabs

Explore whole community

lizfam.ps <- subset_samples(pfs.fam.ps, type != "insect")
form.fam.ps <- subset_samples(lizfam.ps, type != "24.hour")
form.fam.ps <- subset_samples(form.fam.ps, type != "6.hour")
form.fam.ps <- subset_samples(form.fam.ps, type != "T0")
set.seed(1)
cc.form<- differentialTest(formula = ~ form,
                                 phi.formula = ~ form,
                                 formula_null = ~ 1,
                                 phi.formula_null = ~ form,
                                 test = "LRT", boot = FALSE,
                                 data = form.fam.ps,
                                 fdr_cutoff = 0.05)
cc.form$significant_taxa
## [1] "ASV_1" "ASV_6"
otu_to_taxonomy(cc.form$significant_taxa, form.fam.ps, level = "Family")
##                ASV_1                ASV_6 
## "Enterobacteriaceae"    "Lachnospiraceae"

Check ASVs

enter.cc <- bbdml(formula = ASV_1 ~ form, 
                  phi.formula = ~ form,
                  data = form.fam.ps)
enter.cc
## 
## Call:
## bbdml(formula = ASV_1 ~ form, phi.formula = ~form, data = form.fam.ps)
## 
## 
## Coefficients associated with abundance:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.7874     0.3784   4.724  0.00023 ***
## formfeces    -4.1316     0.5808  -7.114 2.46e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Coefficients associated with dispersion:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -1.2769     0.4382  -2.914   0.0101 *
## formfeces    -0.6660     0.7120  -0.935   0.3635  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Log-likelihood: -159.98
lachno.cc <- bbdml(formula = ASV_6 ~ form, 
                  phi.formula = ~ form,
                  data = form.fam.ps)
lachno.cc
## 
## Call:
## bbdml(formula = ASV_6 ~ form, phi.formula = ~form, data = form.fam.ps)
## 
## 
## Coefficients associated with abundance:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.0606     0.4741  -8.566 2.26e-07 ***
## formfeces     3.5943     0.5160   6.965 3.19e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Coefficients associated with dispersion:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.0733     0.5821  -5.279 7.48e-05 ***
## formfeces     0.6191     0.7552   0.820    0.424    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Log-likelihood: -138.59