Chapter 7 Is the GM of cold- vs warm-adapted populations similar in the wild?
load("data/data_03122025.Rdata")
load("data/objects_03122025.Rdata")
load("data/alpha_16092025.Rdata")
load("data/beta_16092025.Rdata")7.1 Shared and unique MAGs
locationcolors=c('#c4d7d1','#e08683')
subset_meta <- sample_metadata %>%
filter(time_point=="Wild")
genome_counts_rel_fil<- genome_counts_filt %>%
select(one_of(c("genome", subset_meta$Tube_code))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where( ~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames(., "genome")
genome_counts_rel_pa=1*(genome_counts_rel_fil>0)
table_upset_analysis_cont=t(aggregate(t(genome_counts_rel_pa),by=list(subset_meta$Population),FUN=sum)[,-1])
colnames(table_upset_analysis_cont)=levels(as.factor(subset_meta$Population))
table_upset_analysis=(table_upset_analysis_cont>0)*1
table_upset_analysis=data.frame(table_upset_analysis)
table_upset_analysis=apply(table_upset_analysis,2,as.integer)
rownames(table_upset_analysis) <- rownames(genome_counts_rel_pa)
upset(as.data.frame(table_upset_analysis),
keep.order = T,
sets = rev(c("Cold","Warm")),
sets.bar.color= rev(locationcolors),
mb.ratio = c(0.55, 0.45), order.by = "freq")
7.2 Shapiro test
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Wild") %>%
filter(metric=="richness") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.272924alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Wild") %>%
filter(metric=="neutral") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.6043408alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Wild") %>%
filter(metric=="phylogenetic") %>%
summarize(shapiro_p_value = shapiro.test(value)$p.value) %>%
pull(shapiro_p_value)[1] 0.88765387.3 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point=="Wild") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(alpha=0.5) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080', "#d57d2c")) +
scale_fill_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#00808050', "#d57d2c50")) +
facet_wrap(. ~ metric,scales = "free") +
coord_cartesian(xlim = c(1, NA)) +
stat_compare_means(method="t.test",size=3, label.x=.58) +
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_text(size=10),
axis.text.x = element_blank(),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")
7.4 Beta diversity
Number of samples used
samples_to_keep <- sample_metadata %>%
filter(time_point == "Wild") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point == "Wild")
length(samples_to_keep)[1] 27Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$species) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.0000 0.0000002 0 999 0.997
Residuals 25 0.2559 0.0102362
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.997
Podarcis_muralis 0.99627 adonis2(richness ~ species,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.559800 | 0.2109623 | 6.684166 | 0.001 |
| Residual | 25 | 5.833935 | 0.7890377 | NA | NA |
| Total | 26 | 7.393735 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$species) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000033 0.000033 0.003 999 0.955
Residuals 25 0.270374 0.010815
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.955
Podarcis_muralis 0.95642 adonis2(neutral ~ species,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.926602 | 0.2614502 | 8.850119 | 0.001 |
| Residual | 25 | 5.442304 | 0.7385498 | NA | NA |
| Total | 26 | 7.368906 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$species) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03686 0.036864 2.539 999 0.124
Residuals 25 0.36298 0.014519
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis 0.131
Podarcis_muralis 0.12363 adonis2(phylo ~ species,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3229671 | 0.2158407 | 6.881277 | 0.001 |
| Residual | 25 | 1.1733546 | 0.7841593 | NA | NA |
| Total | 26 | 1.4963217 | 1.0000000 | NA | NA |
NMDS
richness %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta, by = join_by(sample == Tube_code)) %>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species, shape=type)) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080',"#d57d2c")) +
scale_shape_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold-intervention"),
labels=c("Cold-control", "Warm-control", "Cold-intervention"),
values=c("circle","square","triangle")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
theme_classic() +
theme(legend.position="none")
neutral %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta, by = join_by(sample == Tube_code))%>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species, shape=type)) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080',"#d57d2c")) +
scale_shape_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold-intervention"),
labels=c("Cold-control", "Warm-control", "Cold-intervention"),
values=c("circle","square","triangle")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
theme_classic() +
theme(legend.position="none")
phylo %>%
metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
left_join(subset_meta, by = join_by(sample == Tube_code)) %>%
group_by(species) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(., aes(x=NMDS1,y=NMDS2, color=species, shape=type)) +
scale_color_manual(name="species",
breaks=c("Podarcis_muralis","Podarcis_liolepis"),
labels=c("Podarcis muralis","Podarcis liolepis"),
values=c('#008080',"#d57d2c")) +
scale_shape_manual(name="Type",
breaks=c("Cold_control", "Warm_control", "Cold-intervention"),
labels=c("Cold-control", "Warm-control", "Cold-intervention"),
values=c("circle","square","triangle")) +
geom_point(size=2) +
geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
labs(y = element_blank (), x = "NMDS1 \n Phylogenetic beta diversity") +
theme_classic() +
theme(legend.position = "none") 
7.5 Functional capacity
7.5.1 Metabolic capacity index at functional level
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Wild") %>%
group_by(species) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
species MCI sd
<chr> <dbl> <dbl>
1 Podarcis_liolepis 0.327 0.0245
2 Podarcis_muralis 0.346 0.0195MCI_element_wild <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
filter(time_point == "Wild")
shapiro.test(MCI_element_wild$value) #normality test
Shapiro-Wilk normality test
data: MCI_element_wild$value
W = 0.96381, p-value = 0.4494
F test to compare two variances
data: value by species
F = 1.5711, num df = 8, denom df = 17, p-value = 0.412
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.5132788 6.3689920
sample estimates:
ratio of variances
1.571133
Two Sample t-test
data: value by species
t = -2.2537, df = 25, p-value = 0.03323
alternative hypothesis: true difference in means between group Podarcis_liolepis and group Podarcis_muralis is not equal to 0
95 percent confidence interval:
-0.037427955 -0.001685206
sample estimates:
mean in group Podarcis_liolepis mean in group Podarcis_muralis
0.3265678 0.3461243 7.5.2 Differences in functional pathways
element_gift_wild <- element_gift %>%
filter(time_point == "Wild") %>%
select(-time_point, -type) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(Tube_code, where( ~ is.numeric(.) && sum(.) > 0)) %>%
left_join(sample_metadata[c(1, 4)], by = "Tube_code")significant_elements_wild <- element_gift_wild %>%
pivot_longer(-c(Tube_code, Population),
names_to = "trait",
values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ Population)$p.value) %>%
mutate(p_adjust = p.adjust(p_value, method = "BH")) %>%
filter(p_adjust < 0.05) %>%
rownames_to_column(., "Elements") %>%
select(-Elements)
element_gift_sig_wild <- element_gift_wild %>%
select(Tube_code, all_of(intersect(
significant_elements_wild$trait,
colnames(element_gift_wild)
))) %>%
left_join(., sample_metadata[c(1, 4)], by = join_by(Tube_code == Tube_code))
difference_table_wild <- element_gift_sig_wild %>%
select(-Tube_code) %>%
group_by(Population) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(., uniqueGIFT_db[c(1, 3, 4)], by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference = Cold - Warm) %>%
mutate(group_color = ifelse(Difference < 0, "Warm", "Cold")) difference_table_wild %>%
ggplot(aes(x = forcats::fct_reorder(Function, Difference),y = Difference,fill = group_color)) +
geom_col() +
scale_fill_manual(values = c('#008080', "#d57d2c")) +
geom_hline(yintercept = 0) +
coord_flip() +
theme(
axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
legend.title = element_blank(),
panel.background = element_blank(),
panel.grid.major = element_line(
size = 0.15,
linetype = 'solid',
colour = "grey"
)
) +
xlab("Function") +
ylab("Mean difference")