library(tidyverse)
library(dplyr)
library(phyloseq)
library(corncob)
library(lme4)
library(decontam)
library(ggplot2)
library(vegan)
library(egg)
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")
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
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
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
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)
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
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
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