Review Analyses
Marie Bunker
2022-06-21
Packages
library(tidyverse)
library(dplyr)
library(phyloseq)
library(corncob)
library(lme4)
library(ggplot2)
library(vegan)
library(msa)
library(DECIPHER)
library(phangorn)
library(ape)
library(picante)
library(ggtext)
library(pairwiseAdonis)
library(MicrobiotaProcess)
library(VennDiagram)
library(ggVennDiagram)
library(RVenn)
library(ggvenn)
library(microbiome)Prep
counts <- read.csv("R_files/counts.csv",row.names = 1)
meta <- read.csv("R_files/meta.csv", row.names = 1)
tax <- read.csv("R_files/tax.csv", row.names = 1)
tree <- read.tree("R_files/tree")
counts <- as.data.frame(t(counts))offseason.females <- c("F178","F193","F38","F39","MR10","F154")meta$toe.clip <- as.factor(meta$toe.clip)
meta$sex <- as.factor(meta$sex)
meta$RS <- as.factor(meta$RS)
meta$season <- as.factor(meta$season)
meta$year <- as.factor(meta$year)
meta$half.month <- as.factor(meta$half.month)
meta <- meta[order(rownames(meta)),]
counts <- counts[order(rownames(counts)),]
meta$season <- factor(meta$season, levels = c("pre","breeding","post"))
meta$RS <- factor(meta$RS, levels =c("pre","V","G","PostOvi"))
str(meta)## 'data.frame': 258 obs. of 9 variables:
## $ toe.clip : Factor w/ 209 levels "","1011","1020",..: 60 61 25 67 199 202 69 70 124 25 ...
## $ sex : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
## $ svl : num 61 58 70 65 66 65 64 67.5 63 69 ...
## $ RS : Factor w/ 4 levels "pre","V","G",..: 2 2 2 2 2 2 2 2 1 1 ...
## $ season : Factor w/ 3 levels "pre","breeding",..: 2 2 2 2 2 2 2 2 1 1 ...
## $ year : Factor w/ 3 levels "2017","2018",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ half.month: Factor w/ 12 levels "April_Early",..: 10 10 10 10 10 10 10 10 1 1 ...
## $ month : chr "May" "May" "May" "May" ...
## $ bmass : num 6.9 5.4 9.3 7.8 8.3 9.1 7.5 9 6.6 7.7 ...counts <- counts[,colSums(counts) > 10]
tax <- tax[rownames(tax) %in% colnames(counts),]
counts.ps <- otu_table(counts, taxa_are_rows = FALSE)
colnames(counts.ps) <- colnames(counts)
tax.ps <- phyloseq::tax_table(tax)
taxa_names(tax.ps) <- taxa_names(counts.ps)
colnames(tax.ps) <- c("Kingdom","Phylum","Class","Order","Family","Genus")
meta.ps <- sample_data(meta)
ps <- phyloseq(counts.ps, tax.ps, meta.ps, tree)
ps.real <- psrichness <- estimate_richness(ps, measures = "Observed")ps@otu_table <- transform_sample_counts(ps@otu_table,
function(x) log10(x +1))shannon<- estimate_richness(ps, measures = "Shannon")
meta <- cbind(meta, shannon,richness)
meta.ps <- sample_data(meta)
ps@sam_data <- meta.pspick_new_outgroup <- function(tree.unrooted){
require("magrittr")
require("data.table")
require("ape") # ape::Ntip
# tablify parts of tree that we need.
treeDT <-
cbind(
data.table(tree.unrooted$edge),
data.table(length = tree.unrooted$edge.length)
)[1:Ntip(tree.unrooted)] %>%
cbind(data.table(id = tree.unrooted$tip.label))
# Take the longest terminal branch as outgroup
new.outgroup <- treeDT[which.max(length)]$id
return(new.outgroup)
}
new.outgroup = pick_new_outgroup(ps@phy_tree)
ps@phy_tree <- root(ps@phy_tree, outgroup=new.outgroup, resolve.root=TRUE)
is.rooted(ps@phy_tree)## [1] TRUEfam.ps <- tax_glom(ps.real, "Family", NArm = FALSE)
phy.ps <- tax_glom(ps.real, "Phylum", NArm = FALSE)pd<- pd(ps@otu_table, ps@phy_tree)
meta$PD <- pd$PDmeta %>%
count(sex, season)## sex season n
## 1 F pre 30
## 2 F breeding 66
## 3 F post 43
## 4 M pre 23
## 5 M breeding 64
## 6 M post 32meta %>%
group_by(sex) %>%
summarise(mean.svl = mean(svl))## # A tibble: 2 × 2
## sex mean.svl
## <fct> <dbl>
## 1 F 63.9
## 2 M 57.7samples.f <-
meta[meta$sex == "F",] %>%
group_by(RS) %>%
summarize(mean.svl = mean(svl))
samples.f## # A tibble: 4 × 2
## RS mean.svl
## <fct> <dbl>
## 1 pre 64.1
## 2 V 64.8
## 3 G 62.9
## 4 PostOvi 63.8meta[meta$sex == "F",] %>%
count(RS)## RS n
## 1 pre 31
## 2 V 32
## 3 G 30
## 4 PostOvi 46rc1 <- rarecurve(ps.real@otu_table,label = FALSE)
rc.zoom <- rarecurve(ps.real@otu_table,label = FALSE, xlim=c(0, 2000))ps.full <- ps
ps <- subset_samples(ps, !sample_names(ps) %in% offseason.females)
meta.full <- meta
meta <- meta[!rownames(meta) %in% offseason.females,]
counts.full <- counts
counts <- counts[!rownames(counts) %in% offseason.females,]Year Analysis
Alpha
year.meta <- meta[meta$half.month == "June_Late" | meta$half.month == "May_Late",]year.aov.shan <- aov(Shannon ~ year * sex * svl, data = year.meta)
plot(density(year.aov.shan$residuals))
qqnorm(year.aov.shan$residuals)
qqline(year.aov.shan$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(year.aov.shan$residuals~year.aov.shan$fitted.values)
lines(lowess(year.aov.shan$fitted.values,year.aov.shan$residuals), col="blue")
year.aov.ob <- aov(log10(Observed) ~ year * sex * svl, data = year.meta)
plot(density(year.aov.ob$residuals))
qqnorm(year.aov.ob$residuals)
qqline(year.aov.ob$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(year.aov.ob$residuals~year.aov.ob$fitted.values)
lines(lowess(year.aov.ob$fitted.values,year.aov.ob$residuals), col="blue")
year.aov.pd <- aov(log10(PD) ~ year * sex * svl, data = year.meta)
plot(density(year.aov.pd$residuals))
qqnorm(year.aov.pd$residuals)
qqline(year.aov.pd$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(year.aov.pd$residuals~year.aov.pd$fitted.values)
lines(lowess(year.aov.pd$fitted.values,year.aov.pd$residuals), col="blue")
aov.year.shan <- aov(Shannon ~ year * sex * svl + Error(toe.clip), data = year.meta)
aov.year.ob <- aov(log10(Observed) ~ year * sex * svl + Error(toe.clip), data = year.meta)
aov.year.PD <- aov(log10(PD) ~ year * sex * svl + Error(toe.clip), data = year.meta)
summary(aov.year.shan)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.709 0.3543 2.002 0.1439
## sex 1 0.494 0.4939 2.790 0.1000 .
## svl 1 0.046 0.0459 0.259 0.6124
## year:sex 2 0.948 0.4739 2.678 0.0768 .
## year:svl 2 0.251 0.1254 0.708 0.4964
## sex:svl 1 0.009 0.0089 0.050 0.8238
## year:sex:svl 2 0.670 0.3352 1.894 0.1592
## Residuals 61 10.797 0.1770
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq
## year 2 0.10303 0.05151
## svl 1 0.12540 0.12540
## year:sex 1 0.01413 0.01413summary(aov.year.ob)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.1405 0.07026 2.086 0.133
## sex 1 0.0912 0.09119 2.708 0.105
## svl 1 0.0109 0.01089 0.323 0.572
## year:sex 2 0.1792 0.08959 2.660 0.078 .
## year:svl 2 0.0472 0.02361 0.701 0.500
## sex:svl 1 0.0017 0.00167 0.050 0.824
## year:sex:svl 2 0.1313 0.06564 1.949 0.151
## Residuals 61 2.0543 0.03368
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq
## year 2 0.019169 0.009584
## svl 1 0.024701 0.024701
## year:sex 1 0.003439 0.003439summary(aov.year.PD)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.391 0.19555 2.268 0.112
## sex 1 0.106 0.10551 1.224 0.273
## svl 1 0.188 0.18762 2.176 0.145
## year:sex 2 0.253 0.12670 1.470 0.238
## year:svl 2 0.002 0.00104 0.012 0.988
## sex:svl 1 0.009 0.00924 0.107 0.745
## year:sex:svl 2 0.353 0.17661 2.049 0.138
## Residuals 61 5.259 0.08622
##
## Error: Within
## Df Sum Sq Mean Sq
## year 2 0.13054 0.06527
## svl 1 0.22951 0.22951
## year:sex 1 0.05218 0.05218aov.year.shan <- aov(Shannon ~ year + sex + svl + Error(toe.clip), data = year.meta)
aov.year.ob <- aov(log10(Observed) ~ year + sex + svl + Error(toe.clip), data = year.meta)
aov.year.PD <- aov(log10(PD) ~ year + sex + svl + Error(toe.clip), data = year.meta)
summary(aov.year.shan)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.709 0.3543 1.901 0.157
## sex 1 0.494 0.4939 2.649 0.108
## svl 1 0.046 0.0459 0.246 0.621
## Residuals 68 12.675 0.1864
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.10303 0.05151 3.645 0.347
## svl 1 0.12540 0.12540 8.872 0.206
## Residuals 1 0.01413 0.01413summary(aov.year.ob)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.1405 0.07026 1.979 0.146
## sex 1 0.0912 0.09119 2.569 0.114
## svl 1 0.0109 0.01089 0.307 0.582
## Residuals 68 2.4136 0.03549
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.019169 0.009584 2.787 0.390
## svl 1 0.024701 0.024701 7.183 0.227
## Residuals 1 0.003439 0.003439summary(aov.year.PD)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.391 0.19555 2.263 0.112
## sex 1 0.106 0.10551 1.221 0.273
## svl 1 0.188 0.18762 2.171 0.145
## Residuals 68 5.877 0.08643
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## year 2 0.13054 0.06527 1.251 0.534
## svl 1 0.22951 0.22951 4.399 0.283
## Residuals 1 0.05218 0.05218year.ps <- subset_samples(ps, sample_names(ps) %in% rownames(year.meta))dist.year <- vegdist(year.ps@otu_table, method = "bray")
year.a <- UniFrac(year.ps, weighted = TRUE)
year.b <- UniFrac(year.ps, weighted = FALSE)anova(betadisper(dist.year, year.ps@sam_data[["year"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.00817 0.0040847 0.2816 0.7554
## Residuals 74 1.07347 0.0145063anova(betadisper(dist.year, year.ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.06082 0.06082 6.3687 0.01373 *
## Residuals 75 0.71623 0.00955
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(dist.year, year.ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 1.7808 0.089039 2.177 0.01172 *
## Residuals 56 2.2904 0.040901
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(dist.year ~ year * sex * svl, data = year.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.year ~ year * sex * svl, data = year.meta)
## Df SumOfSqs R2 F Pr(>F)
## year 2 0.7575 0.03173 1.2430 0.256
## sex 1 0.3598 0.01507 1.1808 0.303
## svl 1 0.7252 0.03038 2.3800 0.039 *
## year:sex 2 0.5455 0.02285 0.8950 0.515
## year:svl 2 0.5587 0.02340 0.9168 0.503
## sex:svl 1 0.1308 0.00548 0.4291 0.812
## year:sex:svl 2 0.9899 0.04146 1.6242 0.142
## Residual 65 19.8072 0.82963
## Total 76 23.8747 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(dist.year ~ year + sex + svl, data = year.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.year ~ year + sex + svl, data = year.meta)
## Df SumOfSqs R2 F Pr(>F)
## year 2 0.7575 0.03173 1.2378 0.257
## sex 1 0.3598 0.01507 1.1759 0.304
## svl 1 0.7252 0.03038 2.3701 0.041 *
## Residual 72 22.0320 0.92282
## Total 76 23.8747 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(year.a, year.ps@sam_data[["year"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.04860 0.024300 2.3181 0.1056
## Residuals 74 0.77571 0.010483anova(betadisper(year.a, year.ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.01562 0.015625 1.4305 0.2355
## Residuals 75 0.81921 0.010923anova(betadisper(year.a, year.ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 0.26849 0.013425 1.2305 0.2656
## Residuals 56 0.61093 0.010910set.seed(1)
adonis2(year.a ~ year * sex * svl, data = year.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = year.a ~ year * sex * svl, data = year.meta)
## Df SumOfSqs R2 F Pr(>F)
## year 2 0.09181 0.03739 1.3951 0.251
## sex 1 0.02112 0.00860 0.6418 0.575
## svl 1 0.05443 0.02216 1.6541 0.172
## year:sex 2 0.04176 0.01701 0.6345 0.668
## year:svl 2 0.02456 0.01000 0.3731 0.886
## sex:svl 1 0.00741 0.00302 0.2252 0.889
## year:sex:svl 2 0.07573 0.03084 1.1507 0.328
## Residual 65 2.13876 0.87098
## Total 76 2.45556 1.00000set.seed(1)
adonis2(year.a ~ year + sex + svl, data = year.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = year.a ~ year + sex + svl, data = year.meta)
## Df SumOfSqs R2 F Pr(>F)
## year 2 0.09181 0.03739 1.4444 0.230
## sex 1 0.02112 0.00860 0.6645 0.563
## svl 1 0.05443 0.02216 1.7126 0.162
## Residual 72 2.28821 0.93185
## Total 76 2.45556 1.00000anova(betadisper(year.b, year.ps@sam_data[["year"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.01711 0.0085557 0.4682 0.628
## Residuals 74 1.35228 0.0182741anova(betadisper(year.b, year.ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.00366 0.0036645 0.2082 0.6495
## Residuals 75 1.31998 0.0175998anova(betadisper(year.b, year.ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 1.2073 0.060365 1.97 0.0243 *
## Residuals 56 1.7160 0.030643
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(year.b ~ year * sex * svl, data = year.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = year.b ~ year * sex * svl, data = year.meta)
## Df SumOfSqs R2 F Pr(>F)
## year 2 0.6271 0.03834 1.5379 0.089 .
## sex 1 0.2259 0.01382 1.1083 0.317
## svl 1 1.0476 0.06406 5.1390 0.001 ***
## year:sex 2 0.3172 0.01940 0.7779 0.703
## year:svl 2 0.3227 0.01973 0.7915 0.690
## sex:svl 1 0.0956 0.00584 0.4687 0.909
## year:sex:svl 2 0.4665 0.02853 1.1443 0.267
## Residual 65 13.2509 0.81028
## Total 76 16.3535 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(year.b ~ year + sex + svl, data = year.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = year.b ~ year + sex + svl, data = year.meta)
## Df SumOfSqs R2 F Pr(>F)
## year 2 0.6271 0.03834 1.5619 0.082 .
## sex 1 0.2259 0.01382 1.1256 0.305
## svl 1 1.0476 0.06406 5.2190 0.001 ***
## Residual 72 14.4529 0.88378
## Total 76 16.3535 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hibernation
hib.meta <- meta[meta$month == "September" | meta$month == "April",]
hib.ps <- subset_samples(ps, month == "September" | month == "April")aov.hib <- aov(Shannon ~ month * sex * svl, data = hib.meta)
plot(density(aov.hib$residuals))
qqnorm(aov.hib$residuals)
qqline(aov.hib$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.hib$residuals~aov.hib$fitted.values)
lines(lowess(aov.hib$fitted.values,aov.hib$residuals), col="blue")
aov.hib <- aov(Shannon ~ month * sex * svl + Error(toe.clip), data = hib.meta)
summary(aov.hib)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.073 0.0729 0.433 0.51414
## sex 1 0.450 0.4497 2.672 0.10995
## svl 1 1.529 1.5285 9.084 0.00446 **
## month:sex 1 0.003 0.0032 0.019 0.89088
## month:svl 1 0.329 0.3287 1.953 0.16994
## sex:svl 1 0.000 0.0002 0.001 0.97532
## month:sex:svl 1 0.045 0.0452 0.269 0.60712
## Residuals 40 6.731 0.1683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.29538 0.29538 2.039 0.289
## svl 1 0.02852 0.02852 0.197 0.701
## month:sex 1 0.06316 0.06316 0.436 0.577
## month:svl 1 0.05803 0.05803 0.401 0.591
## sex:svl 1 0.00298 0.00298 0.021 0.899
## month:sex:svl 1 0.15169 0.15169 1.047 0.414
## Residuals 2 0.28967 0.14484aov.hib.ob <- aov(log10(Observed) ~ month * sex * svl, data = hib.meta)
plot(density(aov.hib.ob$residuals))
qqnorm(aov.hib.ob$residuals)
qqline(aov.hib.ob$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.hib.ob$residuals~aov.hib.ob$fitted.values)
lines(lowess(aov.hib.ob$fitted.values,aov.hib.ob$residuals), col="blue")
aov.hib.ob <- aov(log10(Observed) ~ month * sex * svl + Error(toe.clip), data = hib.meta)
summary(aov.hib.ob)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.0150 0.01501 0.476 0.4944
## sex 1 0.0807 0.08069 2.557 0.1177
## svl 1 0.2845 0.28449 9.014 0.0046 **
## month:sex 1 0.0010 0.00101 0.032 0.8592
## month:svl 1 0.0651 0.06507 2.062 0.1588
## sex:svl 1 0.0000 0.00002 0.001 0.9796
## month:sex:svl 1 0.0090 0.00901 0.285 0.5961
## Residuals 40 1.2625 0.03156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.05344 0.05344 2.112 0.283
## svl 1 0.00618 0.00618 0.244 0.670
## month:sex 1 0.01222 0.01222 0.483 0.559
## month:svl 1 0.01113 0.01113 0.440 0.575
## sex:svl 1 0.00042 0.00042 0.017 0.909
## month:sex:svl 1 0.03129 0.03129 1.237 0.382
## Residuals 2 0.05061 0.02531aov.hib.pd <- aov(log10(PD) ~ month * sex * svl, data = hib.meta)
plot(density(aov.hib.pd$residuals))
qqnorm(aov.hib.pd$residuals)
qqline(aov.hib.pd$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.hib.pd$residuals~aov.hib.pd$fitted.values)
lines(lowess(aov.hib.pd$fitted.values,aov.hib.pd$residuals), col="blue")
aov.hib.pd <- aov(log10(PD) ~ month * sex * svl + Error(toe.clip), data = hib.meta)
summary(aov.hib.pd)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.0046 0.0046 0.121 0.72961
## sex 1 0.0426 0.0426 1.129 0.29439
## svl 1 0.3489 0.3489 9.244 0.00415 **
## month:sex 1 0.0140 0.0140 0.370 0.54659
## month:svl 1 0.1354 0.1354 3.587 0.06547 .
## sex:svl 1 0.0524 0.0524 1.388 0.24576
## month:sex:svl 1 0.0198 0.0198 0.526 0.47257
## Residuals 40 1.5097 0.0377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.07799 0.07799 2.867 0.232
## svl 1 0.01253 0.01253 0.461 0.567
## month:sex 1 0.03046 0.03046 1.120 0.401
## month:svl 1 0.03083 0.03083 1.133 0.399
## sex:svl 1 0.00539 0.00539 0.198 0.700
## month:sex:svl 1 0.09753 0.09753 3.586 0.199
## Residuals 2 0.05439 0.02720aov.hib <- aov(Shannon ~ month + sex + svl + Error(toe.clip), data = hib.meta)
aov.hib.ob <- aov(log10(Observed) ~ month + sex + svl + Error(toe.clip), data = hib.meta)
aov.hib.pd <- aov(log10(PD) ~ month + sex + svl + Error(toe.clip), data = hib.meta)
summary(aov.hib)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.073 0.0729 0.451 0.5052
## sex 1 0.450 0.4497 2.784 0.1023
## svl 1 1.529 1.5285 9.462 0.0036 **
## Residuals 44 7.108 0.1615
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.2954 0.29538 3.134 0.127
## svl 1 0.0285 0.02852 0.303 0.602
## Residuals 6 0.5655 0.09426summary(aov.hib.ob)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.0150 0.01501 0.494 0.48599
## sex 1 0.0807 0.08069 2.654 0.11040
## svl 1 0.2845 0.28449 9.358 0.00377 **
## Residuals 44 1.3376 0.03040
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.05344 0.05344 3.034 0.132
## svl 1 0.00618 0.00618 0.351 0.575
## Residuals 6 0.10568 0.01761summary(aov.hib.pd)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.0046 0.0046 0.116 0.73480
## sex 1 0.0426 0.0426 1.083 0.30375
## svl 1 0.3489 0.3489 8.867 0.00471 **
## Residuals 44 1.7312 0.0393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## month 1 0.07799 0.07799 2.140 0.194
## svl 1 0.01253 0.01253 0.344 0.579
## Residuals 6 0.21861 0.03643set.seed(1)
dist.hib <- vegdist(hib.ps@otu_table, method = "bray")
a.dist.hib <- UniFrac(hib.ps, weighted = TRUE)
b.dist.hib <- UniFrac(hib.ps, weighted = FALSE)anova(betadisper(dist.hib, hib.ps@sam_data[["month"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.00178 0.001784 0.0683 0.7948
## Residuals 54 1.41050 0.026120anova(betadisper(dist.hib, hib.ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.07803 0.078028 2.8678 0.09613 .
## Residuals 54 1.46924 0.027208
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(dist.hib, hib.ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 16 1.3448 0.084051 2.2825 0.01804 *
## Residuals 39 1.4361 0.036824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(dist.hib ~ month * sex * svl, data = hib.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.hib ~ month * sex * svl, data = hib.meta)
## Df SumOfSqs R2 F Pr(>F)
## month 1 0.1312 0.00863 0.4894 0.791
## sex 1 0.8486 0.05584 3.1644 0.012 *
## svl 1 0.4707 0.03097 1.7554 0.104
## month:sex 1 0.1382 0.00910 0.5155 0.776
## month:svl 1 0.0561 0.00369 0.2091 0.955
## sex:svl 1 0.3810 0.02507 1.4209 0.206
## month:sex:svl 1 0.2998 0.01973 1.1179 0.382
## Residual 48 12.8718 0.84697
## Total 55 15.1974 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(dist.hib ~ month + sex + svl, data = hib.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.hib ~ month + sex + svl, data = hib.meta)
## Df SumOfSqs R2 F Pr(>F)
## month 1 0.1312 0.00863 0.4964 0.762
## sex 1 0.8486 0.05584 3.2098 0.009 **
## svl 1 0.4707 0.03097 1.7806 0.124
## Residual 52 13.7469 0.90456
## Total 55 15.1974 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(a.dist.hib, hib.ps@sam_data[["month"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.001101 0.0011009 0.2411 0.6254
## Residuals 54 0.246577 0.0045662anova(betadisper(a.dist.hib, hib.ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.009066 0.0090664 1.8728 0.1768
## Residuals 54 0.261421 0.0048411anova(betadisper(a.dist.hib, hib.ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 16 0.1766 0.0110377 2.7453 0.005109 **
## Residuals 39 0.1568 0.0040206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(a.dist.hib ~ month * sex * svl, data = hib.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = a.dist.hib ~ month * sex * svl, data = hib.meta)
## Df SumOfSqs R2 F Pr(>F)
## month 1 0.01101 0.00781 0.5317 0.663
## sex 1 0.13501 0.09573 6.5215 0.002 **
## svl 1 0.10639 0.07544 5.1391 0.009 **
## month:sex 1 0.00974 0.00690 0.4704 0.719
## month:svl 1 0.01486 0.01054 0.7177 0.514
## sex:svl 1 0.09031 0.06404 4.3623 0.015 *
## month:sex:svl 1 0.04930 0.03496 2.3813 0.071 .
## Residual 48 0.99372 0.70460
## Total 55 1.41034 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(a.dist.hib ~ month + sex + svl, data = hib.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = a.dist.hib ~ month + sex + svl, data = hib.meta)
## Df SumOfSqs R2 F Pr(>F)
## month 1 0.01101 0.00781 0.4943 0.688
## sex 1 0.13501 0.09573 6.0630 0.006 **
## svl 1 0.10639 0.07544 4.7778 0.011 *
## Residual 52 1.15793 0.82103
## Total 55 1.41034 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(b.dist.hib, hib.ps@sam_data[["month"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.01813 0.018135 1.1809 0.282
## Residuals 54 0.82922 0.015356anova(betadisper(b.dist.hib, hib.ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.00100 0.0009992 0.0596 0.8081
## Residuals 54 0.90582 0.0167745anova(betadisper(b.dist.hib, hib.ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 16 0.91054 0.056908 3.1598 0.001698 **
## Residuals 39 0.70240 0.018010
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(b.dist.hib ~ month * sex * svl, data = hib.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = b.dist.hib ~ month * sex * svl, data = hib.meta)
## Df SumOfSqs R2 F Pr(>F)
## month 1 0.1889 0.01778 1.0252 0.367
## sex 1 0.2467 0.02322 1.3387 0.191
## svl 1 0.3857 0.03631 2.0928 0.023 *
## month:sex 1 0.2251 0.02119 1.2214 0.242
## month:svl 1 0.2279 0.02145 1.2365 0.254
## sex:svl 1 0.3611 0.03399 1.9593 0.031 *
## month:sex:svl 1 0.1419 0.01335 0.7697 0.663
## Residual 48 8.8459 0.83270
## Total 55 10.6232 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(b.dist.hib ~ month + sex + svl, data = hib.meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = b.dist.hib ~ month + sex + svl, data = hib.meta)
## Df SumOfSqs R2 F Pr(>F)
## month 1 0.1889 0.01778 1.0023 0.404
## sex 1 0.2467 0.02322 1.3088 0.202
## svl 1 0.3857 0.03631 2.0461 0.027 *
## Residual 52 9.8018 0.92269
## Total 55 10.6232 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Multivariate model
all.aov1 <- aov(Shannon ~ season * sex * svl, data = meta)
plot(density(all.aov1$residuals))
qqnorm(all.aov1$residuals)
qqline(all.aov1$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(all.aov1$residuals~all.aov1$fitted.values)
lines(lowess(all.aov1$fitted.values,all.aov1$residuals), col="blue")
summary(aov(Shannon ~ season * sex * svl + Error(toe.clip), data = meta))##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 1.27 0.635 3.650 0.0278 *
## sex 1 3.25 3.252 18.702 2.44e-05 ***
## svl 1 0.08 0.084 0.483 0.4880
## season:sex 2 0.00 0.000 0.000 0.9997
## season:svl 2 0.29 0.144 0.828 0.4386
## sex:svl 1 0.02 0.015 0.087 0.7689
## season:sex:svl 2 0.34 0.169 0.971 0.3807
## Residuals 194 33.73 0.174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.095 0.04759 0.190 0.827
## svl 1 0.120 0.11978 0.479 0.493
## season:sex 2 0.008 0.00390 0.016 0.985
## season:svl 2 0.146 0.07284 0.292 0.749
## sex:svl 1 0.041 0.04150 0.166 0.686
## season:sex:svl 2 0.621 0.31037 1.242 0.301
## Residuals 36 8.994 0.24983summary(aov(Shannon ~ season + sex + svl + Error(toe.clip), data = meta))##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 1.27 0.635 3.712 0.0261 *
## sex 1 3.25 3.252 19.016 2.07e-05 ***
## svl 1 0.08 0.084 0.491 0.4843
## Residuals 201 34.37 0.171
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.095 0.04759 0.209 0.813
## svl 1 0.120 0.11978 0.525 0.473
## Residuals 43 9.810 0.22813all.aov2 <- aov(Observed ~ season * sex * svl, data = meta)
plot(density(all.aov2$residuals))
qqnorm(all.aov2$residuals)
qqline(all.aov2$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(all.aov2$residuals~all.aov2$fitted.values)
lines(lowess(all.aov2$fitted.values,all.aov2$residuals), col="blue")
all.aov2 <- aov(log10(Observed + 1) ~ season * sex * svl, data = meta)
plot(density(all.aov2$residuals))
qqnorm(all.aov2$residuals)
qqline(all.aov2$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(all.aov2$residuals~all.aov2$fitted.values)
lines(lowess(all.aov2$fitted.values,all.aov2$residuals), col="blue")
summary(aov(log10(Observed + 1) ~ season * sex * svl + Error(toe.clip), data = meta))##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.237 0.1187 3.729 0.0258 *
## sex 1 0.591 0.5912 18.565 2.61e-05 ***
## svl 1 0.016 0.0163 0.511 0.4755
## season:sex 2 0.000 0.0000 0.000 0.9999
## season:svl 2 0.050 0.0248 0.779 0.4604
## sex:svl 1 0.003 0.0035 0.110 0.7409
## season:sex:svl 2 0.059 0.0294 0.923 0.3989
## Residuals 194 6.178 0.0318
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.0162 0.00812 0.177 0.839
## svl 1 0.0241 0.02408 0.524 0.474
## season:sex 2 0.0010 0.00050 0.011 0.989
## season:svl 2 0.0283 0.01413 0.307 0.737
## sex:svl 1 0.0075 0.00750 0.163 0.689
## season:sex:svl 2 0.1128 0.05641 1.227 0.305
## Residuals 36 1.6548 0.04597summary(aov(log10(Observed + 1) ~ season + sex + svl + Error(toe.clip), data = meta))##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.237 0.1187 3.794 0.0241 *
## sex 1 0.591 0.5912 18.892 2.2e-05 ***
## svl 1 0.016 0.0163 0.520 0.4716
## Residuals 201 6.290 0.0313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.0162 0.00812 0.194 0.825
## svl 1 0.0241 0.02408 0.574 0.453
## Residuals 43 1.8044 0.04196all.aov3 <- aov(log10(PD + 1) ~ season * sex * svl, data = meta)
plot(density(all.aov3$residuals))
qqnorm(all.aov3$residuals)
qqline(all.aov3$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(all.aov3$residuals~all.aov3$fitted.values)
lines(lowess(all.aov3$fitted.values,all.aov3$residuals), col="blue")
summary(aov(log10(PD + 1) ~ season * sex * svl + Error(toe.clip), data = meta))##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.113 0.05636 1.934 0.14742
## sex 1 0.297 0.29716 10.194 0.00164 **
## svl 1 0.044 0.04352 1.493 0.22325
## season:sex 2 0.023 0.01129 0.387 0.67953
## season:svl 2 0.002 0.00121 0.042 0.95930
## sex:svl 1 0.024 0.02383 0.818 0.36700
## season:sex:svl 2 0.024 0.01201 0.412 0.66296
## Residuals 194 5.655 0.02915
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.0007 0.00033 0.008 0.992
## svl 1 0.0361 0.03612 0.898 0.350
## season:sex 2 0.0427 0.02133 0.530 0.593
## season:svl 2 0.0349 0.01743 0.433 0.652
## sex:svl 1 0.0041 0.00412 0.102 0.751
## season:sex:svl 2 0.1691 0.08453 2.101 0.137
## Residuals 36 1.4481 0.04023summary(aov(log10(PD + 1) ~ season + sex +svl + Error(toe.clip), data = meta))##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.113 0.05636 1.978 0.14106
## sex 1 0.297 0.29716 10.427 0.00145 **
## svl 1 0.044 0.04352 1.527 0.21799
## Residuals 201 5.728 0.02850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## season 2 0.0007 0.00033 0.008 0.992
## svl 1 0.0361 0.03612 0.914 0.344
## Residuals 43 1.6989 0.03951set.seed(1)
dist <- vegdist(ps@otu_table, method = "bray")
a.dist <- UniFrac(ps, weighted = TRUE)
b.dist <- UniFrac(ps, weighted = FALSE)anova(betadisper(dist, ps@sam_data[["season"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.0821 0.041044 2.4919 0.08481 .
## Residuals 249 4.1012 0.016471
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(dist, ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.3335 0.33349 19.957 1.201e-05 ***
## Residuals 250 4.1776 0.01671
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(dist, ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 1.7781 0.088905 4.075 7.252e-08 ***
## Residuals 231 5.0397 0.021817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(dist ~ season * sex * svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist ~ season * sex * svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 1.546 0.02127 2.7498 0.007 **
## sex 1 0.795 0.01093 2.8265 0.032 *
## svl 1 0.890 0.01225 3.1661 0.014 *
## season:sex 2 0.412 0.00567 0.7323 0.666
## season:svl 2 0.929 0.01278 1.6525 0.101
## sex:svl 1 0.317 0.00436 1.1270 0.363
## season:sex:svl 2 0.321 0.00442 0.5710 0.815
## Residual 240 67.484 0.92832
## Total 251 72.695 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1removing interaction term
set.seed(1)
adonis2(dist ~ season + sex + svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist ~ season + sex + svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 1.546 0.02127 2.7494 0.007 **
## sex 1 0.795 0.01093 2.8260 0.032 *
## svl 1 0.890 0.01225 3.1657 0.014 *
## Residual 247 69.463 0.95555
## Total 251 72.695 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(a.dist, ps@sam_data[["season"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.02152 0.0107583 1.626 0.1988
## Residuals 249 1.64750 0.0066165anova(betadisper(a.dist, ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.0379 0.037900 6.0239 0.0148 *
## Residuals 250 1.5729 0.006292
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(a.dist, ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 0.33118 0.016559 2.4089 0.0009663 ***
## Residuals 231 1.58789 0.006874
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(a.dist ~ season * sex * svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = a.dist ~ season * sex * svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.1204 0.01643 2.1057 0.071 .
## sex 1 0.0890 0.01213 3.1108 0.034 *
## svl 1 0.0835 0.01139 2.9214 0.045 *
## season:sex 2 0.0586 0.00800 1.0252 0.370
## season:svl 2 0.0572 0.00780 0.9993 0.369
## sex:svl 1 0.0211 0.00288 0.7387 0.496
## season:sex:svl 2 0.0387 0.00528 0.6766 0.664
## Residual 240 6.8628 0.93609
## Total 251 7.3314 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(a.dist ~ season + sex + svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = a.dist ~ season + sex + svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.1204 0.01643 2.1130 0.071 .
## sex 1 0.0890 0.01213 3.1216 0.034 *
## svl 1 0.0835 0.01139 2.9316 0.044 *
## Residual 247 7.0384 0.96005
## Total 251 7.3314 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(b.dist, ps@sam_data[["season"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.0323 0.016158 1.0246 0.3604
## Residuals 249 3.9266 0.015769anova(betadisper(b.dist, ps@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.0056 0.0056164 0.3573 0.5505
## Residuals 250 3.9295 0.0157182anova(betadisper(b.dist, ps@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 0.9502 0.047512 2.8651 7.779e-05 ***
## Residuals 231 3.8308 0.016583
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(b.dist ~ season * sex * svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = b.dist ~ season * sex * svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.990 0.01859 2.4028 0.002 **
## sex 1 0.601 0.01129 2.9173 0.006 **
## svl 1 1.042 0.01955 5.0539 0.001 ***
## season:sex 2 0.283 0.00531 0.6866 0.848
## season:svl 2 0.401 0.00753 0.9737 0.460
## sex:svl 1 0.206 0.00387 0.9996 0.393
## season:sex:svl 2 0.280 0.00526 0.6797 0.864
## Residual 240 49.465 0.92859
## Total 251 53.269 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(b.dist ~ season + sex + svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = b.dist ~ season + sex + svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.990 0.01859 2.4157 0.002 **
## sex 1 0.601 0.01129 2.9330 0.006 **
## svl 1 1.042 0.01955 5.0811 0.001 ***
## Residual 247 50.635 0.95056
## Total 251 53.269 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Common Core
ps.core <- core(ps.real, detection = 0, prevalence = .5)
ps.core@tax_table## Taxonomy Table: [9 taxa by 6 taxonomic ranks]:
## Kingdom Phylum Class
## ASV_19 "Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_6 "Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_22 "Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_18 "Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_23 "Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_4 "Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_10 "Bacteria" "Proteobacteria" "Gammaproteobacteria"
## ASV_7 "Bacteria" "Epsilonbacteraeota" "Campylobacteria"
## ASV_12 "Bacteria" "Epsilonbacteraeota" "Campylobacteria"
## Order Family Genus
## ASV_19 "Enterobacteriales" "Enterobacteriaceae" "Izhakiella"
## ASV_6 "Enterobacteriales" "Enterobacteriaceae" NA
## ASV_22 "Enterobacteriales" "Enterobacteriaceae" "Izhakiella"
## ASV_18 "Enterobacteriales" "Enterobacteriaceae" "Izhakiella"
## ASV_23 "Enterobacteriales" "Enterobacteriaceae" "Izhakiella"
## ASV_4 "Enterobacteriales" "Enterobacteriaceae" NA
## ASV_10 "Enterobacteriales" "Enterobacteriaceae" "Izhakiella"
## ASV_7 "Campylobacterales" "Helicobacteraceae" "Helicobacter"
## ASV_12 "Campylobacterales" "Helicobacteraceae" "Helicobacter"f.ps <- subset_samples(ps.real, sex == "F")
f.ps.core <- microbiome::core(f.ps, detection = 0, prevalence = .4)m.ps <- subset_samples(ps.real, sex == "M")
m.ps.core <- microbiome::core(m.ps, detection = 0, prevalence = .4)f.ps.fam <- subset_samples(fam.ps, sex == "F")
f.ps.core.fam <- microbiome::core(f.ps.fam , detection = 0, prevalence = .5)m.ps.fam <- subset_samples(fam.ps, sex == "M")
m.ps.core.fam <- microbiome::core(m.ps.fam, detection = 0, prevalence = .5)pre.core <- subset_samples(ps, season == "pre")
pre.core <- core(pre.core, detection = 0, prevalence = 0.5)
breeding.core <- subset_samples(ps, season == "breeding")
breeding.core <- core(breeding.core, detection = 0, prevalence = 0.5)
post.core <- subset_samples(ps, season == "post")
post.core <- core(post.core, detection = 0, prevalence = 0.5)pre.core.fam <- subset_samples(fam.ps, season == "pre")
pre.core.fam <- core(pre.core.fam, detection = 0, prevalence = 0.5)
breeding.core.fam <- subset_samples(fam.ps, season == "breeding")
breeding.core.fam <- core(breeding.core.fam, detection = 0, prevalence = 0.5)
post.core.fam <- subset_samples(fam.ps, season == "post")
post.core.fam <- core(post.core.fam, detection = 0, prevalence = 0.5)Temporal Core
year.venn <- get_vennlist(year.ps, "year")
shared.year <- Reduce(intersect, year.venn)
sex.venn <- get_vennlist(ps, "sex")
shared.sex <- Reduce(intersect, sex.venn)
season.venn <- get_vennlist(ps, "season")
shared.season <- Reduce(intersect, season.venn)
core.taxa <- colnames(counts)[colnames(counts) %in% shared.year & colnames(counts) %in% shared.sex & colnames(counts) %in% shared.season]ps.core <- prune_taxa(core.taxa, ps)n.core <- vegdist(ps.core@otu_table, method = "bray")
a.core <- UniFrac(ps.core, weighted = TRUE)
b.core <- UniFrac(ps.core, weighted = FALSE)anova(betadisper(n.core, ps.core@sam_data[["season"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.1651 0.082572 4.1715 0.01652 *
## Residuals 249 4.9288 0.019794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(n.core, ps.core@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.4398 0.43979 20.86 7.763e-06 ***
## Residuals 250 5.2708 0.02108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(n.core, ps.core@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 1.8631 0.093153 3.3812 4.077e-06 ***
## Residuals 231 6.3640 0.027550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(n.core ~ season * sex * svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = n.core ~ season * sex * svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 1.483 0.02204 2.8517 0.012 *
## sex 1 0.809 0.01202 3.1114 0.033 *
## svl 1 0.890 0.01322 3.4223 0.017 *
## season:sex 2 0.304 0.00452 0.5852 0.730
## season:svl 2 0.783 0.01164 1.5060 0.164
## sex:svl 1 0.319 0.00474 1.2255 0.334
## season:sex:svl 2 0.303 0.00451 0.5833 0.761
## Residual 240 62.425 0.92732
## Total 251 67.318 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(n.core ~ season + sex + svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = n.core ~ season + sex + svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 1.483 0.02204 2.8566 0.011 *
## sex 1 0.809 0.01202 3.1168 0.032 *
## svl 1 0.890 0.01322 3.4282 0.015 *
## Residual 247 64.135 0.95272
## Total 251 67.318 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(a.core, ps.core@sam_data[["season"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.0443 0.022146 0.8689 0.4207
## Residuals 249 6.3462 0.025487anova(betadisper(a.core, ps.core@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.3168 0.31676 14.621 0.000166 ***
## Residuals 250 5.4162 0.02166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(a.core, ps.core@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 1.2626 0.063129 2.2897 0.001827 **
## Residuals 231 6.3689 0.027571
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(a.core ~ season * sex * svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = a.core ~ season * sex * svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.2704 0.01086 1.3920 0.267
## sex 1 0.3100 0.01245 3.1908 0.068 .
## svl 1 0.4630 0.01859 4.7657 0.033 *
## season:sex 2 0.1049 0.00421 0.5402 0.602
## season:svl 2 0.1932 0.00776 0.9943 0.364
## sex:svl 1 0.0404 0.00162 0.4161 0.598
## season:sex:svl 2 0.2046 0.00822 1.0530 0.380
## Residual 240 23.3141 0.93629
## Total 251 24.9006 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(a.core ~ season + sex + svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = a.core ~ season + sex + svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.2704 0.01086 1.3999 0.264
## sex 1 0.3100 0.01245 3.2091 0.066 .
## svl 1 0.4630 0.01859 4.7931 0.030 *
## Residual 247 23.8573 0.95810
## Total 251 24.9006 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(b.core, ps.core@sam_data[["season"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.0671 0.033533 1.5253 0.2196
## Residuals 249 5.4742 0.021985anova(betadisper(b.core, ps.core@sam_data[["sex"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 1 0.1413 0.141283 6.059 0.01451 *
## Residuals 250 5.8295 0.023318
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(b.core, ps.core@sam_data[["svl"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 20 1.1031 0.055154 2.2534 0.002213 **
## Residuals 231 5.6540 0.024476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(b.core ~ season * sex * svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = b.core ~ season * sex * svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.560 0.01417 1.8139 0.107
## sex 1 0.413 0.01044 2.6738 0.044 *
## svl 1 0.612 0.01547 3.9604 0.007 **
## season:sex 2 0.194 0.00490 0.6270 0.719
## season:svl 2 0.266 0.00673 0.8622 0.520
## sex:svl 1 0.083 0.00209 0.5342 0.711
## season:sex:svl 2 0.355 0.00897 1.1485 0.336
## Residual 240 37.067 0.93724
## Total 251 39.550 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(1)
adonis2(b.core ~ season + sex + svl, data = meta)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = b.core ~ season + sex + svl, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## season 2 0.560 0.01417 1.8227 0.104
## sex 1 0.413 0.01044 2.6867 0.043 *
## svl 1 0.612 0.01547 3.9796 0.006 **
## Residual 247 37.965 0.95993
## Total 251 39.550 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Female reproductive State
meta.f <- meta.full[meta.full$sex == "F",]
ps.f <- subset_samples(ps.full, sex == "F")aov.rs.shan <- aov(Shannon ~ RS, data = meta.f)
plot(density(aov.rs.shan$residuals))
qqnorm(aov.rs.shan$residuals)
qqline(aov.rs.shan$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.rs.shan$residuals~aov.rs.shan$fitted.values)
lines(lowess(aov.rs.shan$fitted.values,aov.rs.shan$residuals), col="blue")
aov.rs.ob <- aov(Observed ~ RS, data = meta.f)
plot(density(aov.rs.ob$residuals))
qqnorm(aov.rs.ob$residuals)
qqline(aov.rs.ob$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.rs.ob$residuals~aov.rs.ob$fitted.values)
lines(lowess(aov.rs.ob$fitted.values,aov.rs.ob$residuals), col="blue")
aov.rs.ob <- aov(log10(Observed) ~ RS, data = meta.f)
plot(density(aov.rs.ob$residuals))
qqnorm(aov.rs.ob$residuals)
qqline(aov.rs.ob$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.rs.ob$residuals~aov.rs.ob$fitted.values)
lines(lowess(aov.rs.ob$fitted.values,aov.rs.ob$residuals), col="blue")
aov.rs.pd <- aov(PD ~ RS, data = meta.f)
plot(density(aov.rs.pd$residuals))
qqnorm(aov.rs.pd$residuals)
qqline(aov.rs.pd$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.rs.pd$residuals~aov.rs.pd$fitted.values)
lines(lowess(aov.rs.pd$fitted.values,aov.rs.pd$residuals), col="blue")
aov.rs.pd <- aov(log10(PD) ~ RS, data = meta.f)
plot(density(aov.rs.pd$residuals))
qqnorm(aov.rs.pd$residuals)
qqline(aov.rs.pd$residuals, datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75))
plot(aov.rs.pd$residuals~aov.rs.pd$fitted.values)
lines(lowess(aov.rs.pd$fitted.values,aov.rs.pd$residuals), col="blue")
aov.rs.shan <- aov(Shannon ~ RS + Error(toe.clip), data = meta.f)
summary(aov.rs.shan)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## RS 3 1.093 0.3642 1.869 0.139
## Residuals 103 20.068 0.1948
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## RS 3 0.102 0.03405 0.146 0.932
## Residuals 29 6.778 0.23373aov.rs.ob <- aov(log10(Observed) ~ RS + Error(toe.clip), data = meta.f)
summary(aov.rs.ob)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## RS 3 0.211 0.07045 1.892 0.136
## Residuals 103 3.836 0.03724
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## RS 3 0.0191 0.00635 0.142 0.934
## Residuals 29 1.2966 0.04471aov.rs.pd <- aov(log10(PD) ~ RS + Error(toe.clip), data = meta.f)
summary(aov.rs.pd)##
## Error: toe.clip
## Df Sum Sq Mean Sq F value Pr(>F)
## RS 3 0.363 0.12111 1.612 0.191
## Residuals 103 7.738 0.07513
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## RS 3 0.0709 0.02362 0.273 0.844
## Residuals 29 2.5110 0.08659set.seed(1)
dist.f <- vegdist(ps.f@otu_table, method = "bray")
a.dist.f <- UniFrac(ps.f, weighted = TRUE)
b.dist.f <- UniFrac(ps.f, weighted = FALSE)anova(betadisper(dist.f, ps.f@sam_data[["RS"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.05634 0.018782 0.9517 0.4177
## Residuals 135 2.66427 0.019735set.seed(1)
adonis2(dist.f ~ RS, data = meta.f)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.f ~ RS, data = meta.f)
## Df SumOfSqs R2 F Pr(>F)
## RS 3 1.407 0.03952 1.8518 0.032 *
## Residual 135 34.189 0.96048
## Total 138 35.595 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(betadisper(a.dist.f, ps.f@sam_data[["RS"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.01966 0.0065539 1.3326 0.2664
## Residuals 135 0.66392 0.0049180set.seed(1)
adonis2(a.dist.f ~ RS, data = meta.f)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = a.dist.f ~ RS, data = meta.f)
## Df SumOfSqs R2 F Pr(>F)
## RS 3 0.1009 0.02953 1.3691 0.201
## Residual 135 3.3153 0.97047
## Total 138 3.4162 1.00000anova(betadisper(b.dist.f, ps.f@sam_data[["RS"]]))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 3 0.07637 0.025457 1.3964 0.2467
## Residuals 135 2.46108 0.018230set.seed(1)
adonis2(b.dist.f ~ RS, data = meta.f)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = b.dist.f ~ RS, data = meta.f)
## Df SumOfSqs R2 F Pr(>F)
## RS 3 0.9114 0.03069 1.4247 0.07 .
## Residual 135 28.7874 0.96931
## Total 138 29.6987 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1