Alpha diversities
#Calculate Hill numbers
richness <- genome_counts %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=0) %>%
t() %>%
as.data.frame() %>%
rename(richness=1) %>%
rownames_to_column(var="sample")
neutral <- genome_counts %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1) %>%
t() %>%
as.data.frame() %>%
rename(neutral=1) %>%
rownames_to_column(var="sample")
phylogenetic <- genome_counts %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1,tree=genome_tree) %>%
t() %>%
as.data.frame() %>%
rename(phylogenetic=1) %>%
rownames_to_column(var="sample")
# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method="gower")
functional <- genome_counts %>%
column_to_rownames(var="genome") %>%
select(where(~!all(. == 0))) %>%
hilldiv(.,q=1,dist=dist) %>%
t() %>%
as.data.frame() %>%
rename(functional=1) %>%
rownames_to_column(var="sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
#Merge into a single table
#alpha_div <- cbind(sample=colnames(genome_counts[-1]),richness=q0n,neutral=round(q1n,3),phylo=round(q1p,3),func=round(q1f,3)) %>%
alpha_div <- richness %>%
full_join(neutral,by=join_by(sample==sample)) %>%
full_join(phylogenetic,by=join_by(sample==sample)) %>%
full_join(functional,by=join_by(sample==sample)) %>%
pivot_longer(-sample, names_to = "metric", values_to = "diversity") %>%
left_join(., sample_metadata, by = join_by(sample == sample)) %>%
filter(cage %in% c("C03","C04","C05","C06","C10","C11","C12","C13")) %>%
mutate(diversity = as.numeric(diversity)) %>%
mutate(metric = factor(metric, levels = c("richness","neutral","phylogenetic","functional"))) #sort metrics
alpha_div %>%
filter(metric == "richness") %>%
lmerTest::lmer(diversity ~ dominance + group + (1 | animal) + (1 | cage), data = ., REML = FALSE) %>%
summary()
boundary (singular) fit: see help('isSingular')
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: diversity ~ dominance + group + (1 | animal) + (1 | cage)
Data: .
AIC BIC logLik deviance df.resid
2128.2 2148.4 -1058.1 2116.2 210
Scaled residuals:
Min 1Q Median 3Q Max
-2.94196 -0.42123 0.09392 0.54104 2.55048
Random effects:
Groups Name Variance Std.Dev.
animal (Intercept) 71.62 8.463
cage (Intercept) 0.00 0.000
Residual 990.74 31.476
Number of obs: 216, groups: animal, 39; cage, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 84.779 6.559 41.660 12.926 3.6e-16 ***
dominance 17.870 8.793 48.704 2.032 0.0476 *
groupvariable -15.273 5.702 38.397 -2.679 0.0108 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) domnnc
dominance -0.671
groupvaribl -0.634 0.001
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
alpha_div %>%
filter(metric == "neutral") %>%
lmerTest::lmer(diversity ~ dominance + group + (1 | animal) + (1 | cage), data = ., REML = FALSE) %>%
summary()
boundary (singular) fit: see help('isSingular')
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: diversity ~ dominance + group + (1 | animal) + (1 | cage)
Data: .
AIC BIC logLik deviance df.resid
1693.9 1714.1 -840.9 1681.9 210
Scaled residuals:
Min 1Q Median 3Q Max
-2.33081 -0.59501 0.03916 0.64883 2.56233
Random effects:
Groups Name Variance Std.Dev.
animal (Intercept) 0 0.00
cage (Intercept) 0 0.00
Residual 141 11.87
Number of obs: 216, groups: animal, 39; cage, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 28.686 2.088 216.000 13.739 <2e-16 ***
dominance 4.100 2.834 216.000 1.446 0.1495
groupvariable -4.028 1.804 216.000 -2.233 0.0266 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) domnnc
dominance -0.679
groupvaribl -0.626 0.003
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
alpha_div %>%
filter(metric == "phylogenetic") %>%
lmerTest::lmer(dominance ~ diversity + group + (1 | animal) + (1 | cage), data = ., REML = FALSE) %>%
summary()
boundary (singular) fit: see help('isSingular')
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: dominance ~ diversity + group + (1 | animal) + (1 | cage)
Data: .
AIC BIC logLik deviance df.resid
-310.9 -290.6 161.4 -322.9 210
Scaled residuals:
Min 1Q Median 3Q Max
-4.8063 -0.3472 0.0156 0.2933 2.7098
Random effects:
Groups Name Variance Std.Dev.
animal (Intercept) 0.074930 0.27373
cage (Intercept) 0.000000 0.00000
Residual 0.006122 0.07824
Number of obs: 216, groups: animal, 39; cage, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.745e-01 9.055e-02 4.525e+01 5.241 4.06e-06 ***
diversity 5.659e-03 5.422e-03 1.785e+02 1.044 0.298
groupvariable 2.169e-03 1.011e-01 3.894e+01 0.021 0.983
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) dvrsty
diversity -0.271
groupvaribl -0.829 -0.005
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
alpha_div %>%
filter(metric == "functional") %>%
lmerTest::lmer(dominance ~ diversity + group + (1 | animal) + (1 | cage), data = ., REML = FALSE) %>%
summary()
boundary (singular) fit: see help('isSingular')
Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: dominance ~ diversity + group + (1 | animal) + (1 | cage)
Data: .
AIC BIC logLik deviance df.resid
-311.8 -291.5 161.9 -323.8 210
Scaled residuals:
Min 1Q Median 3Q Max
-4.7854 -0.3173 0.0196 0.3444 2.6510
Random effects:
Groups Name Variance Std.Dev.
animal (Intercept) 0.07500 0.27387
cage (Intercept) 0.00000 0.00000
Residual 0.00609 0.07804
Number of obs: 216, groups: animal, 39; cage, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.849e-01 1.195e-01 1.171e+02 3.22 0.00166 **
diversity 8.293e-02 5.880e-02 1.779e+02 1.41 0.16018
groupvariable 3.035e-03 1.011e-01 3.894e+01 0.03 0.97621
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) dvrsty
diversity -0.684
groupvaribl -0.631 0.003
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
alpha_div %>%
filter(metric == "richness") %>%
filter(treatment != "AN") %>%
ggplot(aes(x=dominance, y=diversity, group=cage)) +
geom_point()+
geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
facet_nested(. ~ group + cage) +
theme_minimal() +
theme(axis.text.x = element_blank()) +
labs(y="Richness")
alpha_div %>%
filter(metric == "neutral") %>%
filter(treatment != "AN") %>%
ggplot(aes(x=dominance, y=diversity, group=cage)) +
geom_point()+
geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
facet_nested(. ~ group + cage) +
theme_minimal() +
theme(axis.text.x = element_blank()) +
labs(y="Neutral diversity")
alpha_div %>%
filter(metric == "phylogenetic") %>%
filter(treatment != "AN") %>%
ggplot(aes(x=dominance, y=diversity, group=cage)) +
geom_point()+
geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
facet_nested(. ~ group + cage) +
theme_minimal() +
theme(axis.text.x = element_blank()) +
labs(y="Phylogenetic diversity")
alpha_div %>%
filter(metric == "functional") %>%
filter(treatment != "AN") %>%
ggplot(aes(x=dominance, y=diversity, group=cage)) +
geom_point()+
geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
facet_nested(. ~ group + cage) +
theme_minimal() +
theme(axis.text.x = element_blank()) +
labs(y="Functional diversity")
#Plot diversities
alpha_plot_neutral <- alpha_div %>%
filter(metric=="neutral") %>%
mutate(sample = factor(sample, levels = sample_metadata %>%
arrange(animal,match(treatment,rev(c("OP","HT","HR","CT","CR","DT","DR","AN","T1","T2","T3")))) %>%
select(sample) %>%
filter(!is.na(sample)) %>%
pull() )) %>% # sort samples
filter(!is.na(animal)) %>%
ggplot(aes(x=diversity, y=sample)) +
geom_bar(stat='identity', fill="#6c9ebc") +
facet_nested(group + cage + animal ~ metric, scales="free") + #facet per treatment and animal
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line( size=.1, color="grey" ),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text.y = element_blank()
)
alpha_plot_phylo <- alpha_div %>%
filter(metric=="phylogenetic") %>%
mutate(sample = factor(sample, levels = sample_metadata %>%
arrange(animal,match(treatment,rev(c("OP","HT","HR","CT","CR","DT","DR","AN","T1","T2","T3")))) %>%
select(sample) %>%
filter(!is.na(sample)) %>%
pull() )) %>% # sort samples
filter(!is.na(animal)) %>%
ggplot(aes(x=diversity, y=sample)) +
geom_bar(stat='identity', fill="#6c9ebc") +
facet_nested(group + cage + animal ~ metric, scales="free") + #facet per treatment and animal
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line( size=.1, color="grey" ),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text.y = element_blank()
)
alpha_plot_func <- alpha_div %>%
filter(metric=="functional") %>%
mutate(sample = factor(sample, levels = sample_metadata %>%
arrange(animal,match(treatment,rev(c("OP","HT","HR","CT","CR","DT","DR","AN","T1","T2","T3")))) %>%
select(sample) %>%
filter(!is.na(sample)) %>%
pull() )) %>% # sort samples
filter(!is.na(animal)) %>%
ggplot(aes(x=diversity, y=sample)) +
geom_bar(stat='identity', fill="#6c9ebc") +
facet_nested(group + cage + animal ~ metric, scales="free") + #facet per treatment and animal
coord_cartesian(xlim = c(1, NA)) +
theme_classic() +
theme(strip.text.y = element_text(angle = 0),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)
)
grid.arrange(alpha_plot_neutral, alpha_plot_phylo, alpha_plot_func, nrow = 1)
#Calculate Hill numbers
beta_neutral <- genome_counts %>%
column_to_rownames(var="genome") %>%
hillpair(.,q=1, metric="C")