---
title: 'Rmarkdown report for: Diurnal oscillations in gut bacterial load and composition eclipse seasonal and lifetime dynamics in wild meerkats, Suricata suricatta'
author: "Alice Risely"
date: "2nd September 2021"
output:
  html_document:
    toc: yes
    toc_depth: '3'
    theme: journal
    df_print: paged
  pdf_document:
    toc: yes
    toc_depth: 3
---

```{r echo = F}
knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      # record the current time before each chunk
      now <<- Sys.time()
    } else {
      # calculate the time difference after a chunk
      res <- difftime(Sys.time(), now, units = "mins")
      # return a character string to show the time
     paste("This chunk took", round(res, 2), "minutes")
    }
  }
}))

```

\newpage

# READ ME

This code processes the dataset 'processed_data_phyloseq.RDS', which is a phyloseq object that can be downloaded from <https://zenodo.org/deposit/5337076>. This dataset contains 16S microbiome data from 1109 meerkat faecal samples. The dataset is not normalised, but lab and sand contaminants have been filtered out, and only samples used in the final analysis are included. Original data contained an internal standard of Allobacillus and Imtechella - these have already been removed. Samples with low sequencing depth and any replicates have been removed. The code contains general data exploration and model validation that are additional to those presented in the paper. At the end it also analyzes the technical replicate data ('technical_replicate_data_phyloseq.RDS') and the pilot study data on captive meerkats ('pilot_study_data_phyloseq.RDS'). Note that this code takes several hours to run (~12 hours). To do a test run on a reduced dataset, there is the option to run a chunk after data import that generates a small but representative dataset, which takes about an hour to process, although some of the sensitivity analyses will throw up an error for this dataset due to small sample sizes. This Rmarkdown script allows errors to occur without halting the report. To turn this option off, remove "error = TRUE" from the global options. Each code chunk ends with elapsed time, and the total time taken to run the script is as the end of the document. For details on the study population and study design, please refer to the manuscript.

# Manuscript information

Diurnal oscillations in gut bacterial load and composition eclipse seasonal and lifetime dynamics in wild meerkats, Suricata suricatta

Risely, A.1\*, Wilhelm, K.1, Clutton-Brock, T.2, 3, 4, Manser, M.B. 3, 4, 5, and Sommer, S.1

1.  Institute for Evolutionary Ecology and Conservation Genomics, Ulm, Germany
2.  Large Animal Research Group, Department of Zoology, University of Cambridge, Downing Street, CB2 3EJ Cambridge, United Kingdom
3.  University of Pretoria, Mammal Research Institute, Pretoria, ZA
4.  Kalahari Research Trust, Kuruman River Reserve, Northern Cape, ZA
5.  Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zurich

Corresponding author: Alice Risely, riselya\@gmail.com

#### Abstract

Circadian rhythms in gut microbiota composition are crucial for metabolic function, yet the extent to which they govern microbial dynamics compared to seasonal and lifetime processes remains unknown. Here we investigate gut bacterial dynamics in wild meerkats (Suricata suricatta) over a 20-year period to compare diurnal, seasonal, and lifetime processes in concert, applying ratios of absolute abundance. We found that diurnal oscillations in bacterial load and composition eclipsed seasonal and lifetime dynamics. Diurnal oscillations were characterised by a peak in Clostridium abundance at dawn, were associated with temperature-constrained foraging schedules, and did not decay with age. Some genera exhibited seasonal fluctuations, whilst others developed with age, although we found little support for microbial senescence in very old meerkats. Strong microbial circadian rhythms in this species may reflect the extreme daily temperature fluctuations typical of arid-zone climates. Our findings demonstrate that accounting for circadian rhythms is essential for future gut microbiome research.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, out.width = '100%',  fig.align = "center", message = FALSE, warning = FALSE, fig.height = 4, fig.width =6, time_it = TRUE, error = TRUE)
```

\newpage

# R setup

#### Load packages

```{r, message = FALSE, warning = FALSE}

# data analysis packages
library(phyloseq)
library(mgcv)
library(gratia)
library(vegan)
library(tidyverse)
library(caret)
library(metagMisc)
library(expss)
library(scales)
library(broom)
library(lme4)
library(performance)

# viz packages
library(ggplot2)
library(magick)
library(RColorBrewer)
library(ggsci)
library(viridis)
library(ggpubr)
library(sjPlot)
library(gridExtra)

memory.limit(1000000) # set large memory limit

```

### Session info

```{r}
sessionInfo()

start_time<-Sys.time()
```

# Data import and normalisation

#### Import data and scale to internal standard

-   Data is in the format of a phyloseq object.
-   The scaling factor per sample is included as a metadata column called 'scaling factor'. This is used to normalise the data.
-   Varables are defined below.

```{r, message = FALSE, warning = FALSE}
meerkat_final<-readRDS("processed_data_phyloseq.RDS") # import data as phyloseq object

meerkat_final #data summary 

## scale to internal standard
asv_table<-data.frame(meerkat_final@otu_table@.Data)
names(asv_table)<-sample_data(meerkat_final)$feature.id
meerkat_normalized<-sweep(asv_table, MARGIN=2, sample_data(meerkat_final)$scaling_factor, `*`)
meerkat_normalized<-round(meerkat_normalized, digits = 0)
meerkat_normalized<-otu_table(meerkat_normalized, taxa_are_rows = TRUE)
meerkat_scaled<-meerkat_final
otu_table(meerkat_scaled)<-meerkat_normalized

### add total abundance of reads to dataframe
sample_data(meerkat_scaled)$TotalAbundance<-sample_sums(meerkat_scaled)

paste("Number of individuals sampled:",length(unique(sample_data(meerkat_scaled)$IndividID)))

# 1109 samples from 235 meerkats

## how many samples per individual?

summary(data.frame(table(sample_data(meerkat_scaled)$IndividID))$Freq)

```


#### Optional test dataset

To speed up analysis to see if the script works, run this code which generates a smaller dataset (~300) samples) to work with for testing (change "eval = FALSE" to "eval = TRUE" option at the start of this chunk to use this dataset). This test dataset needs to be quite specific for models to work - the samples need to be distributed across limited numbers of individuals and social groups so random effects do not break models. Some validation code might not work with this reduced dataset, therefore I have allowed errors so that these do not halt the script.

```{r testdata, eval = FALSE}

#not run

samplecount<-data.frame(sample_data(meerkat_scaled)) %>% 
  group_by(IndividID) %>% 
  summarise(NoSamples = n())

samplecount<-samplecount%>%arrange(IndividID, -NoSamples)

sample_data(meerkat_scaled)$SampleCount<-vlookup(sample_data(meerkat_scaled)$IndividID, samplecount, lookup_column = "IndividID", result_column = "NoSamples")


test_data<-subset_samples(meerkat_scaled, SampleCount >5)
samplenames<-sample_names(test_data)
subsamples<-sample(samplenames, 300) # choose 150

# also need some samples from sequencing run 4 to that some of the code works
run4<-subset_samples(meerkat_scaled, Seq_run == "RUN4")
samplenames_run4<-sample_names(run4)
subsamples_final<-c(subsamples, samplenames_run4)

test_data<-prune_samples(subsamples_final, meerkat_scaled)
test_data<-microbiome::core(test_data, detection = 10, prevalence =0.01)

# remove rare social groups
majorgroups<-data.frame(table(sample_data(test_data)$GroupAtSampling))
majorgroups<-majorgroups[order(-majorgroups$Freq),]
majorgroups<-subset(majorgroups, Freq >10)
majorgroups<-as.character(majorgroups$Var1)
sample_data(test_data)$Test2<-sample_data(test_data)$GroupAtSampling %in% majorgroups

test_data<-subset_samples(test_data, Test2==TRUE) 

meerkat_scaled<-test_data

```


# Description of sample metadata

-   Structure of metadata and all variables collated.

```{r, message = FALSE, warning = FALSE}
str(data.frame(sample_data(meerkat_final)))
```

[**VARIABLE DEFINITIONS**]{.ul}

[1] **IndividID** = Meerkat ID

[2] **feature.id** = DNA sample ID

[3] **SampleNo** = faecal sample ID

[4] **Storage** = Sample storage (freezedried or frozen)

[5] **Seq_run** = Sequencing run

[6] **Seq_depth_nospikein** = Sequencing depth after internal standard removed

[7] **HoursTilFrozen** = Number of hours sample was in field before frozen

[8] **scaling_factor** = Scaling factor to scale sample to internal standard

[9] **SampleDate** = Date sample was collected

[10] **SampleTime** = Time sample was collected

[11] **SampleTimeCat** = Categorical time variable (morning/afternoon)

[12] **HoursAfterSunrise** = Number of hours after sunrise at sample collection

[13] **month** = the month the sample was collected

[14] **Season** = Wet or dry season

[15] **Year** = Year sample was collected

[16] **AgeAtSamplingYears** = Meerkat age in years at sample collection

[17] **AgeCat** = Categorical age variable (young/adult/old)

[18] **Sex** = Meerkat sex

[19] **GroupAtSampling** = Social group at sampling

[20] **Condition_weight** = Residual body mass from GAM (mass ~ age)

[21] **SocialStatus** = Social status at time of sample collection (dom/sub)

[22] **HoursSinceForagingStart** = Time meerkat was foraging before sample collection

[23] **Temp_max** = Maximum temperature on day of sample collection

[24] **Temp_min** = Minimum temperature on day of sample collection

[25] **HourlyTemp** = Temperature (to the nearest hour) at sample collection

[26] **sum_rainfall_month** = Cumulative rainfall in mm over the previous month

# Figure 1: Sampling distribution and study design

-   Order meerkat ID by median sampling date for Fig. 1a.

```{r}
metadata <- data.frame(sample_data(meerkat_scaled)) #extract metadata as DF

# order meerkat ID my their median sampling date
metadata1 <- metadata %>%
  group_by(IndividID) %>%
  summarize(first = mean(SampleDate))

metadata1 <- metadata1[order(metadata1$first), ]
ID_order <- as.character(metadata1$IndividID)

# order ID by median sampling date
sample_data(meerkat_scaled)$IndividID <-
  factor(sample_data(meerkat_scaled)$IndividID, levels = ID_order) 
```

#### Figure 1a: Sample timeline

```{r, message = FALSE, warning = FALSE}

fig1a<-ggplot(data=sample_data(meerkat_scaled),aes(x=SampleDate,y=IndividID,group=IndividID))+
  geom_line(size=0.5,alpha=0.5)+
  geom_point(aes(fill=HoursAfterSunrise),colour="black",pch=21,size=1.5)+
  theme_bw(base_size=14)+
  theme(panel.grid.minor=element_blank(),
        panel.background=element_blank(),axis.line=element_line(colour="black"))+
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank())+
  scale_fill_viridis_c(option="plasma",direction=-1,name="Time of day")+
  scale_x_date(date_breaks = "1 year")+
  theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1))+
  ylab("Meerkat ID")+
  xlab("")+
  scale_y_discrete(expand=expansion(add=c(0.2,1.5)))+
  theme(legend.position="none")

```

#### Figure 1e: Sampling distributions across scales

-   Summarise sampling distribution across the three temporal scales.

```{r, message = FALSE, warning = FALSE}


hist1<-ggplot(sample_data(meerkat_scaled), aes(x = HoursAfterSunrise))+
  geom_histogram(fill = "grey", col = "black", binwidth = 1, boundary = 0)+
  theme_bw(base_size = 22)+
  xlab("Hours after sunrise")+
  ylab("Frequency")

hist2<-ggplot(sample_data(meerkat_scaled), aes(x = as.factor(month)))+
  geom_histogram(fill = "grey", col = "black", binwidth = 1, stat = "count")+
  theme_bw(base_size = 22)+
  xlab("Month")+
  ylab("Frequency")+
  scale_x_discrete(breaks = c("2","4", "6","8","10","12"))+
  theme(axis.title.y = element_blank())

hist3<-  ggplot(sample_data(meerkat_scaled), aes(x = AgeAtSamplingYears))+
geom_histogram(fill = "grey", col = "black", binwidth = 1, boundary = 0)+
  theme_bw(base_size = 22)+
  xlab("Age (years)")+
  ylab("Frequency")+
  scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10))+
  theme(axis.title.y = element_blank())


Fig1e<- ggarrange(hist1, hist2, hist3, ncol = 3, widths = c(1.2,1,1))
#ggarrange(hist1, hist2, hist3, ncol = 3, widths = c(1.2,1,1))

# save as pdf
#ggsave("FIGURES/fig1e.pdf")

```

#### Build figure 1

```{r figure1, fig.height=6, fig.width=12}
## import Fig1bcd

fig1bcd <- magick::image_read_pdf("C:/Users/risel/Dropbox/Sommer postdoc/Meerkat project/PROJECTS/MeerkatSpikeInData/Analysis updated/TEMPORAL DYNAMICS ANALYSIS/FIGURES/Fig1abc.pdf")

## give it a border

fig1bcd_border <- fig1bcd  %>%
  magick::image_border("black", "15 x 15")

## import fig1e that i saved earlier

fig1e <- magick::image_read_pdf("C:/Users/risel/Dropbox/Sommer postdoc/Meerkat project/PROJECTS/MeerkatSpikeInData/Analysis updated/TEMPORAL DYNAMICS ANALYSIS/FIGURES/fig1e.pdf")

# give border

fig1e_border <- fig1e %>%
  magick::image_border("black", "15 x 15")

# this figure won't wont fully work with reduced datset.

fig1a +annotation_raster(fig1bcd_border, xmin = as.Date("1997-01-01"), xmax = as.Date("2010-01-01"), ymin =  "ID 134", ymax = "ID 231") +
annotation_raster(fig1e_border, xmin = as.Date("2010-01-01"), xmax = as.Date("2020-01-01"), ymin =  "ID 39", ymax = "ID 61") 

# ggsave("FIGURES/Figure_1.pdf", width = 12, height = 6, units = c("in"))
```

\newpage

# Temporal dynamics of bacterial load, alpha and beta diversity

## Bacterial load

We model log10-transformed bacterial load (i.e. total 16S copy number per each scaled sample) against time of day, month, and meerkat age using Generalized Additive Mixed Models whilst accounting for methodological variables (storage method, sampling depth, sequencing run, time in field prior to freezing) as fixed linear effects, individual ID and social group as random effects, and including an auto-regression term to account for temporal autocorrelation (Wood 2017). The auto-regression term doesn't really change models but I included it anyway as a precaution.

### Fit GAMM predicting bacterial load

-   Extract metadata as dataframe.
-   Assess data distribution.
-   Scale variables around the mean and log-transform if neccessary.
-   Fit GAMM predicting load.
-   Load is log-transformed because fitting other distributions (Poisson/neg binomial/zero inflated neg binomial) performed worse (although overall results were the same).

```{r fig.height=3, fig.width=6, message=FALSE, warning=FALSE}
metadata<-data.frame(sample_data(meerkat_scaled))

### bacterial load distribution
ggplot(metadata, aes(x = TotalAbundance))+geom_histogram(fill = "grey", col = "black", boundary = 0) + ggtitle("Untransformed")+
ggplot(metadata, aes(x = log10(TotalAbundance)))+geom_histogram(fill = "grey", col = "black", boundary = 0) + ggtitle("Log10-transformed")
```

```{r load-model, fig.height=4, fig.width=8}
#scale variables
metadata$Seq_depth_nospikein<-as.numeric(scale(metadata$Seq_depth_nospikein))
metadata$HoursTilFrozen<-as.numeric(scale(metadata$HoursTilFrozen))


## GAMM MODEL

abundance_model <- mgcv::gamm(log10(TotalAbundance)~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)

# GAMM summary

sjPlot::tab_model(abundance_model$gam,  
                  show.stat = TRUE, 
                  show.se = TRUE, 
                  title = "Table S2: GAMM predicting bacterial load. Note that the estimate for smooth terms is actually EDF (estimated degrees of freedom), it is not an 'estimate'. EDF reflects the degree of non-linearity of a curve (1 = linear).",
                  dv.labels = "Bacterial load (log-transformed)",
                   pred.labels = c("Intercept", "Storage (frozen)", "Seq. run (2)", "Seq run (3)", "Seq. run (4)", "Time in field", "Hours after sunrise (non-linear)", "Month (non-linear)", "Age in years (non-linear)", "Seq. depth (non-linear)", "Meerkat ID (random)", "Social group (random)"))

# Temporal effects (without residuals)
gratia::draw(abundance_model$gam, select = c(1:3), ncol = 3, scales = "fixed")
```

#### Check model

```{r check-load-model, fig.height=8, fig.width=10}

model_check_load<-gratia::appraise(abundance_model$gam)

annotate_figure(model_check_load, top = text_grob("Model check of bacterial load")) # check model

```


#### Plot estimates for bacterial load

-   Extract confidence intervals per term from GAMM.
-   Generate rug data (for sample distribution layer).
-   Generate partial residuals dataframe per term.
-   Backtransform estimates into bacterial load.
-   Add lines where knots are placed

```{r, echo = FALSE, message = FALSE, warning = FALSE}
# function for tick breaks
int_breaks <- function(x, n = 5) {
  l <- pretty(x, n)
  l[abs(l %% 1) < .Machine$double.eps ^ 0.5] 
}
```

```{r load-fig, message = FALSE, warning = FALSE}
confint_df_time<-confint(abundance_model$gam, parm = "s(HoursAfterSunrise)",  type = "confidence")
confint_df_season<-confint(abundance_model$gam, parm = "s(month)",  type = "confidence")
confint_df_age<-confint(abundance_model$gam, parm = "s(AgeAtSamplingYears)",  type = "confidence")

confint_df_time$Variable<-"HoursAfterSunrise"
names(confint_df_time)[3]<-"Time"
confint_df_season$Variable<-"month"
names(confint_df_season)[3]<-"Time"
confint_df_age$Variable<-"AgeAtSamplingYears"
names(confint_df_age)[3]<-"Time"

confint_df_load<-rbind(confint_df_time, confint_df_season, confint_df_age)

###########add geom_rug data ##################

rugdata<- data.frame(cbind(metadata$HoursAfterSunrise, metadata$month, metadata$AgeAtSamplingYears))
names(rugdata)<-c("HoursAfterSunrise", "month", "AgeAtSamplingYears")
rugdata<-gather(rugdata, Variable, Time, HoursAfterSunrise:AgeAtSamplingYears, factor_key=TRUE)
rugdata$est = 0
rugdata$Variable<-factor(rugdata$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))

################# add residuals #######

fv <- predict(abundance_model$gam,type="terms") ## get term estimates

prsd_time <- data.frame(as.numeric(residuals(abundance_model$gam,type="working")) + fv[,3])
names(prsd_time)<-"Residuals"
prsd_time$Time<-metadata$HoursAfterSunrise
prsd_time$Variable<-"HoursAfterSunrise"
prsd_time$predicted<-predict(abundance_model$gam,type="response")

prsd_month <-data.frame(as.numeric(residuals(abundance_model$gam,type="working")) + fv[,4])
names(prsd_month)<-"Residuals"
prsd_month$Time<-metadata$month
prsd_month$Variable<-"month"
prsd_month$predicted<-predict(abundance_model$gam,type="response")

prsd_age <- data.frame(as.numeric(residuals(abundance_model$gam,type="working")) + fv[,5])
names(prsd_age)<-"Residuals"
prsd_age$Time<-metadata$AgeAtSamplingYears
prsd_age$Variable<-"AgeAtSamplingYears"
prsd_age$predicted<-predict(abundance_model$gam,type="response")

resids_df<-rbind(prsd_time, prsd_month, prsd_age)
resids_df$est<-resids_df$predicted+resids_df$Residuals

## rearrange level order 
resids_df$Variable<-factor(resids_df$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))
confint_df_load$Variable<-factor(confint_df_load$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))

#### make y axis the real abunance values, not estimates

#10^summary(abundance_model$gam)$p.coeff[1] # mean

confint_df_load$est<-confint_df_load$est+summary(abundance_model$gam)$p.coeff[1] 
confint_df_load$lower<-confint_df_load$lower+summary(abundance_model$gam)$p.coeff[1] 
confint_df_load$upper<-confint_df_load$upper+summary(abundance_model$gam)$p.coeff[1] 

confint_df_load$est<-10^confint_df_load$est
confint_df_load$lower<-10^confint_df_load$lower 
confint_df_load$upper<-10^confint_df_load$upper

resids_df$est<-10^resids_df$est


resids_load<-resids_df

####


P1_resids<-ggplot(subset(confint_df_load, Variable == "HoursAfterSunrise"), aes(x = Time, y = est))+
  theme_light(base_size = 14)+
  geom_point(data = subset(resids_load, Variable =="HoursAfterSunrise"),alpha = 0.5, aes(fill=Time), col = "white", pch = 21, size = 2)+
#  geom_line( size = 3,aes(col = Time))+
  geom_line( size = 1,col = "black", alpha = 0.7)+
  scale_color_viridis_c(option = "plasma", direction= -1)+
  scale_fill_viridis_c(option = "plasma", direction= -1)+
  geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.3)+
  theme( strip.background = element_blank(),
    strip.text.x = element_blank())+
  ylab("Bacterial abundance")+
  xlab("Hours after sunrise")+

  scale_x_continuous(breaks = int_breaks)+
  theme(legend.position = "none")+
  coord_cartesian(ylim = c(1000,1000000))+
  theme( panel.grid = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x)))+
  geom_vline(xintercept = c(as.numeric(abundance_model$gam$smooth[[1]]$xp)), col = "darkgrey", linetype = "dashed")


P2_resids<-ggplot(subset(confint_df_load, Variable == "month"), aes(x = Time, y = est))+
  theme_light(base_size = 14)+
  geom_jitter(data = subset(resids_load, Variable =="month"),alpha = 0.5, aes(fill=Time), col = "white", pch = 21, size = 2, width = 0.15)+
  
  geom_line( size = 1,col = "black", alpha = 0.7)+
  scale_color_gradientn(colours = c( "seagreen","seagreen","seagreen","seagreen","khaki","khaki","khaki","khaki","khaki","seagreen",  "seagreen"))+
  scale_fill_gradientn(colours = c( "seagreen","seagreen","seagreen","seagreen","khaki","khaki","khaki","khaki","khaki","seagreen",  "seagreen"))+
  geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.3)+
  theme( strip.background = element_blank(),
    strip.text.x = element_blank())+
  ylab("Model estimate")+
  xlab("Month")+
  coord_cartesian(ylim = c(1000,1000000))+
  theme(axis.title.y = element_blank(), 
    axis.text.y = element_blank())+
  scale_x_continuous(breaks = int_breaks)+
  theme(legend.position = "none")+
  # ylim(-2,1)+
  theme( panel.grid = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x)))+
  geom_vline(xintercept = c(as.numeric(abundance_model$gam$smooth[[2]]$xp)), col = "darkgrey", linetype = "dashed")


P3_resids<-ggplot(subset(confint_df_load, Variable == "AgeAtSamplingYears"), aes(x = Time, y = est))+
  theme_light(base_size = 14)+
  geom_point(data = subset(resids_load, Variable =="AgeAtSamplingYears"),alpha = 0.5, aes(fill=Time), col = "white", pch = 21, size = 2)+
  geom_line( size = 1,col = "black", alpha = 0.7)+
  scale_color_gradientn(colours = c( "blue","blue","darkgray","darkgray","darkgray","red", "red", "red", "red", "red", "red"))+
  scale_fill_gradientn(colours = c( "blue","blue","darkgray","darkgray","darkgray","red", "red", "red", "red", "red", "red"))+
  geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.3)+
  theme( strip.background = element_blank(),
    strip.text.x = element_blank())+
  ylab("Model estimate")+
  xlab("Age (years)")+
  coord_cartesian(ylim = c(1000,1000000))+
  theme(axis.title.y = element_blank(), 
    axis.text.y = element_blank())+
  scale_x_continuous(breaks = int_breaks)+
  theme(legend.position = "none")+
  # ylim(-2,1)+
  theme( panel.grid = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"))+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
    labels = trans_format("log10", math_format(10^.x)))+
  geom_vline(xintercept = c(as.numeric(abundance_model$gam$smooth[[3]]$xp)), col = "darkgrey", linetype = "dashed")

Fig2abc<- ggarrange(P1_resids, P2_resids, P3_resids, ncol = 3, align = "h", widths = c(1.2,1,1))

```

#### Partial R\^2 of terms

-   Remove variables and examine drop in R\^2.

```{r partial-r-load}

# just methods 

abundance_model_methods <- mgcv::gamm(log10(TotalAbundance)~
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen,
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)

paste("The partial R2 for methods is", summary(abundance_model_methods$gam)$r.sq)


## no methods
abundance_model_nomethods <- mgcv::gamm(log10(TotalAbundance)~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
   #    s(Seq_depth_nospikein, bs = "cr")+
    #Storage+
  #  Seq_run+
   # HoursTilFrozen,
    s(IndividID, bs = "re")+
   s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)



### minus time of day

load_r2_time<-mgcv::gamm(log10(TotalAbundance)~
   # s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)


paste("The partial R2 for time of day is", summary(abundance_model$gam)$r.sq - summary(load_r2_time$gam)$r.sq)

# minus month

load_r2_month<-mgcv::gamm(log10(TotalAbundance)~
    s(HoursAfterSunrise, bs = "cr")+
  #  s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)


paste("The partial R2 for month is", summary(abundance_model$gam)$r.sq - summary(load_r2_month$gam)$r.sq)


# minus age

load_r2_age<-mgcv::gamm(log10(TotalAbundance)~
    s(HoursAfterSunrise, bs = "cr")+
   s(month, bs = "cc")+
   # s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)


paste("The partial R2 for age is", summary(abundance_model$gam)$r.sq - summary(load_r2_age$gam)$r.sq)



## id and social group
abundance_model_ID <- mgcv::gamm(log10(TotalAbundance)~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
       s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen,
  #  s(IndividID, bs = "re")+
  #  s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)


paste("The partial R2 for meerkat ID is", summary(abundance_model$gam)$r.sq - summary(abundance_model_ID$gam)$r.sq)

```

### Sensitivity analyses for bacterial load (Fig. S10)

-   Select 20 samples per hour and repeat model to test for effects of unequal sampling distribution.

```{r sensitivity-load,  fig.height=10, fig.width=10}
metadata$TimeInterval <- cut(metadata$HoursAfterSunrise, breaks = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15))

table(metadata$TimeInterval)

# remove intervals with less than 20 samples

newdata<-subset(metadata, TimeInterval != "(6,7]" & TimeInterval != "(7,8]" & TimeInterval != "(14,15]" )
table(newdata$TimeInterval)
newdata$TimeInterval<-factor(newdata$TimeInterval)

# sample 20 rows per time inetrval
set.seed(999)
newdata1<-newdata%>% group_by(TimeInterval) %>% sample_n(20)

## repeat model


reduced_model_load<-mgcv::gamm(log10(TotalAbundance)~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen,
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=newdata1,
  family = gaussian)



figS10a<-draw(reduced_model_load$gam, select = c(1:3), ncol = 3, scales = "fixed") 
figS10a<-annotate_figure(figS10a, top = text_grob("a) 20 samples per hour", size = 16, face = "bold"))
figS10b<-draw(abundance_model_nomethods$gam, scales = "fixed", ncol = 3, select = c(1:3)) 
figS10b<-annotate_figure(figS10b, top = text_grob("b) Methods excluded", size = 16, face = "bold"))


ggarrange(figS10a, figS10b, ncol = 1)

```

#### Model validation

-   Create test data (70% of data) and training data (30%) and calculate R\^2.
-   Do this 100 times and get average R\^2.

```{r validatation-load, fig.height=6, fig.width=8}

# Test predictive power of model
# Create the training and test datasets

rsq_list<-list()

for(i in 1:100){
  # Step 1: Get row numbers for the training data
  trainRowNumbers <- createDataPartition(metadata$Storage, p=0.7, list=FALSE)
  
  # Step 2: Create the training  dataset
  trainData <- metadata[trainRowNumbers,]
  dim(trainData)
  
  # Step 3: Create the test dataset
  testData <- metadata[-trainRowNumbers,]
  dim(testData)
  
  # new model of train data
  
abundance_model_train <- mgcv::gamm(log10(TotalAbundance)~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=trainData,
  family = gaussian)
  
  
  # Predict on testData
  testData$predicted <- mgcv::predict.gam(abundance_model_train$gam, testData)
  ggplot(testData, aes(x = predicted, y = log10(TotalAbundance)))+geom_point()+geom_smooth(method = "lm")
  rsq_list[[i]]<-summary(lm(log10(TotalAbundance) ~ predicted, data =testData))$r.sq
  
}

abundance_rsqs<-do.call(c, rsq_list)

summary(abundance_rsqs)

# R2 distribution of 100 test models

# black line is mean for test datasets
# red line is true R2 from full dataset

ggplot()+aes(abundance_rsqs)+ 
  geom_histogram(fill = "grey", col = "black")+
  theme_light(base_size = 10)+
  xlab("R^2 distribution for test dataset - bacterial load")+
  ylab("Frequency")+
  geom_vline(xintercept = mean(abundance_rsqs), size = 2, linetype = "dashed")+
  geom_vline(xintercept = summary(abundance_model$gam)$r.sq , size = 2, linetype = "dashed", col = "red")

# one sided t-test 
t.test(abundance_rsqs, mu = summary(abundance_model$gam)$r.sq, alternative = "two.sided")

```

\newpage

## Alpha diversity

We model observed ASV richness against time of day, month, and meerkat age using Generalized Additive Mixed Models whilst accounting for methodological variables (storage method, sampling depth, and sequencing run) as fixed linear effects, individual ID and social group as random effects, and including a autoregression term to account for temporal autocorrelation (Wood 2017).

#### Estimating alpha diversity

-   Estimate ASV richness and Shannon diversity
-   Scale variables around the mean and log10 transform if necessary.

```{r alpha-calculate, fig.height=3, fig.width=6, message=FALSE, warning=FALSE}

# estimate alpha diversity
sample_data(meerkat_scaled)$Observed<-phyloseq::estimate_richness(meerkat_scaled, measures="Observed")
sample_data(meerkat_scaled)$Observed<-sample_data(meerkat_scaled)$Observed$Observed

sample_data(meerkat_scaled)$Shannon<-phyloseq::estimate_richness(meerkat_scaled, measures="Shannon")
sample_data(meerkat_scaled)$Shannon<-sample_data(meerkat_scaled)$Shannon$Shannon

## extract metadata as dataframe
metadata<-data.frame(sample_data(meerkat_scaled))

alpha_histograms<-ggarrange(ggplot(metadata, aes(x = Observed))+geom_histogram(fill = "grey", col = "black", boundary = 0) + ggtitle("ASV richness"),
  ggplot(metadata, aes(x = Shannon))+geom_histogram(fill = "grey", col = "black", boundary = 0) + ggtitle("Shannon"))

annotate_figure(alpha_histograms, top = text_grob("Distribution of alpha diversity measures", size = 16, face = "bold"))
                                                                                                      
# scale variables around the mean

metadata$Seq_depth_nospikein<-as.numeric(scale(metadata$Seq_depth_nospikein))
metadata$HoursTilFrozen<-as.numeric(scale(metadata$HoursTilFrozen))

```

### Fit GAMM to observed ASV richness

```{r alpha-gamm, fig.height=4, fig.width=8, message=FALSE, warning=FALSE}

m_alpha <- mgcv::gamm(Observed~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata)


sjPlot::tab_model(m_alpha$gam,  show.stat = TRUE, show.se = TRUE, title = "Table S3: GAMM predicting observed ASV richness. Note that the estimate for smooth terms is actually EDF (estimated degrees of freedom), it is not an 'estimate'. EDF reflects the degree of non-linearity of a curve (1 = linear)",
                   dv.labels = "Observed ASV richness",
                   pred.labels = c("Intercept", "Storage (frozen)", "Seq. run (2)", "Seq run (3)", "Seq. run (4)", "Time in field", "Hours after sunrise (non-linear)", "Month (non-linear)", "Age in years (non-linear)", "Seq. depth (non-linear)", "Meerkat ID (random)", "Social group (random)"))

draw(m_alpha$gam, select = c(1:3), ncol = 3, scales = "fixed")
```

#### Check model

```{r fig.height=8, fig.width=10}
model_check_alpha<-gratia::appraise(m_alpha$gam)

annotate_figure(model_check_alpha, top = text_grob("Model check of observed ASV richness")) # check model

```

#### Plot effects for each temporal variable with partial residuals

-   Extract confidence intervals for each term usinging gratia's confint() function.
-   Extract residuals using 'working' residuals
-   rbind into a dataframe.

```{r, message = FALSE, warning = FALSE}

# generate table of model estimates
confint_df_time<-confint(m_alpha$gam, parm = "s(HoursAfterSunrise)",  type = "confidence")
confint_df_season<-confint(m_alpha$gam, parm = "s(month)",  type = "confidence")
confint_df_age<-confint(m_alpha$gam, parm = "s(AgeAtSamplingYears)",  type = "confidence")

confint_df_time$Variable<-"HoursAfterSunrise"
names(confint_df_time)[3]<-"Time"
confint_df_season$Variable<-"month"
names(confint_df_season)[3]<-"Time"
confint_df_age$Variable<-"AgeAtSamplingYears"
names(confint_df_age)[3]<-"Time"

confint_df_alpha<-rbind(confint_df_time, confint_df_season, confint_df_age)
confint_df_alpha$Variable<-factor(confint_df_alpha$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))

### make y axis the real values, not estimates, by adding the model intercept (mean)
confint_df_alpha$est<-confint_df_alpha$est+summary(m_alpha$gam)$p.coeff[1] 
confint_df_alpha$lower<-confint_df_alpha$lower+summary(m_alpha$gam)$p.coeff[1] 
confint_df_alpha$upper<-confint_df_alpha$upper+summary(m_alpha$gam)$p.coeff[1] 

## parital residuals

fv <- predict(m_alpha$gam,type="terms") ## get term estimates

prsd_time <- data.frame(as.numeric(residuals(m_alpha$gam,type="working")) + fv[,3])
names(prsd_time)<-"Residuals"
prsd_time$Time<-metadata$HoursAfterSunrise
prsd_time$Variable<-"HoursAfterSunrise"
prsd_time$predicted<-predict(m_alpha$gam,type="response")

prsd_month <-data.frame(as.numeric(residuals(m_alpha$gam,type="working")) + fv[,4])
names(prsd_month)<-"Residuals"
prsd_month$Time<-metadata$month
prsd_month$Variable<-"month"
prsd_month$predicted<-predict(m_alpha$gam,type="response")


prsd_age <- data.frame(as.numeric(residuals(m_alpha$gam,type="working")) + fv[,5])
names(prsd_age)<-"Residuals"
prsd_age$Time<-metadata$AgeAtSamplingYears
prsd_age$Variable<-"AgeAtSamplingYears"
prsd_age$predicted<-predict(m_alpha$gam,type="response")


resids_df<-rbind(prsd_time, prsd_month, prsd_age)
resids_df$est<-resids_df$predicted+resids_df$Residuals

## rearrange level order 
resids_df$Variable<-factor(resids_df$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))

resids_alpha<-resids_df

```

-   Plot Figure 2b-d
-   Add lines where knots are placed

```{r alpha-plot-resids, message = FALSE, warning = FALSE}
alpha1_resids<-ggplot(subset(confint_df_alpha, Variable == "HoursAfterSunrise"), aes(x = Time, y = est))+
  theme_light(base_size = 14)+
  geom_point(data = subset(resids_alpha, Variable =="HoursAfterSunrise"),alpha = 0.5, aes(fill=Time), col = "white", pch = 21, size = 2)+
  geom_line( size = 1, col = "black", alpha = 0.7)+
  geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.3)+
  theme( strip.background = element_blank(),
    strip.text.x = element_blank())+
  ylab("Observed ASV richness")+
  xlab("Hours after sunrise")+
  ylim(0,400)+
  scale_x_continuous(breaks = int_breaks)+
  theme(legend.position = "none")+
  scale_color_viridis_c(option = "plasma", direction= -1)+
  scale_fill_viridis_c(option = "plasma", direction= -1)+
  theme( panel.grid = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"))+
    geom_vline(xintercept = c(as.numeric(m_alpha$gam$smooth[[1]]$xp)), col = "darkgrey", linetype = "dashed")


alpha2_resids<-ggplot(subset(confint_df_alpha, Variable == "month"), aes(x = Time, y = est))+
  theme_light(base_size = 14)+
  geom_jitter(data = subset(resids_alpha, Variable =="month"),alpha = 0.5, aes(fill=Time), col = "white", pch = 21, size = 2, width = 0.15)+
  geom_line( size = 1, col = "black", alpha = 0.7)+
  geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.3)+
  ylim(0,400)+
  xlab("Month")+
  theme( strip.background = element_blank(),
    strip.text.x = element_blank())+
  theme( axis.title.y= element_blank(), axis.text.y = element_blank())+
  scale_x_continuous(breaks = int_breaks)+
  theme(legend.position = "none")+
  scale_colour_gradientn(colours = c( "seagreen","seagreen","seagreen","seagreen","khaki","khaki","khaki","khaki","khaki","seagreen",  "seagreen"))+
  scale_fill_gradientn(colours = c( "seagreen","seagreen","seagreen","seagreen","khaki","khaki","khaki","khaki","khaki","seagreen",  "seagreen"))+
  theme( panel.grid = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"))+
      geom_vline(xintercept = c(as.numeric(m_alpha$gam$smooth[[2]]$xp)), col = "darkgrey", linetype = "dashed")



alpha3_resids<-ggplot(subset(confint_df_alpha, Variable == "AgeAtSamplingYears"), aes(x = Time, y = est))+
  theme_light(base_size = 14)+
  geom_point(data = subset(resids_alpha, Variable =="AgeAtSamplingYears"),alpha = 0.5, aes(fill=Time), col = "white", pch = 21, size = 2)+
  geom_line( size = 1, col = "black", alpha = 0.7)+
  geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.3)+
  ylim(0,400)+
  xlab("Age (years)")+
  theme( strip.background = element_blank(),
    strip.text.x = element_blank())+
  theme(axis.title.y= element_blank(), axis.text.y = element_blank())+
  scale_x_continuous(breaks = int_breaks)+
  theme(legend.position = "none")+
  scale_colour_gradientn(colours = c( "blue","blue","darkgray","darkgray","darkgray","darkgray", "red", "red", "red", "red", "red", "red"))+
  scale_fill_gradientn(colours = c( "blue","blue","darkgray","darkgray","darkgray","darkgray", "red", "red", "red", "red", "red", "red"))+
  theme( panel.grid = element_blank(),
    panel.background = element_blank(), axis.line = element_line(colour = "black"))+
      geom_vline(xintercept = c(as.numeric(m_alpha$gam$smooth[[3]]$xp)), col = "darkgrey", linetype = "dashed")



Fig2def<-ggarrange(alpha1_resids, alpha2_resids, alpha3_resids, ncol = 3, widths = c(1.2,1,1))

```

#### Partial R\^2 of terms

-   Remove variables and examine drop in R\^2.

```{r partial-r-alpha}

## methods

partial_r2_methods <- mgcv::gamm(Observed~
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen,
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata)

paste("The partial R2 for methods is", summary(partial_r2_methods$gam)$r.sq)



## no methods
partial_r2_nomethods <- mgcv::gamm(Observed~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(IndividID, bs = "re")+
   s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata,
  family = gaussian)


## Partial R2 time of day

partial_r2_time <- mgcv::gamm(Observed~
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata)

paste("The partial R2 for time of day is", summary(m_alpha$gam)$r.sq - summary(partial_r2_time$gam)$r.sq)


## Partial R2 month

partial_r2_month <- mgcv::gamm(Observed~
    s(HoursAfterSunrise, bs = "cr")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata)

paste("The partial R2 for month is", summary(m_alpha$gam)$r.sq - summary(partial_r2_month$gam)$r.sq)

## Partial R2 age

partial_r2_age <- mgcv::gamm(Observed~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata)

paste("The partial R2 for age is", summary(m_alpha$gam)$r.sq - summary(partial_r2_age$gam)$r.sq)

## id and social group

partial_r2_ID <- mgcv::gamm(Observed~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen,
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata)

paste("The partial R2 for ID and social group is", summary(m_alpha$gam)$r.sq - summary(partial_r2_ID$gam)$r.sq)


```

### Sensitivity analyses for alpha diversity (Fig. S11)

-   Select 20 samples per hour and repeat model to test for effects of unequal sampling distribution.

```{r subsampling-alpha, fig.height=8, fig.width=10}

metadata$TimeInterval <- cut(metadata$HoursAfterSunrise, breaks = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15))

table(metadata$TimeInterval)
unique(metadata$TimeInterval)

newdata<-subset(metadata, TimeInterval != "(6,7]" & TimeInterval != "(7,8]" & TimeInterval != "(14,15]" )
table(newdata$TimeInterval)
newdata$TimeInterval<-factor(newdata$TimeInterval)

set.seed(999)
newdata1<-newdata%>% group_by(TimeInterval) %>% sample_n(20)

## repeat model

dim(newdata1)

reduced_model_alpha<-mgcv::gamm(Observed~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen,
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=newdata1,
  family = gaussian)

figS11a<-draw(reduced_model_alpha$gam, select = c(1:3), ncol = 3, scales = "fixed") # Fig S11a
figS11a<-annotate_figure(figS11a, top = text_grob("a) 20 samples per hour", size = 14, face = "bold"))

figS11b<-draw(partial_r2_nomethods$gam, scales = "fixed", ncol = 3, select = c(1:3)) # Fig S11b
figS11b<-annotate_figure(figS11b, top = text_grob("b) Methods excluded", size = 14, face = "bold"))


ggarrange(figS11a, figS11b, ncol = 1)

```

#### Model validation

-   Create test data (70% of data) and training data (30%) and calculate R\^2.
-   Do this 100 times and get average.

```{r validation-alpha, fig.height=6, fig.width=8}

# Test predictive power of model
# Create the training and test datasets

rsq_list<-list()

for(i in 1:100){
  
  # Step 1: Get row numbers for the training data
  trainRowNumbers <- createDataPartition(metadata$Storage, p=0.7, list=FALSE)
  
  # Step 2: Create the training  dataset
  trainData <- metadata[trainRowNumbers,]
  
  # Step 3: Create the test dataset
  testData <- metadata[-trainRowNumbers,]
  dim(testData)
  
  m_alpha_train <- mgcv::gamm(Observed~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corAR1(), #residual correlation structure
  data=trainData)
  
  
  # Predict on testData
  testData$predicted <- mgcv::predict.gam(m_alpha_train$gam, testData)
  ggplot(testData, aes(x = predicted, y = Observed))+geom_point()+geom_smooth(method = "lm")
  rsq_list[[i]]<-summary(lm(Observed ~ predicted, data =testData))$r.sq
  
}

alpha_rsqs<-do.call(c, rsq_list)

summary(alpha_rsqs)

# black line is mean R2 for test data
# red line is true R2 from all data

ggplot()+aes(alpha_rsqs)+ 
  geom_histogram(fill = "grey", col = "black")+
  theme_light(base_size = 10)+
  xlab("R square distribution")+
  ylab("Frequency")+
  geom_vline(xintercept = mean(alpha_rsqs), size = 2, linetype = "dashed")+
  geom_vline(xintercept = summary(m_alpha$gam)$r.sq , size = 2, linetype = "dashed", col = "red")

# one sided t-test 
t.test(alpha_rsqs, mu = summary(m_alpha$gam)$r.sq, alternative = "two.sided")
```

### Fitting GAMM to Shannon diversity

```{r fig.height=4, fig.width=8}
m_shannon <- mgcv::gamm(Shannon~
    s(HoursAfterSunrise, bs = "cr")+
    s(month, bs = "cc")+
    s(AgeAtSamplingYears, bs = "cr")+
    s(Seq_depth_nospikein, bs = "cr")+
    Storage+
    Seq_run+
    HoursTilFrozen+
    s(IndividID, bs = "re")+
    s(GroupAtSampling, bs = "re"),
  correlation = corARMA(form = ~ 1|Year, p = 3),
  data=metadata)

## model summary

sjPlot::tab_model(m_shannon$gam,  show.stat = TRUE, show.se = TRUE, title = "Extra table: GAMM predicting Shannon diversity. Note that the estimate for smooth terms is actually EDF (estimated degrees of freedom), it is not an 'estimate'. EDF reflects the degree of non-linearity of a curve (1 = linear)",
                   dv.labels = "Shannon diversity",
                   pred.labels = c("Intercept", "Storage (frozen)", "Seq. run (2)", "Seq run (3)", "Seq. run (4)", "Time in field", "Hours after sunrise (non-linear)", "Month (non-linear)", "Age in years (non-linear)", "Seq. depth (non-linear)", "Meerkat ID (random)", "Social group (random)"))

# model estimates (plot)
shannon_effects<-draw(m_shannon$gam, select = c(1:3), ncol = 3, scales = "fixed")
annotate_figure(shannon_effects, top = text_grob("GAMM estimates for Shannon diversity", size = 14, face = "bold")) # check model


```

\newpage

## Beta diversity

We next explore shifts in beta diversity across the three temporal scales by ordinating taxa composition using a Multi-Dimensional Scaling (MDS) using Weighted Unifrac distance. We also present results for unweighted unifrac, and results of PERMANOVAS.

#### Ordinate with Weighted Unifrac

-   Filter dataset to only ASVs with over 50 reads
-   Ordinate and check relative contribution of axes

```{r beta-diversity, message = FALSE, warning = FALSE}

scaled_core<-microbiome::core(meerkat_scaled, detection = 50, prevalence = 0)

paste("This filtering step removed", 1-round(sum(sample_sums(scaled_core))/sum(sample_sums(meerkat_scaled)), 3), "% of reads")

```

```{r}
ord_wunifrac_mds<-ordinate(scaled_core, "MDS", "wunifrac")
ord_unifrac_mds<-ordinate(scaled_core, "MDS", "unifrac")

dist_wunifrac<-distance(scaled_core, method="wunifrac")
dist_unifrac<-distance(scaled_core, method="unifrac")

```

#### Axis contributions to variation

```{r fig.height=3, fig.width=6}
qplot(y = ord_wunifrac_mds$values$Relative_eig[1:10], x = factor(1:10),  geom="col", ylab = "Relative eingenvalue", xlab = "Axis",  main = "Scree plot")+ ggtitle("Weighted Unifrac axis contributions")+ylim(0,0.65)+
qplot(y = ord_unifrac_mds$values$Relative_eig[1:10], x = factor(1:10),  geom="col", ylab = "Relative eingenvalue", xlab = "Axis",  main = "Scree plot")+ ggtitle("Unweighted Unifrac axis contributions")+ylim(0,0.65)

```

#### Plotting beta diversity

-   Add MDS scores to sample metadata.
-   Split continous variables into categories for plotting.
-   For each plot, generate centroid ('cent') and segement ('segs') dataframes.

```{r fig.height=9, fig.width=9, warning=FALSE, , message=FALSE}
### ADD COLUMNS TO METADATA

sample_data(meerkat_scaled)$MDS1<-data.frame(ord_wunifrac_mds$vectors)$Axis.1
sample_data(meerkat_scaled)$MDS2<-data.frame(ord_wunifrac_mds$vectors)$Axis.2
sample_data(meerkat_scaled)$MDS3<-data.frame(ord_wunifrac_mds$vectors)$Axis.3
sample_data(meerkat_scaled)$MDS4<-data.frame(ord_wunifrac_mds$vectors)$Axis.4


## unweighted unifrac

sample_data(meerkat_scaled)$MDS1_UU<-data.frame(ord_unifrac_mds$vectors)$Axis.1
sample_data(meerkat_scaled)$MDS2_UU<-data.frame(ord_unifrac_mds$vectors)$Axis.2


### convert to dataframe 

metadata<-data.frame(sample_data(meerkat_scaled))


########

cent <- aggregate(cbind(MDS1, MDS2) ~ SampleTimeCat, data = metadata, FUN = mean)

beta1<-ggplot(metadata, aes(x = MDS1, y = MDS2)) +
  theme_bw(base_size = 14)+
  stat_ellipse(geom = "polygon",aes(col = SampleTimeCat), alpha = 0, size = 1)+
  geom_point(aes(fill = HoursAfterSunrise), pch = 21, size = 2, alpha = 0.8, col = "white") + # spiders
  #geom_point(aes(fill = HoursAfterSunrise), pch = 21, size = 1.5) + 
  scale_fill_viridis_c(option = "plasma", direction= -1)+
  scale_color_manual(values = c("darkmagenta", "darkgoldenrod1"))+
  geom_point(data = cent, size = 6, aes(col = SampleTimeCat), stroke = 1, alpha = 0.8) +  
  geom_point(data = cent, size = 6, col = "black", stroke = 1, pch = 21, alpha = 0.8) +  
  theme(legend.position = "none")+
  theme(plot.margin = unit(c(0.2,0.3,0.3,0.4), "cm"))+
  ylab("Axis 2 (17.1%)")+
  xlab("Axis 1 (25.8%)")+
    theme( panel.grid = element_blank())


#### season
cent <- aggregate(cbind(MDS1, MDS2) ~ Season, data = metadata, FUN = mean)


beta2<-ggplot(metadata, aes(x = MDS1, y = MDS2)) +
  theme_bw(base_size = 14)+
    stat_ellipse(geom = "polygon",aes(col = Season), alpha = 0, size = 1)+
 geom_point(aes(fill = month), pch = 21, size = 2, alpha = 0.8, col = "white") + 
  scale_fill_gradientn(colours = c( "seagreen","seagreen","seagreen","seagreen","khaki","khaki","khaki","khaki","khaki","seagreen",  "seagreen"))+
  # geom_point(aes(fill = Season), pch = 21, size = 1.5) + 
#  scale_fill_manual(values = c( "seagreen","khaki" ))+
  scale_color_manual(values = c( "seagreen","khaki"))+
 geom_point(data = cent, size = 6, aes(col = Season), stroke = 1, alpha = 0.8) +  
  geom_point(data = cent, size = 6, col = "black", stroke = 1, pch = 21, alpha = 0.7) +   
  theme(axis.text.y = element_blank(), axis.title.y = element_blank())+
  theme(legend.position = "none")+
   xlab("Axis 1 (25.8%)")+
    theme( panel.grid = element_blank())


#### age
cent <- aggregate(cbind(MDS1, MDS2) ~ AgeCat, data = metadata, FUN = mean)


beta3<-ggplot(metadata, aes(x = MDS1, y = MDS2)) +
  theme_bw(base_size = 14)+
      stat_ellipse(geom = "polygon",aes(col = AgeCat), alpha = 0, size = 1)+
 geom_point(aes(fill = AgeAtSamplingYears),  pch = 21, size = 2, alpha = 0.8, col = "white") + # spiders
  #geom_point(aes(fill = AgeCat), pch = 21, size = 1.5) + 
  scale_fill_gradientn(colours = c( "blue","blue","darkgray","darkgray","darkgray","darkgray", "red", "red", "red", "red", "red", "red"))+
  scale_color_manual(values = c( "red", "blue", "darkgrey"))+
 geom_point(data = cent, size = 6, aes(col = AgeCat), stroke = 1, alpha = 0.8) +  
  geom_point(data = cent, size = 6, col = "black", stroke = 1, pch = 21, alpha = 0.7) +   
  theme(axis.text.y = element_blank(), axis.title.y = element_blank())+
  theme(legend.position = "none")+
   xlab("Axis 1 (25.8%)")+
    theme( panel.grid = element_blank())


Fig2ghi<-ggarrange(beta1, beta2, beta3, ncol = 3, widths = c(1.2,1,1))

```

#### Unweighted unifrac

```{r fig.height=4, fig.width=12}

cent <- aggregate(cbind(MDS1_UU, MDS2_UU) ~ SampleTimeCat, data = metadata, FUN = mean)

beta1_UU<-ggplot(metadata, aes(x = MDS1_UU, y = MDS2_UU)) +
  theme_bw(base_size = 12)+
  stat_ellipse(geom = "polygon",aes(col = SampleTimeCat), alpha = 0, size = 1)+
  geom_point(aes(fill = HoursAfterSunrise), pch = 21, size = 2, alpha = 0.8, col = "white") + # spiders
  #geom_point(aes(fill = HoursAfterSunrise), pch = 21, size = 1.5) + 
  scale_fill_viridis_c(option = "plasma", direction= -1)+
  scale_color_manual(values = c("darkmagenta", "darkgoldenrod1"))+
  geom_point(data = cent, size = 6, aes(col = SampleTimeCat), stroke = 1, alpha = 0.8) +  
  geom_point(data = cent, size = 6, col = "black", stroke = 1, pch = 21, alpha = 0.8) +  
  #theme(legend.position = "none")+
  theme(plot.margin = unit(c(0.2,0.3,0.3,0.4), "cm"))+
  ylab("Axis 2 (9%)")+
  xlab("Axis 1 (16%)")+
    theme( panel.grid = element_blank())


#### season
cent <- aggregate(cbind(MDS1_UU, MDS2_UU) ~ Season, data = metadata, FUN = mean)


beta2_UU<-ggplot(metadata, aes(x = MDS1_UU, y = MDS2_UU)) +
  theme_bw(base_size = 12)+
    stat_ellipse(geom = "polygon",aes(col = Season), alpha = 0, size = 1)+
 geom_point(aes(fill = month), pch = 21, size = 2, alpha = 0.8, col = "white") + 
  scale_fill_gradientn(colours = c( "seagreen","seagreen","seagreen","seagreen","khaki","khaki","khaki","khaki","khaki","seagreen",  "seagreen"))+
  scale_color_manual(values = c( "seagreen","khaki"))+
 geom_point(data = cent, size = 6, aes(col = Season), stroke = 1, alpha = 0.8) +  
  geom_point(data = cent, size = 6, col = "black", stroke = 1, pch = 21, alpha = 0.7) +   
  theme(axis.text.y = element_blank(), axis.title.y = element_blank())+
 # theme(legend.position = "none")+
   xlab("Axis 1 (16%)")+
    theme( panel.grid = element_blank())


#### age
cent <- aggregate(cbind(MDS1_UU, MDS2_UU) ~ AgeCat, data = metadata, FUN = mean)


beta3_UU<-ggplot(metadata, aes(x = MDS1_UU, y = MDS2_UU)) +
  theme_bw(base_size = 12)+
      stat_ellipse(geom = "polygon",aes(col = AgeCat), alpha = 0, size = 1)+
 geom_point(aes(fill = AgeAtSamplingYears),  pch = 21, size = 2, alpha = 0.8, col = "white") + # spiders
  scale_fill_gradientn(colours = c( "blue","blue","darkgray","darkgray","darkgray","darkgray", "red", "red", "red", "red", "red", "red"))+
  scale_color_manual(values = c( "red", "blue", "darkgrey"))+
 geom_point(data = cent, size = 6, aes(col = AgeCat), stroke = 1, alpha = 0.8) +  
  geom_point(data = cent, size = 6, col = "black", stroke = 1, pch = 21, alpha = 0.7) +   
  theme(axis.text.y = element_blank(), axis.title.y = element_blank())+
  #theme(legend.position = "none")+
   xlab("Axis 1 (16%)")+
    theme( panel.grid = element_blank())

ggarrange(beta1_UU, beta2_UU, beta3_UU, ncol = 3, widths = c(1.2,1,1))

```

### Marginal PERMANOVAS

```{r}

metadata<-data.frame(sample_data(scaled_core))
dim(metadata)

perm_wunifrac_marginal<-adonis2(dist_wunifrac~
    Storage+
    Seq_run+
    HoursTilFrozen+
    Seq_depth_nospikein + 
    month+
    AgeAtSamplingYears+
    HoursAfterSunrise+
      GroupAtSampling+
      IndividID, 
  data = metadata,
  by = "margin")


perm_unifrac_marginal<-adonis2(dist_unifrac~
    Storage+
    Seq_run+
    HoursTilFrozen+
    Seq_depth_nospikein + 
    month+
    AgeAtSamplingYears+
    HoursAfterSunrise+
      GroupAtSampling+
      IndividID, 
  data = metadata,
  by = "margin")


```

#### Table S4: Beta diversity permanovas

```{r}

# weighted unifrac
perm_wunifrac_marginal

# unweighted unifrac
perm_unifrac_marginal

```

# Figure 2: Temporal dynamics of load, alpha diversity, and beta diversity

```{r fig.height=10, fig.width=12}

ggarrange(Fig2abc,Fig2def,Fig2ghi ,ncol = 1)

#ggsave("FIGURES/figure2.pdf",width = 9, height = 9, units = c("in"))

```

# Figure 3: Identifying genera driving temporal shifts in beta diversity

-   Aim is to generate barplot and colour ordinations by dominant genera/enterotype per sample
-   First we need to replace missing genera names with 'Unknown [family]' so that these are not excluded.
-   Merge ASVs by genus
-   Add column for dominant genus per sample

```{r}

# replace unknown genera with eg 'unknown Lachnospiraceae genus'
taxtable<-data.frame(scaled_core@tax_table@.Data)
taxtable$Family<-as.character(taxtable$Family)
taxtable$Genus<-as.character(taxtable$Genus)

taxtable$Genus<-ifelse(is.na(taxtable$Genus) , paste("Unknown", taxtable$Family,  sep = " "), taxtable$Genus)

taxtable$Genus<-ifelse(taxtable$Genus== "uncultured" , paste("Unknown", taxtable$Family,  sep = " "), taxtable$Genus)
taxtable$Genus<-ifelse(taxtable$Genus== "uncultured bacterium" , paste("Unknown", taxtable$Family, sep = " "), taxtable$Genus)
row.names(taxtable)<-taxa_names(scaled_core)
taxtable<-as.matrix(taxtable)
tax_table(scaled_core)<-taxtable

#scaled_core<-tip_glom(scaled_core, h = 0.2)
scaled_core_genus<-speedyseq::tax_glom(scaled_core, taxrank = "Genus")

#### remove some unknown or strange genera

scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "Unknown uncultured bacterium")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured organism")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured bacterium")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "Unknown uncultured")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "gut metagenome")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured Chloroflexi bacterium")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "metagenome")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured Firmicutes bacterium")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured rumen bacterium")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured soil bacterium")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "Unknown Unknown Family")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "BCf9-17 termite group")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "Blattella germanica (German cockroach)")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "CPla-4 termite group" )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "DSSF69" )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "E1B-B3-114"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "Family XIII AD3011 group"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "FFCH7168" )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "GCA-900066225"   )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "hydrothermal vent metagenome"   )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "Incertae Sedis"   )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "M2PB4-61 termite group" )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=   "OLB13"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "OM27 clade"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "OM27 clade"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "CHKCI001"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "CHKCI001"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "Termite planctomycete cluster" )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "uncultured Lachnospiraceae bacterium")
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "Unknown gir-aah93h0" )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "Unknown JG30-KF-CM45"  )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "Unknown Rs-E47 termite group"   )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus != "UBA1819"   )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "ZOR0006"   )
scaled_core_genus<-subset_taxa(scaled_core_genus, Genus !=  "GCA-900066575" )

taxtable<-data.frame(scaled_core_genus@tax_table@.Data)$Genus

taxa_names(scaled_core_genus)<-taxtable

#taxa_names(scaled_core_genus)=="Unknown Bacillaceae"

```

-   Add dominant taxa as column

```{r}

dominant_phylo<-scaled_core_genus
sample_data(dominant_phylo)<-sample_data(dominant_phylo)[,"feature.id"]

# melt data
psmelt_df<-speedyseq::psmelt(dominant_phylo)
head(psmelt_df)

dominant_taxa_df<-psmelt_df %>% 
  group_by(feature.id) %>% 
  arrange(feature.id, -Abundance) %>% 
  dplyr::mutate(rank=rank(-Abundance))

dominant_taxa_df<-subset(dominant_taxa_df, rank ==1)
dominant_taxa_df$Genus<-factor(dominant_taxa_df$Genus)

# add to metadata 
sample_data(scaled_core_genus)$DomGenus <-vlookup(sample_data(scaled_core_genus)$feature.id, dominant_taxa_df, lookup_column = "Sample", result_column = "Genus")

## colour by focus genera (we justify this list later)
# but for now want the colours to match later figures of the focus genera

top_genera<-c( "Clostridium sensu stricto 1" ,
"Blautia" ,
"Unknown Lachnospiraceae",
 "Enterococcus"  ,
"[Ruminococcus] torques group",
"Bacteroides",  
"Paeniclostridium",
"Romboutsia" ,
"Lachnoclostridium",  
"Lactococcus"  ,
"Peptococcus",
"Helicobacter" ,
 "Unknown Bacillaceae",
 "Raoultibacter"  ,
 "Cellulomonas", 
"Eubacterium"  ,
"Other" )

sample_data(scaled_core_genus)$DomGenus_plot<-ifelse(sample_data(scaled_core_genus)$DomGenus %in% top_genera, as.character(sample_data(scaled_core_genus)$DomGenus), "Other")
sample_data(scaled_core_genus)$DomGenus_plot<-factor(sample_data(scaled_core_genus)$DomGenus_plot, levels =top_genera)

```

-   Generate colour palette

```{r}

## colours scale

mypaly<-brewer.pal(12,"Paired")

mypaly[16]<-"darkgreen"
mypaly[5]<-"red"
mypaly[8]<-"darkorange"
mypaly[7]<-"darkblue"
mypaly[9]<-"magenta"
mypaly[10]<-"cyan"
mypaly[11]<-"#ffe20b"
mypaly[6]<-"#ae0001"
mypaly[12]<-"#795548"
mypaly[13]<-"peachpuff"
mypaly[14]<-"blue"
mypaly[15]<-"mediumorchid"
mypaly[15]<-"#8B008B"
mypaly[4]<-"forestgreen"
mypaly[17]<-"grey"

mypaly<-mypaly[c(1,4,3,2,5:8, 10,9,11:17)]

```

-   Barplot
-   Merge samples by half hour time intervals

```{r}

#Merge by genus
meerkat_genus <- scaled_core_genus %>%
 speedyseq::tax_glom(taxrank = "Genus") %>%                     # agglomerate at phylum level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt() %>%                                         # Melt to long format
  arrange(Genus)                                      # Sort data frame alphabetically by family

head(meerkat_genus)
dim(meerkat_genus)

# generate column for time interval
meerkat_genus$TimeInterval<-cut(meerkat_genus$HoursAfterSunrise, breaks = c(0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12,12.5,13,13.5,14,14.5))

#head(meerkat_genus)
meerkat_genus$Genus_plot<-ifelse(meerkat_genus$Genus %in% top_genera == TRUE, as.character(meerkat_genus$Genus), "Other")


meerkat_genus$Genus_plot<-factor(meerkat_genus$Genus_plot, levels = top_genera)

# merge
merged_samples<-meerkat_genus %>% group_by(TimeInterval, Genus_plot) %>% dplyr::summarize(mean_rel_abundance = sum(Abundance))

# barplot by half an hour intervals

barplot_by_time<- ggplot( merged_samples, aes(x = TimeInterval, y = mean_rel_abundance, fill = Genus_plot)) + 
  geom_col( width = 1, position = "fill", col = "black") +
    scale_fill_manual(values = mypaly)+
  ylab("Relative abundance")+
  xlab("30 minute time intervals after sunrise")+
  theme_light(base_size = 14)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position = "none")+
   theme(plot.margin=unit(c(0,0.2,0.2,0.2),"cm"))

sample_data(meerkat_scaled)$TimeInterval<- cut(sample_data(meerkat_scaled)$HoursAfterSunrise, breaks = c(0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12,12.5,13,13.5,14,14.5))

hist_plot<-data.frame(table(sample_data(meerkat_scaled)$TimeInterval))

histogram<-ggplot(hist_plot, aes(x = Var1, y = Freq))+
  geom_bar(stat = "identity")+
  geom_text(aes(label=Freq), position=position_dodge(width=0.9), vjust=-0.25)+
  theme_bw(base_size = 14)+
  theme(axis.text.x = element_blank())+
  theme(axis.title.x = element_blank())+
  ylab("No. samples")+
  ylim(0,150)+
   theme(plot.margin=unit(c(0.5,0.2,0,0.2),"cm"))

Fig3_2<-ggarrange(histogram, barplot_by_time, ncol = 1, heights = c(0.3,1), 
                  legend = FALSE, labels = c("a)"),   font.label = list(size = 20), hjust = 1)

```

### Statistics using envfit()

-   Keep only relevant biological variables
-   Use envfit() to model centroids using continuous variables for time of day and age
-   Month has to be categorical (dry/wet) since it can't have a linear relationship with microbiome
-   For plotting, use categorical variables so the direction of the effects are clear
-   Extract vector centroids
-   Keep only variables that are significant/have detectable effect sizes
-   Plot ordination with arrows pointing towards centroids of variables of interest

#### Axes 1 and 2 using continuous variables

```{r arrows, message = FALSE, warning = FALSE}
metadata<-data.frame(sample_data(scaled_core_genus))

vectors_df<-ord_wunifrac_mds$vectors

## model using continuous variables for time of day and age (to get partial R2 for each)

# axes 1 and 2 using continuous variables
(fit_continuous <- envfit(vectors_df, metadata[,c("HoursAfterSunrise",  "AgeAtSamplingYears", "Season",  "Storage",  "Seq_run", "Seq_depth_nospikein")], perm = 999, na.rm = TRUE,  choices=c(1,2)))
```

#### Axes 3 and 4 using continuous variables

```{r}
# axes 3 and 4 nusing continuous variables

(fit_continuous <- envfit(vectors_df, metadata[,c("HoursAfterSunrise",  "AgeAtSamplingYears", "Season",  "Storage",  "Seq_run", "Seq_depth_nospikein")], perm = 999, na.rm = TRUE,  choices=c(3,4)))
```


* For plots we use categorical variables for arrows since this is more intuitive

```{r}
metadata<-metadata[,c("SampleTimeCat",  "AgeCat", "Season",  "Storage",  "Seq_run", "Seq_depth_nospikein")]
```

#### Axes 1 and 2 using categorical variables

```{r}

fit <- envfit(vectors_df, metadata, perm = 999, na.rm = TRUE,  choices=c(1,2))
df<-data.frame(vectors_df)
df$Dom_taxa<-sample_data(scaled_core_genus)$DomGenus_plot
df$TotalAbundance<-sample_data(scaled_core_genus)$TotalAbundance

#Get the vectors for bioenv.fit
df_biofit<-vegan::scores(fit,display=c("factors"), choices=1:2)

## scaling factor for arrows to fill 30% of plot
mul <- ordiArrowMul(df_biofit, fill = 0.3)
df_biofit <- df_biofit * mul                    # Scale the biplot scores
df_biofit<-as.data.frame(df_biofit)


# remove methodological variables from biofit dataframe, 
#as I am just interested in biological variables after controling for methods
#head(df_biofit)


df_biofit<- df_biofit[c("SampleTimeCatAFTERNOON","SampleTimeCatMORNING"),]
row.names(df_biofit)<-c("Afternoon", "Morning")



### spider plot
cent <- aggregate(cbind(Axis.1, Axis.2) ~ Dom_taxa, data = df, FUN = mean)

segs <- merge(df, setNames(cent, c('Dom_taxa','oNMDS1','oNMDS2')),
  by = 'Dom_taxa', sort = FALSE)

df1<-df
df_biofit1<-df_biofit
cent1<-cent
segs1<-segs

# order taxa so that they at least slightly match later figures!
df1$Dom_taxa<-factor(df1$Dom_taxa,levels=top_genera)
segs1$Dom_taxa<-factor(segs1$Dom_taxa,levels=top_genera)


P1<-ggplot(df1, aes(x = Axis.1, y = Axis.2)) +
  theme_light(base_size = 14)+
 geom_point(aes(fill = Dom_taxa), size = 3, pch = 21, alpha = 0.8, col = "white") + # spiders
  scale_fill_manual(values = mypaly[c(-16)])+
  scale_color_manual(values = mypaly[c(-16)])+
  geom_jitter(data = cent1, size = 7, aes(fill = Dom_taxa), col = "black", pch = 21,  stroke = 1.5, width = 0.02, height = 0.02) +  
  geom_segment(data=df_biofit1, aes(x = 0, y = 0, xend = Axis.1*0.7, yend = Axis.2*0.7),
    arrow = arrow(length = unit(0.5, "cm"), type = "closed"), lwd = 1,  col = "black")+
  geom_label(data=df_biofit1,  
    alpha = 0.9, col = "black", size = 4, fill = "yellow", 
    hjust = c(0,1), 
    vjust = c(1,1),
    aes(Axis.1*0.7,Axis.2*0.7, label = rownames(df_biofit1)))+
  xlab("Axis 1 (25.8%)")+
  ylab("Axis 2 (17.1%)")+
  guides(fill = guide_legend(override.aes = list(size = 6), ncol = 5))+
  theme(legend.title = element_blank())+
  annotate("text", x = 0, y = 0, label = "X", size = 7)+
  theme(plot.margin = margin(0.2,0.2,1,0.2, "cm"))+
   annotate(geom="text", x=-0.3, y=-0.27, label="Hours after sunrise: R^2 = 0.27, p < 0.001",
             color="black", size = 5, hjust = 0)+
  xlim(-0.3,0.3)

```

#### Axes 3 and 4 using categorical variables

```{r fig.height=12, fig.width=10}

metadata<-data.frame(sample_data(scaled_core_genus))

metadata<-metadata[,c("SampleTimeCat",  "AgeCat", "Season", "Storage",  "Seq_run", "Seq_depth_nospikein")]

vectors_df<-ord_wunifrac_mds$vectors

# Model variables contributing to axis 3 and 4

fit <- envfit(vectors_df, metadata, perm = 999, na.rm = TRUE,  choices=c(3,4))# axes 3 and 4
df<-data.frame(vectors_df)
df$Dom_taxa<-sample_data(scaled_core_genus)$DomGenus_plot
df$TotalAbundance<-sample_data(scaled_core_genus)$TotalAbundance

#Get the vectors for bioenv.fit
df_biofit<-vegan::scores(fit,display=c("factors"))

## scaling factor for arrows to fill 80% of plot
mul <- ordiArrowMul(df_biofit, fill = 0.3)
df_biofit <- df_biofit * mul                    # Scale the biplot scores
df_biofit<-as.data.frame(df_biofit)

# remove methodological variables from biofit dataframe, 
#as I am just interested in biological variables after controling for methods
df_biofit<- df_biofit[c("AgeCatMIDDLEAGED", "AgeCatOLD", "AgeCatYOUNG", "SeasonDRY", "SeasonWET"),]
row.names(df_biofit)<-c("Adult", "Old", "Young", "Dry season", "Wet season" )

### spider plot
cent <- aggregate(cbind(Axis.3, Axis.4) ~ Dom_taxa, data = df, FUN = mean)

segs <- merge(df, setNames(cent, c('Dom_taxa','oNMDS1','oNMDS2')),
  by = 'Dom_taxa', sort = FALSE)

df2<-df
df_biofit2<-df_biofit
cent2<-cent
segs2<-segs


# order taxa so that they at least slightly match later figures!
df2$Dom_taxa<-factor(df2$Dom_taxa,levels=top_genera)
segs2$Dom_taxa<-factor(segs2$Dom_taxa,levels=top_genera)


P2<-ggplot(df2, aes(x = Axis.3, y = Axis.4)) +
  theme_light(base_size = 14)+
  geom_point(aes(fill = Dom_taxa), size = 3, pch = 21, alpha = 0.8, col = "white") + # spiders
   scale_fill_manual(values = mypaly[c(-16)])+
  scale_color_manual(values = mypaly[c(-16)])+
   annotate("text", x = 0, y = 0, label = "X", size = 6)+
 geom_point(data = cent2, size = 7, aes(fill = Dom_taxa), col = "black", pch = 21,  stroke = 1.5) +  
  geom_segment(data=df_biofit2, aes(x = 0, y = 0, xend = Axis.3*0.5, yend = Axis.4*0.5),
    arrow = arrow(length = unit(0.5, "cm"), type = "closed"), lwd = 0.5,  col = "black")+
  geom_label(data=as.data.frame(df_biofit2), #position=position_jitter(width=0.01, height = 0.01),#labels
    alpha = 0.9, col = "black", size = 4, fill = "yellow", 
    vjust = c(-0.2, 1, 1, -0.5, 1.5),
    hjust = c(1,0,0,0.5,0.7),
    aes(Axis.3*0.5,Axis.4*0.5, label = rownames(df_biofit2)))+
  xlab("Axis 3 (11.7%)")+
  ylab("Axis 4 (7.1%)")+
  guides(fill = guide_legend(override.aes = list(size = 7), ncol = 5))+
  theme(legend.title = element_blank())+
  theme(plot.margin = margin(0.2,0.2,1,0.2, "cm"))+
  annotate(geom="text", x=-0.3, y=-0.19, label="Age: R^2 = 0.043, p < 0.001 \nSeason: R^2 = 0.007, p = 0.003",
             color="black", size = 5, hjust = 0)




Fig3_1<-ggarrange(P1,P2, ncol = 2, nrow = 1, common.legend = T, legend = "bottom", labels = c("b)","c)"),   
                  font.label = list(size = 20))


ggarrange(Fig3_2,Fig3_1, ncol = 1, common.legend = T, heights = c(1,1.3))+
  theme(plot.margin = margin(0.1,0.1,0.1,1, "cm")) 

#ggsave("FIGURES/Figure3.pdf", width = 11, height = 12.5, units = "in")

```

\newpage

# Differential abundance analysis for all genera with \> 15% prevalence

#### Differential abundance analysis for morning/afternoon

```{r}

scaled_core_genus_0.2_original<- microbiome::core(scaled_core_genus, detection = 0, prevalence =0.16) 

scaled_core_genus_0.2<-scaled_core_genus_0.2_original

sample_data(scaled_core_genus_0.2)<-sample_data(scaled_core_genus_0.2)[,c("feature.id", "SampleTimeCat")]

table(sample_data(scaled_core_genus)$SampleTimeCat)


genus0.2_df<-speedyseq::psmelt(scaled_core_genus_0.2)

head(genus0.2_df)


## non-parametric

models_by_genus <- genus0.2_df %>%
  group_by(OTU) %>%
  do(model1 = tidy(wilcox.test(log10(Abundance+1)~ SampleTimeCat, data = ., exact=FALSE, conf.int=TRUE, conf.level=0.95), conf.int=TRUE)) %>%   
  gather(model_name, model, -OTU) %>%                        
  unnest()

models_by_genus$p.adj<- p.adjust(models_by_genus$p.value, method = p.adjust.methods, n = length(models_by_genus$p.value))

models_by_genus$Sig<-models_by_genus$p.adj<0.05

models_by_genus_time<-models_by_genus


P_time<-ggplot(models_by_genus_time, aes(y = reorder(OTU, -estimate), x = estimate))+
    geom_vline(xintercept = 0, col = "black")+
  geom_point(aes(col = Sig), size = 3)+  
   geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height = 0, col = Sig), size = 1.5)+
  theme_bw(base_size = 16)+
  labs(color = "Significant")+
  theme(axis.title = element_blank())+
 theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
   theme(plot.margin=unit(c(2,0.2,0.2,0.2),"cm"))+
  xlim(c(min(models_by_genus_time$conf.low),1.3))+
  ggtitle("a) Morning <--> Afternoon")+
 theme(plot.title = element_text(hjust = 1, size = 14, face = "bold"))

```

#### Differential abundance analysis for dry/wet season

```{r}

scaled_core_genus_0.2<-scaled_core_genus_0.2_original

sample_data(scaled_core_genus_0.2)<-sample_data(scaled_core_genus_0.2)[,c("feature.id", "Season")]
table(sample_data(scaled_core_genus)$Season)

genus0.2_df<-speedyseq::psmelt(scaled_core_genus_0.2)

head(genus0.2_df)

genus0.2_df$Season<-factor(genus0.2_df$Season, levels = c("WET SUMMER", "DRY WINTER"))



## non-parametric

models_by_genus <- genus0.2_df %>%
  group_by(OTU) %>%
  do(model1 = tidy(wilcox.test(log10(Abundance+1)~ Season, data = . , exact=TRUE, conf.int=TRUE, conf.level=0.95), conf.int=TRUE)) %>%   
  gather(model_name, model, -OTU) %>%                        
  unnest()

models_by_genus$p.adj<- p.adjust(models_by_genus$p.value, method = p.adjust.methods, n = length(models_by_genus$p.value))

models_by_genus$Sig<-models_by_genus$p.adj<0.05


models_by_genus_season<-models_by_genus

# positive values higher in wet season

P_season<-ggplot(models_by_genus_season, aes(y = reorder(OTU, -estimate), x = estimate))+
    geom_vline(xintercept = 0, col = "black")+
  geom_point(aes(col = Sig), size = 3)+  
 geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height = 0, col = Sig), size = 1.5)+
  theme_bw(base_size = 16)+
  theme(axis.title = element_blank())+
 theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
   theme(plot.margin=unit(c(2,0.2,0.2,0.2),"cm"))+
  xlim(c(min(models_by_genus_time$conf.low), max(models_by_genus_time$conf.high)))+
  ggtitle("b) Dry <--> Wet")+
 theme(plot.title = element_text(hjust = 1.2, size = 14, face = "bold"))

```

#### Differential abundance analysis for young/adult

```{r}

scaled_core_genus_0.2<-scaled_core_genus_0.2_original

## make a comparison group 2-4 years old (to reduce unequal sample sizes)

sample_data(scaled_core_genus_0.2)$AgeCat <- as.factor(ifelse(
  sample_data(scaled_core_genus_0.2)$AgeAtSamplingYears >2 & sample_data(scaled_core_genus_0.2)$AgeAtSamplingYears <4, 
  "ADULT",  as.character(sample_data(scaled_core_genus_0.2)$AgeCat)))


table(sample_data(scaled_core_genus_0.2)$AgeCat)

###

sample_data(scaled_core_genus_0.2)<-sample_data(scaled_core_genus_0.2)[,c("feature.id", "AgeCat")]

genus0.2_df<-speedyseq::psmelt(scaled_core_genus_0.2)

head(genus0.2_df)

genus0.2_df<-subset(genus0.2_df, AgeCat != "MIDDLEAGED" & AgeCat != "OLD" )

genus0.2_df$AgeCat<-factor(genus0.2_df$AgeCat)


## non-parametric

models_by_genus <- genus0.2_df %>%
  group_by(OTU) %>%
  do(model1 = tidy(wilcox.test(log10(Abundance+1)~ AgeCat, data = ., exact=FALSE, conf.int=TRUE, conf.level=0.95 ), conf.int=TRUE)) %>%   
  gather(model_name, model, -OTU) %>%                        
  unnest()

models_by_genus$p.adj<- p.adjust(models_by_genus$p.value, method = p.adjust.methods, n = length(models_by_genus$p.value))

models_by_genus$Sig<-models_by_genus$p.adj<0.05


models_by_genus_young<-models_by_genus


## if positive, more prevalent in adults

P_young<-ggplot(models_by_genus_young, aes(y = reorder(OTU, -estimate), x = estimate))+
  geom_vline(xintercept = 0, col = "black")+
  geom_point(aes(col = Sig), size = 3)+  
  geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height = 0, col = Sig), size = 1.5)+
  theme_bw(base_size = 16)+
  theme(axis.title = element_blank())+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(plot.margin=unit(c(2,0.2,0.2,0.2),"cm"))+
  xlim(c(min(models_by_genus_time$conf.low), 1.3))+
    ggtitle("c) Young <--> Adult")+
  theme(plot.title = element_text(hjust = 0.9, size = 14, face = "bold"))

```

#### Differential abundance analysis for adult/old

```{r fig.height=22, fig.width=20}

scaled_core_genus_0.2<-scaled_core_genus_0.2_original

## make a comparison group 2-4 years old (to reduce unequal sample sizes

sample_data(scaled_core_genus_0.2)$AgeCat <- as.factor(ifelse(
  sample_data(scaled_core_genus_0.2)$AgeAtSamplingYears >2 & sample_data(scaled_core_genus_0.2)$AgeAtSamplingYears <4, 
  "ADULT",  as.character(sample_data(scaled_core_genus_0.2)$AgeCat)))

table(sample_data(scaled_core_genus_0.2)$AgeCat)

###

sample_data(scaled_core_genus_0.2)<-sample_data(scaled_core_genus_0.2)[,c("feature.id", "AgeCat")]

genus0.2_df<-speedyseq::psmelt(scaled_core_genus_0.2)

genus0.2_df<-subset(genus0.2_df, AgeCat != "MIDDLEAGED" & AgeCat != "YOUNG" )

genus0.2_df$AgeCat<-factor(genus0.2_df$AgeCat, levels = c("OLD","ADULT"))

unique(genus0.2_df$AgeCat)


## non-parametric

models_by_genus <- genus0.2_df %>%
  group_by(OTU) %>%
  do(model1 = tidy(wilcox.test(log10(Abundance+1)~ AgeCat, data = ., exact=TRUE, conf.int=TRUE, conf.level=0.95 ), conf.int=TRUE)) %>%   
  gather(model_name, model, -OTU) %>%                        
  unnest()

models_by_genus$p.adj<- p.adjust(models_by_genus$p.value, method = p.adjust.methods, n = length(models_by_genus$p.value))

models_by_genus$Sig<-models_by_genus$p.adj<0.05


models_by_genus_old<-models_by_genus

P_old<-ggplot(models_by_genus_old, aes(y = reorder(OTU, -estimate), x = estimate))+
  geom_vline(xintercept = 0, col = "black")+
  geom_point(aes(col = Sig), size = 3)+  
  geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height = 0, col = Sig), size = 1.5)+
  theme_bw(base_size = 16)+
  theme(axis.title = element_blank())+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(plot.margin=unit(c(2,0.2,0.2,0.2),"cm"))+
  xlim(c(min(models_by_genus_time$conf.low), max(models_by_genus_old$conf.high)))+
  ggtitle("d) Adult <--> Old")+
  theme(plot.title = element_text(hjust = 0.7, size = 14, face = "bold"))


ggarrange(P_time, P_season, P_young, P_old, ncol = 4, common.legend = T, legend = "right")

#ggsave("FIGURES/wilcox_test_estimates.pdf", width = 20, height = 22, units = "in")

```

# Selecting 16 focus genera

-   We applied the 16 focus genera already in Fig. 3, but now we outline how we chose those 16 taxa
-   These are chosen for further analysis on their temporal dynamics
-   These genera are selected because they are the most abundant, but also because they show significant and notable dynamics
-   Because we want to limit the focus taxa to \<16 (for clarity), a small number of quite prevalent/abundant genera are excluded from further analysis, eg Geodermatophilus (likely a soil bacterium ingested during foraging), Fusobacterium, Collinsella, and Catenisphaera. These taxa don't show very strong nor robust temporal dynamics and therefore we do not focus on them, although they are presented in the differential abundance analysis.

```{r, echo = F}
### functions used in core analysis

prevalence <- function(physeq, add_tax = TRUE){
  
  ## Check if taxa are rows
  trows <- taxa_are_rows(physeq)
  
  ## Extract OTU table
  otutab <- as.data.frame(otu_table(physeq))
  
  ## Transpose OTU table (species should be arranged by rows)
  if(trows == FALSE){
    otutab <- t(otutab)
  }
  
  ## Estimate prevalence (number of samples with OTU present)
  prevdf <- apply(X = otutab,
                  #MARGIN = ifelse(trows, yes = 1, no = 2),  # for a non-transposed data
                  MARGIN = 1,
                  FUN = function(x){sum(x > 0)})
  
  ## Add total and average read counts per OTU
  prevdf <- data.frame(Prevalence = prevdf,
                       TotalAbundance = taxa_sums(physeq),
                       MeanAbundance = rowMeans(otutab),
                       MedianAbundance = apply(otutab, 1, median))
  
  ## Add taxonomy table
  if(add_tax == TRUE && !is.null(tax_table(physeq, errorIfNULL = F))){
    prevdf <- cbind(prevdf, tax_table(physeq))
  }
  return(prevdf)
}


```

#### Make table of genus characteristics

-   Calculate prevalence, abundance, and number of samples where that genus is dominant

```{r}

# calculate prevalence per genus

genus_info_df<-data.frame(microbiome::prevalence(scaled_core_genus))
names(genus_info_df)<-"Prevalence"
genus_info_df$Genus<-row.names(genus_info_df)
genus_info_df<-genus_info_df[,c(2,1)]

### add mean abundance as column

genus_info_df$Abundance<-taxa_sums(scaled_core_genus)
genus_info_df$RelAbundance<-taxa_sums(scaled_core_genus)/sum(taxa_sums(scaled_core_genus))


####### add dominant taxa as column
dominant_phylo<-scaled_core_genus
sample_data(dominant_phylo)<-sample_data(dominant_phylo)[,"feature.id"]

# melt data
psmelt_df<-speedyseq::psmelt(dominant_phylo)
#head(psmelt_df)

dominant_taxa_df<-psmelt_df %>% 
  group_by(feature.id) %>% 
  arrange(feature.id, -Abundance) %>% 
  dplyr::mutate(rank=rank(-Abundance))

dominant_taxa_df<-subset(dominant_taxa_df, rank ==1)
dominant_taxa_df$Genus<-factor(dominant_taxa_df$Genus)

# add to metadata 
sample_data(scaled_core)$DomGenus <-vlookup(sample_data(scaled_core)$feature.id, dominant_taxa_df, lookup_column = "Sample", result_column = "Genus")

dominant_genus_freq<-data.frame(table(dominant_taxa_df$OTU))
names(dominant_genus_freq)<-c("Genus", "Dominant")

genus_info_df$Dominant<-vlookup(genus_info_df$Genus, dominant_genus_freq, lookup_column = "Genus", result_column = "Dominant")

genus_info_df<-genus_info_df[order(-genus_info_df$RelAbundance),]
#head(genus_info_df, 10) # ordered by relative abundance
```

-   Keep only genera with over 60% prevalence

```{r}

top_genera_candidates<-subset(genus_info_df, Prevalence > 0.6)
top_genera_candidates<-top_genera_candidates[order(-top_genera_candidates$Abundance),]

DT::datatable(top_genera_candidates)

### based on differential abundance analysis, make a list of genera to  exclude

top_genera<-c( "Clostridium sensu stricto 1" ,
               "Blautia" ,
               "Unknown Lachnospiraceae",
               "Enterococcus"  ,
               "[Ruminococcus] torques group",
               "Bacteroides",  
               "Paeniclostridium",
               "Romboutsia" ,
               "Lachnoclostridium",  
               "Lactococcus"  ,
               "Peptococcus",
               "Helicobacter" ,
               "Unknown Bacillaceae",
               "Raoultibacter"  ,
               "Cellulomonas", 
               "Eubacterium"  ,
               "Other" )

top_genera

focus_taxa_df<-genus_info_df[top_genera[-17],]

focus_taxa_df$Prevalence<-round(focus_taxa_df$Prevalence, 3)
focus_taxa_df$RelAbundance<-round(focus_taxa_df$RelAbundance, 3)
row.names(focus_taxa_df)<-1:16
```

### Prevalence and abundance of 16 focus genera

```{r}

# info on 16 focus genera

DT::datatable(focus_taxa_df)

```

\newpage

# Figure 4: Temporal dynamics of 16 focus genera

-   Fit GAMMs to abundances of each focus genus
-   Repeat models but fit smooths to each storage method seperately to assess robustness

```{r merge-genus, message = FALSE, warning = FALSE}

core_genus <- prune_taxa(top_genera, scaled_core_genus)

paste("Core taxa make up",round(sum(taxa_sums(core_genus))/sum(taxa_sums(meerkat_scaled)),2), "% of reads")

```

-   Fit GAMMs to 16 focus genera and extract estimates using loop.

```{r genus-gamms, message = FALSE, warning = FALSE}
## make vector of genus names
taxonomy<-data.frame(core_genus@tax_table@.Data)
taxa_names(core_genus)<-taxonomy$Genus

## make vector of variables to include in model

variables_all<-c("HoursAfterSunrise", "month", "AgeAtSamplingYears", 
  "IndividID", "GroupAtSampling",
  "Storage", "Seq_depth_nospikein", "Seq_run", "HoursTilFrozen")


### loop to model each genus seperatly and extract estimates
taxanames<-as.character(taxa_names(core_genus))

```

```{r}
base_genus_confint<-list()

for (i in 1:length(taxanames)){
  tryCatch({ #catch errors
    #print(i)
  
    taxa1<-prune_taxa(taxanames[i], core_genus)
   # taxa1<-prune_taxa("Clostridium sensu stricto 1"   , core_genus)
    sample_data(taxa1)$taxa_i_abundance<-sample_sums(taxa1)
    metadata<-data.frame(sample_data(taxa1))
    metadata<-metadata[,c(variables_all, "taxa_i_abundance", "Year")]
    
    #scale variables
    metadata$Seq_depth_nospikein<-as.numeric(scale(metadata$Seq_depth_nospikein))
    
    # fit gamm    
    m_taxa <- mgcv::gamm(log10(taxa_i_abundance+1)~
        s(HoursAfterSunrise, bs = "cr")+
        s(month, bs = "cc")+
        s(AgeAtSamplingYears, bs = "cr")+
        s(Seq_depth_nospikein, bs = "cr")+
       Storage+
        Seq_run+
        HoursTilFrozen+
        s(IndividID, bs = "re")+
        s(GroupAtSampling, bs = "re"),
      correlation = corARMA(form = ~ 1|Year, p = 3),
      data=metadata,
      family = gaussian)
    
    #draw(m_taxa$gam, select = c(1:3), ncol = 3, scales = "fixed")+geom_hline(yintercept = 0)
    
    confint_df_time<-confint(m_taxa$gam, parm = "s(HoursAfterSunrise)",  type = "confidence")
    confint_df_season<-confint(m_taxa$gam, parm = "s(month)",  type = "confidence")
    confint_df_age<-confint(m_taxa$gam, parm = "s(AgeAtSamplingYears)",  type = "confidence")
    
    confint_df_time$Taxa<-taxanames[i]
    confint_df_time$Variable<-"HoursAfterSunrise"
    names(confint_df_time)[3]<-"Time"
    confint_df_time$p.val<-data.frame(summary(m_taxa$gam)$s.table)$p.value[1]
    confint_df_season$Taxa<-taxanames[i]
    confint_df_season$Variable<-"month"
    names(confint_df_season)[3]<-"Time"
   confint_df_season$p.val<-data.frame(summary(m_taxa$gam)$s.table)$p.value[2]

    confint_df_age$Taxa<-taxanames[i]
    confint_df_age$Variable<-"AgeAtSamplingYears"
    names(confint_df_age)[3]<-"Time"
    confint_df_age$p.val<-data.frame(summary(m_taxa$gam)$s.table)$p.value[3]

    
    confint_df<-rbind(confint_df_time, confint_df_season, confint_df_age)
    confint_df$R2<-summary(m_taxa$gam)$r.sq
    
    # convert into real estimate by adding the intercapt (= mean abundance)
    
    confint_df$est_real<-confint_df$est+summary(m_taxa$gam)$p.coeff[1] 
    confint_df$lower_real<- confint_df$lower+summary(m_taxa$gam)$p.coeff[1] 
    confint_df$upper_real<- confint_df$upper+summary(m_taxa$gam)$p.coeff[1] 
    
    confint_df
    base_genus_confint[[i]]<-confint_df
    
  }, error=function(e){})
}

names(base_genus_confint)<-taxanames

base_genus_confint[1]


```

```{r}

confint_summary<-do.call(rbind, base_genus_confint)
confint_summary$Variable<-factor(confint_summary$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))
names(confint_summary)[which(colnames(confint_summary)=="p.val")]<-"p.val.all"
confint_summary$Significant<-confint_summary$p.val.all <0.05

## add abundance for plot

abundance_df<-data.frame(taxa_sums(core_genus))
names(abundance_df)<-"Abundance"
abundance_df$Genus<-row.names(abundance_df)
confint_summary$Abundance<-as.numeric(vlookup(confint_summary$Taxa, abundance_df, lookup_column = "Genus", result_column = "Abundance"))
confint_summary<-data.frame(confint_summary)
confint_summary<-dplyr::arrange(confint_summary, desc(Abundance))
confint_summary$Taxa<-factor(confint_summary$Taxa, levels = rev(top_genera[-17]))

```

-   Repeat models but by storage to assess robustness
-   Any genus-level dynamics that are not robust across storage types are excluded

```{r}
taxanames<-as.character(taxa_names(core_genus))
base_genus_confint_storage<-list()

for (i in 1:length(taxanames)){
  tryCatch({ #catch errors
    #print(i)
    
    taxa1<-prune_taxa(taxanames[i], core_genus)
    # taxa1<-prune_taxa("Clostridium sensu stricto 1"   , core_genus)
    sample_data(taxa1)$taxa_i_abundance<-sample_sums(taxa1)
    metadata<-data.frame(sample_data(taxa1))
    metadata<-metadata[,c(variables_all, "taxa_i_abundance", "Year")]
    
    #scale variables
    metadata$Seq_depth_nospikein<-as.numeric(scale(metadata$Seq_depth_nospikein))
    
    # fit gamm    
    m_taxa <- mgcv::gamm(log10(taxa_i_abundance+1)~
                           s(HoursAfterSunrise, bs = "cr", by = Storage)+
                           s(month, bs = "cc",  by = Storage)+
                           s(AgeAtSamplingYears, bs = "cr",  by = Storage)+
                           s(Seq_depth_nospikein, bs = "cr",  by = Storage)+
                           Seq_run+
                           HoursTilFrozen+
                           s(IndividID, bs = "re")+
                           s(GroupAtSampling, bs = "re"),
                         correlation = corARMA(form = ~ 1|Year, p = 3),
                         data=metadata,
                         family = gaussian)
    
    summary(m_taxa$gam)
    
    #draw(m_taxa$gam, select = c(1:3), ncol = 3, scales = "fixed")+geom_hline(yintercept = 0)
    
    confint_df_time<-confint(m_taxa$gam, parm = "s(HoursAfterSunrise)",  type = "confidence")
    confint_df_season<-confint(m_taxa$gam, parm = "s(month)",  type = "confidence")
    confint_df_age<-confint(m_taxa$gam, parm = "s(AgeAtSamplingYears)",  type = "confidence")
    
    confint_df_time$Taxa<-taxanames[i]
    confint_df_time$Variable<-"HoursAfterSunrise"
    names(confint_df_time)[3]<-"Time"
    ### add p value
    confint_df_time$p.val<-NA
    confint_df_time$p.val[1:200]<-data.frame(summary(m_taxa$gam)$s.table)$p.value[1]
    confint_df_time$p.val[201:400]<-data.frame(summary(m_taxa$gam)$s.table)$p.value[2]
   
    
    confint_df_season$Taxa<-taxanames[i]
    confint_df_season$Variable<-"month"
    names(confint_df_season)[3]<-"Time"
    ### add p value
    confint_df_season$p.val<-NA
    confint_df_season$p.val[1:200]<-data.frame(summary(m_taxa$gam)$s.table)$p.value[3]
    confint_df_season$p.val[201:400]<-data.frame(summary(m_taxa$gam)$s.table)$p.value[4]
    
    confint_df_age$Taxa<-taxanames[i]
    confint_df_age$Variable<-"AgeAtSamplingYears"
    names(confint_df_age)[3]<-"Time"
    ### add p value
    confint_df_age$p.val<-NA
    confint_df_age$p.val[1:200]<-data.frame(summary(m_taxa$gam)$s.table)$p.value[5]
    confint_df_age$p.val[201:400]<-data.frame(summary(m_taxa$gam)$s.table)$p.value[6]
    
    
    confint_df<-rbind(confint_df_time, confint_df_season, confint_df_age)
    confint_df$R2<-summary(m_taxa$gam)$r.sq
    

   base_genus_confint_storage[[i]]<-confint_df
    
  }, error=function(e){})
}

names(base_genus_confint_storage)<-taxanames

#saveRDS(base_genus_confint_storage, "R_OUTPUT/base_genus_confint_storage.RDS")

```


```{r}

confint_genus_storage<-do.call(rbind, base_genus_confint_storage)
confint_genus_storage$Variable<-factor(confint_genus_storage$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))
confint_genus_storage$Taxa<-factor(confint_genus_storage$Taxa, levels = top_genera[-17])
storage_model_summary<-distinct(confint_genus_storage, Storage, Taxa, Variable, .keep_all = T )

storage_model_frozen<-subset(storage_model_summary, Storage=="FROZEN")
storage_model_frozen<-storage_model_frozen[,c("Variable","Taxa","p.val")]
dim(confint_summary)
confint_summary<-merge(confint_summary, storage_model_frozen, by = c("Variable", "Taxa"))
names(confint_summary)[which(colnames(confint_summary)=="p.val")]<-"p.val_frozen"

storage_model_freezedried<-subset(storage_model_summary, Storage=="FREEZEDRIED")
storage_model_freezedried<-storage_model_freezedried[,c("Variable","Taxa","p.val")]
dim(confint_summary)
confint_summary<-merge(confint_summary, storage_model_freezedried, by = c("Variable", "Taxa"))
names(confint_summary)[which(colnames(confint_summary)=="p.val")]<-"p.val_freezedried"

# add column for whether trend is robust
confint_summary$Sig_robust<-ifelse(confint_summary$p.val_frozen < 0.05 & confint_summary$p.val_freezedried < 0.05, TRUE, FALSE)
  
```

#### Figure 4e-g: Back-transformed model estimates for genus-level abundance across temporal scales

-   Add bacterial load to figure

```{r figure4efg, message = FALSE, warning = FALSE,  fig.width = 9, fig.height = 5}

### rugdata
metadata<-data.frame(sample_data(core_genus))
rugdata<- data.frame(cbind(metadata$HoursAfterSunrise, metadata$month, metadata$AgeAtSamplingYears))
names(rugdata)<-c("HoursAfterSunrise", "month", "AgeAtSamplingYears")
rugdata<-gather(rugdata, Variable, Time, HoursAfterSunrise:AgeAtSamplingYears, factor_key=TRUE)
rugdata$est = 0
rugdata$Variable<-factor(rugdata$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))

##########

confint_summary_fig4b<-confint_summary[,c("Variable", "Taxa", "Time", "est_real","Sig_robust", "Abundance")]
confint_df_load_fig4b<-confint_df_load[,c("Variable", "Time", "est")]
confint_df_load_fig4b$Taxa <- "Bacterial load"
confint_df_load_fig4b$Sig_robust <- TRUE
confint_df_load_fig4b$Abundance <- 43124192 # same abundance of clostridium

confint_df_load_fig4b<-confint_df_load_fig4b[,c(1,4,2,3,5,6)]
names(confint_df_load_fig4b)<-c("Variable", "Taxa", "Time", "est_real","Sig_robust", "Abundance")
confint_df_load_fig4b$est_real<-log10(confint_df_load_fig4b$est_real)

confint_summary_fig4b<-rbind(confint_summary_fig4b,confint_df_load_fig4b )

confint_summary_fig4b$Taxa<-factor(confint_summary_fig4b$Taxa, levels = rev(c(top_genera[-17], "Bacterial load")))


Fig4def<-ggplot(confint_summary_fig4b, aes(x = Time, y = 10^est_real, group = Taxa))+
  theme_light(base_size = 14)+
 geom_line( aes(alpha = Sig_robust, col = Taxa, size = Abundance))+
  #geom_ribbon( aes(ymin = 10^lower_real, ymax = 10^upper_real, fill = Taxa, alpha = Sig_robust, size = 1))+
  scale_size(range = c(1, 1.5))+
  facet_wrap(~Variable, scales = "free_x")+
  geom_rug(data = rugdata, aes(x = Time, y = NULL), inherit.aes = F, size = 0.1)+
  scale_color_manual(values = rev(c(mypaly[-17], "darkgrey")),  labels = c(top_genera[-17], "Bacterial load") )+
  ylab("Estimated abundance")+
  scale_x_continuous(breaks = int_breaks)+
  theme(panel.grid = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))+
theme(panel.spacing = unit(1.2, "lines"))+
 theme(plot.margin=unit(c(0.2,0.2,0.2,0.6),"cm"))+
  theme(plot.title = element_text(size=10))+
  guides(size=FALSE, abundance = FALSE, alpha = FALSE)+
  guides(color = guide_legend(override.aes = list(size = 5)))+
  #theme(legend.position = "none")+
  theme(strip.background = element_blank(),
        strip.text.x = element_blank())+
  theme(axis.title.x = element_blank())+
  scale_alpha_manual(values = c(0,0.8))+
  scale_y_log10(limits = c(0.7,1e6), breaks = trans_breaks("log10", function(x) 10^x),
              labels = trans_format("log10", math_format(10^.x)))+
    guides(color = guide_legend(override.aes = list( size = 4, col = mypaly, 
                                                    labels = c(top_genera))))+
   labs(color='Genus') 

```

#### Figure 4e-g: Untransformed model estimates for genus-level abundance across temporal scales

-   Remove bacterial load

```{r fig.height=5, fig.width=10}

confint_sumamry_fig4a<-confint_summary[,c("Variable", "Taxa", "Time", "est","Sig_robust", "Abundance")]
confint_df_load_fig4a<-confint_df_load[,c("Variable", "Time", "est")]
confint_df_load_fig4a$Taxa <- "Bacterial load"
confint_df_load_fig4a$Sig_robust <- FALSE
confint_df_load_fig4a$Abundance <- 43124192 # same abundance of clostridium
str(confint_df_load_fig4a)
str(confint_sumamry_fig4a)
confint_df_load_fig4a<-confint_df_load_fig4a[,c(1,4,2,3,5,6)]
names(confint_df_load_fig4a)<-c("Variable", "Taxa", "Time", "est","Sig_robust", "Abundance")
confint_df_load_fig4a$est<-0

confint_sumamry_fig4a<-rbind(confint_sumamry_fig4a,confint_df_load_fig4a )


confint_sumamry_fig4a$Taxa<-factor(confint_sumamry_fig4a$Taxa, levels = rev(c(top_genera[-17], "Bacterial load")))

confint_sumamry_fig4a$est<-ifelse(confint_sumamry_fig4a$Taxa=="Raoultibacter", confint_sumamry_fig4a$est+0.03, confint_sumamry_fig4a$est)

Fig4abc<-ggplot(confint_sumamry_fig4a, aes(x = Time, y = est, group = Taxa))+
  theme_light(base_size = 14)+
 geom_line( aes(alpha = Sig_robust, col = Taxa, size = Abundance))+
 #geom_line( aes( col = Taxa))+
  scale_size(range = c(1, 1.5))+
  facet_wrap(~Variable, scales = "free_x")+
  geom_rug(data = rugdata, aes(x = Time, y = NULL), inherit.aes = F, size = 0.1)+
  scale_color_manual(values = rev(c(mypaly[-17], "darkgrey")),  labels = c(top_genera[-17], "Bacterial load") )+
  ylab("Model estimate (log10)")+
  scale_x_continuous(breaks = int_breaks)+
  theme(panel.grid = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))+
theme(panel.spacing = unit(1.2, "lines"))+
 theme(plot.margin=unit(c(0.2,0.2,0.2,0.4),"cm"))+
  theme(plot.title = element_text(size=10))+
  guides(size=FALSE, alpha = FALSE, abundance = FALSE)+
  theme(strip.background = element_blank(),
        strip.text.x = element_blank())+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank())+
  geom_hline(yintercept = 0, linetype = "dashed")+
  scale_alpha_manual(values = c(0,0.8))+
    guides(color = guide_legend(override.aes = list( size = 4, col = mypaly, 
                                                    labels = c(top_genera))))+
   labs(color='Genus') 

ggarrange(Fig4abc, Fig4def, ncol = 1, heights = c(1,1), legend = "right", common.legend = T)

#ggsave("FIGURES/figure4.pdf",width = 10, height = 5, units = "in")


```

### Supplementary figures: Dynamics by storage type

```{r fig.height=8, fig.width=10}

confint_genus_storage<-do.call(rbind, base_genus_confint_storage)
confint_genus_storage$Variable<-factor(confint_genus_storage$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))

confint_genus_storage$Taxa<-factor(confint_genus_storage$Taxa, levels = top_genera[-17])

confint_genus_storage$Sig <-confint_genus_storage$p.val < 0.05



### rugdata
metadata<-data.frame(sample_data(core_genus))
rugdata<- data.frame(cbind(metadata$HoursAfterSunrise, metadata$month, metadata$AgeAtSamplingYears))
names(rugdata)<-c("HoursAfterSunrise", "month", "AgeAtSamplingYears")
rugdata<-gather(rugdata, Variable, Time, HoursAfterSunrise:AgeAtSamplingYears, factor_key=TRUE)
rugdata$est = 0
rugdata$Variable<-factor(rugdata$Variable, levels = c("HoursAfterSunrise", "month", "AgeAtSamplingYears"))

### Diurnal

diurnal_storage<-subset(confint_genus_storage, Variable == "HoursAfterSunrise")
head(diurnal_storage)


ggplot(diurnal_storage, aes(x = Time, y = est, group = Storage))+
  geom_line(aes(col = Storage, linetype = Sig) ,size = 1)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = Storage), linetype=2, alpha=0.5)+
  facet_wrap(~Taxa, scales = "free", ncol = 4)+
  theme_bw(base_size = 14)+
  scale_color_manual(values = c("deeppink", "deepskyblue3"))+
  scale_fill_manual(values = c( "deeppink", "deepskyblue3"))+
  ylab("Model estimate")+
  xlab("Hours after sunrise")+
  guides(linetype=guide_legend(title="Significant (p < 0.05)"))

#ggsave("FIGURES/FigS4_diurnal.pdf", width = 12, height = 8, units = "in")

#### seasonal

seasonal_storage<-subset(confint_genus_storage, Variable == "month")


ggplot(seasonal_storage, aes(x = Time, y = est, group = Storage))+
  geom_line(aes(col = Storage, linetype = Sig) ,size = 1)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = Storage), linetype=2, alpha=0.5)+
  facet_wrap(~Taxa, scales = "free", ncol = 4)+
  theme_bw(base_size = 14)+
  scale_color_manual(values = c("deeppink", "deepskyblue3"))+
  scale_fill_manual(values = c( "deeppink", "deepskyblue3"))+
  ylab("Model estimate")+
  xlab("Month")+
  guides(linetype=guide_legend(title="Significant (p < 0.05)"))+
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10, 11, 12 ))

#ggsave("FIGURES/FigS5_seasonal.pdf", width = 12, height = 8, units = "in")


##### age 


age_storage<-subset(confint_genus_storage, Variable == "AgeAtSamplingYears")


ggplot(age_storage, aes(x = Time, y = est, group = Storage))+
  geom_line(aes(col = Storage, linetype = Sig) ,size = 1)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = Storage), linetype=2, alpha=0.5)+
  facet_wrap(~Taxa, scales = "free", ncol = 4)+
  theme_bw(base_size = 14)+
  scale_color_manual(values = c("deeppink", "deepskyblue3"))+
  scale_fill_manual(values = c( "deeppink", "deepskyblue3"))+
  ylab("Model estimate")+
  xlab("Age at sampling (years)")+
  guides(linetype=guide_legend(title="Significant (p < 0.05)"))+
  scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10 ))

#ggsave("FIGURES/FigS6_age.pdf", width = 12, height = 8, units = "in")


```

#### Sample distributions by storage

```{r  fig.height=3, fig.width=5, eval = F, echo = F}

# Time of day
  
ggplot(sample_data(meerkat_scaled), aes(x = HoursAfterSunrise))+
  geom_histogram(aes(fill = Storage),position = "dodge", alpha = 1, binwidth = 0.5, col = "black")+
  scale_fill_manual(values = c( "deeppink", "deepskyblue3"))+
  scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10 ))+
  facet_wrap(~Storage, ncol = 1)+
  theme_bw(base_size = 16)+
  theme(legend.position = "none")+
  ylab("Frequency")+
  xlab ("Hours after sunrise")+
  scale_x_continuous(breaks = c(0,2,4,6,8,10, 12 ,14))


ggplot(sample_data(meerkat_scaled), aes(x = month))+
  geom_histogram(aes(fill = Storage),position = "dodge", alpha = 1, binwidth = 0.5, col = "black")+
  scale_fill_manual(values = c( "deeppink", "deepskyblue3"))+
  scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10 ))+
  facet_wrap(~Storage, ncol = 1)+
  theme_bw(base_size = 16)+
  theme(legend.position = "none")+
  ylab("Frequency")+
  xlab ("Month")


ggplot(sample_data(meerkat_scaled), aes(x = AgeAtSamplingYears))+
  geom_histogram(aes(fill = Storage),position = "dodge", alpha = 1, binwidth = 0.5, col = "black")+
  scale_fill_manual(values = c( "deeppink", "deepskyblue3"))+
  scale_x_continuous(breaks = c(0,1,2,3,4,5,6,7,8,9,10 ))+
  facet_wrap(~Storage, ncol = 1)+
  theme_bw(base_size = 16)+
  theme(legend.position = "none")+
  ylab("Frequency")+
  xlab ("Age (years)")


```

# Figure 5: Biological predictors of diversity and genus-level abundance

-   Fit GAMM to each core genus and diversity measures
-   This is a very long code chunk. Models are run accounting for storage AND when split by storage.
-   Models are also run excluding various terms to partition R2 into temporal, biological and methodological variance.
-   extract_fixed() function needs to be loaded from <https://rdrr.io/github/m-clark/gammit/src/R/extract_fixed.R>.

```{r function, echo = FALSE, message = FALSE, warning = FALSE}
# function to extract estimates

extract_fixed <- function(
  model,
  ci_level = .95,
  digits = 3,
  ...
) {
  fe <- data.frame(summary(model, ...)$p.table)
  colnames(fe) =  c('value', 'se', 't', 'p')
  
  # no confint.gam
  if (ci_level > 0) {
    
    lower <- (1 - ci_level)/2
    upper <- 1 - lower
    nu <- model$df.residual
    mult <- stats::qt(upper, nu)
    
    ci <- data.frame(
      lower = fe$value - mult * fe$se,
      upper = fe$value + mult * fe$se
    )
    
    colnames(ci) <- paste0(c('lower_', 'upper_'), c(lower, upper) * 100)
    
    fe <- data.frame(fe, ci)
  }
  
  fe <- fe %>%
    dplyr::mutate_all(round, digits = digits) %>%
    dplyr::mutate(term = gsub(rownames(fe),
      pattern = '[\\(,\\)]',
      replacement = '')) %>%
    dplyr::select(term, dplyr::everything()) %>%
    dplyr::as_tibble()
  
  fe
  
}

```

### Figure 5a: Base models without biological predictors

This section fits base GAMMs (methodological and temporal variables only) and full GAMMs (including biological variables) to all microbial phenotypes, and extracts model summaries.

-   Make list of phenotypes (genus-level phenotypes and community-wide phenotypes)
-   Add these values as columns to metadata
-   Build 'base models' (models with just methodological and temporal predictors) for each phenotype
-   Extract summary statistics for temporal variables
-   This is kind of a repeat of earlier code, but we extract different information from the models (eg 95% CIs instead of predicted estimates).

```{r}

### add ordination axes
sample_data(core_genus)$MDS1_W<-data.frame(ord_wunifrac_mds$vectors)$Axis.1
sample_data(core_genus)$MDS1_UW<-data.frame(ord_unifrac_mds$vectors)$Axis.1

diversity_names<-c("TotalAbundance", "Observed", "Shannon", "MDS1_W", "MDS1_UW")
taxanames
phenotypes<-c(taxanames, diversity_names)

# make empty list to store results
  
base_model_summary_list<-list()


for (i in 1:length(phenotypes)){
  
  tryCatch({ #catch errors
    
    #print(i)
    
    metadata<-data.frame(sample_data(core_genus))
   # str(metadata)
    asvtable<-data.frame(t(otu_table(core_genus)))
    asvtable_log<-  log10(asvtable+1)
    dim(asvtable_log)
    names(asvtable_log)<-taxanames
    
    metadata_all<-merge(metadata, asvtable_log, by = 0)
    
   row.names(metadata_all)<-metadata_all$feature.id

    
    #scale variables
    metadata_all$Seq_depth_nospikein<-as.numeric(scale(metadata_all$Seq_depth_nospikein))
    metadata_all$HoursTilFrozen<-as.numeric(scale(metadata_all$HoursTilFrozen))
    metadata_all$TotalAbundance<-log10(metadata_all$TotalAbundance+1) #log10 transform load
    metadata_all$Observed<-scale(metadata_all$Observed) # scale observed
    metadata_all$Shannon<-scale(metadata_all$Shannon) # scale shannon

    
    # fit gamm
    m_taxa <- mgcv::gamm(metadata_all[, phenotypes[i]]~
        s(HoursAfterSunrise, bs = "cr")+
        s(month, bs = "cc")+
        s(AgeAtSamplingYears, bs = "cr")+
        s(Seq_depth_nospikein, bs = "cr")+
        Storage+
        Seq_run+
        HoursTilFrozen+
        s(IndividID, bs = "re")+
       s(GroupAtSampling, bs = "re"),
     correlation = corARMA(form = ~ 1|Year, p = 3),
      data=metadata_all,
      family = gaussian)
    
    summary<- summary(m_taxa$gam) 

    
    fixed_terms<-data.frame(extract_fixed(m_taxa$gam))[,c(2,5:7)] ## need to load the function here
    row.names(fixed_terms)<-extract_fixed(m_taxa$gam)$term
    names(fixed_terms)<-c("Effect", "P.val", "lower_CI", "upper_CI")
    
    smoothed_terms<-data.frame(summary$s.table)
    names(smoothed_terms)<-c("lower_CI", "upper_CI", "Effect", "P.val")
    smoothed_terms$lower_CI <-0
    smoothed_terms$upper_CI <-0
    
    results_df<-rbind(fixed_terms, smoothed_terms)
    results_df$R.sq<-summary$r.sq
    results_df$Effect_abs<-abs( results_df$Effect)
    results_df$Variable<-row.names(results_df)
    results_df$Taxa<-phenotypes[i]
    
    base_model_summary_list[[i]]<-results_df
    
  }, error=function(e){})
}

names(base_model_summary_list)<-phenotypes

#saveRDS(base_model_summary_list, "R_OUTPUT/base_models.RDS")
```

```{r}

base_model_summary<-do.call(rbind, base_model_summary_list)
base_model_summary$Significance<-base_model_summary$P.val < 0.01 # add column = whether significant
base_model_summary<-subset(base_model_summary, Variable !="Intercept") # remove intercept row

base_model_summary$Variable<- factor(base_model_summary$Variable, levels = c( # order variables
  "s(HoursAfterSunrise)", "s(month)" , "s(AgeAtSamplingYears)",   
  "s(IndividID)", "s(GroupAtSampling)",
  "StorageFROZEN","s(Seq_depth_nospikein)", "Seq_runRUN2", "Seq_runRUN3","Seq_runRUN4", "HoursTilFrozen"))

```

-   Split base models by storage type to assess how robust associations are

```{r}
# make PS object for frozen and freezedried samples

core_genus_frozen<-subset_samples(core_genus, Storage == "FROZEN")
core_genus_freezedried<-subset_samples(core_genus, Storage == "FREEZEDRIED")

# put in a list

ps_list<-list(core_genus_frozen, core_genus_freezedried)
names(ps_list)<-c("FROZEN", "FREEZEDRIED")

# make empty lists to store results

base_model_storage_i<-list()
base_model_storage_k<-list()

# loops loops loops

for (k in 1:length(ps_list)){
  
  core_genus_storage<- ps_list[[k]]
 # core_genus_storage<- ps_list[[2]]
  
  for (i in 1:length(phenotypes)){
  
  tryCatch({ #catch errors
    
    #print(i)
    
    # generate metadata dataframe
    metadata<-data.frame(sample_data(core_genus_storage))
    asvtable<-data.frame(t(otu_table(core_genus_storage))) # add genus-level abundances to dataframe
    asvtable_log<-  log10(asvtable+1)
    names(asvtable_log)<-taxanames
    metadata_all<-merge(metadata, asvtable_log, by = 0)
    row.names(metadata_all)<-metadata_all$feature.id

    #scale variables
    metadata_all$Seq_depth_nospikein<-as.numeric(scale(metadata_all$Seq_depth_nospikein))
    metadata_all$HoursTilFrozen<-as.numeric(scale(metadata_all$HoursTilFrozen))
    metadata_all$TotalAbundance<-log10(metadata_all$TotalAbundance+1) #log10 transform load
    metadata_all$Observed<-scale(metadata_all$Observed) # scale observed
    metadata_all$Shannon<-scale(metadata_all$Shannon) # scale shannon

    
    # fit gamm
    m_taxa <- mgcv::gamm(metadata_all[, phenotypes[i]]~
   # m_taxa <- mgcv::gamm(metadata_all[, phenotypes[13]]~
        s(HoursAfterSunrise, bs = "cr")+
        s(month, bs = "cc")+
        s(AgeAtSamplingYears, bs = "cr")+
        s(Seq_depth_nospikein, bs = "cr")+
        Seq_run+
        HoursTilFrozen+
        s(IndividID, bs = "re")+
       s(GroupAtSampling, bs = "re"),
       correlation = corARMA(form = ~ 1|Year, p = 3),
      data=metadata_all,
      family = gaussian)
    
    # extract summary statistics
    summary<- summary(m_taxa$gam) 
    fixed_terms<-data.frame(extract_fixed(m_taxa$gam))[,c(2,5:7)] ## need to load the function here
    row.names(fixed_terms)<-extract_fixed(m_taxa$gam)$term
    names(fixed_terms)<-c("Effect", "P.val", "lower_CI", "upper_CI")
    
    smoothed_terms<-data.frame(summary$s.table)
    names(smoothed_terms)<-c("lower_CI", "upper_CI", "Effect", "P.val")
    smoothed_terms$lower_CI <-0
    smoothed_terms$upper_CI <-0
    
    results_df<-rbind(fixed_terms, smoothed_terms)
    results_df$R.sq<-summary$r.sq
    results_df$Effect_abs<-abs( results_df$Effect)
    results_df$Variable<-row.names(results_df)
    results_df$Taxa<-phenotypes[i]
    results_df$Type <- as.character(metadata_all$Storage[1])
    
    base_model_storage_i[[i]]<-results_df
    
  }, error=function(e){})
    
  }
    
    base_model_storage_df<-do.call(rbind, base_model_storage_i)
    
    base_model_storage_k[[k]] <-base_model_storage_df
    
}

base_model_storage_k_df<-do.call(rbind, base_model_storage_k)


#saveRDS(base_model_storage_k_df, "R_OUTPUT/base_models_by_storage.RDS")

```

-   Add whether phenotype is robust

```{r}
base_model_summary_frozen<-subset(base_model_storage_k_df, Type == "FROZEN")
base_model_summary_freezedried<-subset(base_model_storage_k_df, Type == "FREEZEDRIED")

frozen_df<-base_model_summary_frozen[,c("Taxa", "Variable", "P.val")]
names(frozen_df)[3]<-"P.val_frozen"
freezedried_df<-base_model_summary_freezedried[,c("Taxa", "Variable", "P.val")]
names(freezedried_df)[3]<-"P.val_freezedried"


base_model_summary1<- merge(base_model_summary, frozen_df, by = c("Taxa", "Variable"))
base_model_summary2<- merge(base_model_summary1, freezedried_df, by = c("Taxa", "Variable"))

base_model_summary2$Sig_robust<- base_model_summary2$P.val_frozen < 0.05 & base_model_summary2$P.val_freezedried < 0.05


```

### Fig. 5b and c: Fixed-effect models including biological predictors

-   Add biological predictors to models and remove temporal smooths
-   Variance partitioning by comparing R\^2 of full and reduced models.

```{r}

# make empty list to store results
  
full_model_summary_list<-list()


for (i in 1:length(phenotypes)){
  
  tryCatch({ #catch errors
    
    #print(i)
    
    metadata<-data.frame(sample_data(core_genus))
   # str(metadata)
    asvtable<-data.frame(t(otu_table(core_genus)))
    asvtable_log<-  log10(asvtable+1)
    head(asvtable_log)
    names(asvtable_log)<-taxanames
    metadata_all<-merge(metadata, asvtable_log, by = 0)
    row.names(metadata_all)<-metadata_all$feature.id

    
        #scale variables
    metadata_all$Seq_depth_nospikein<-as.numeric(scale(metadata_all$Seq_depth_nospikein))
    metadata_all$HoursTilFrozen<-as.numeric(scale(metadata_all$HoursTilFrozen))
    
    metadata_all$TotalAbundance<-log10(metadata_all$TotalAbundance+1) #log10 transform load
    metadata_all$Observed<-scale(metadata_all$Observed) # scale observed
    metadata_all$Shannon<-scale(metadata_all$Shannon) # scale shannon
    
    metadata_all$Condition_weight<-as.numeric(scale(metadata_all$Condition_weight))
    metadata_all$Temp_min<-as.numeric(scale(metadata_all$Temp_min))
    metadata_all$Temp_max<-as.numeric(scale(metadata_all$Temp_max))
    metadata_all$HourlyTemp<-as.numeric(scale(metadata_all$HourlyTemp))
    metadata_all$HoursSinceForagingStart<-as.numeric(scale(log10(metadata_all$HoursSinceForagingStart+0.1)))
    metadata_all$sum_rainfall_month<-as.numeric(scale(log10(metadata_all$sum_rainfall_month+0.1)))
    metadata_all$SocialStatus<-factor(metadata_all$SocialStatus)
    

# fit gamm minus temporal variables
    m_taxa <- mgcv::gamm(metadata_all[, phenotypes[i]]~
       # s(HoursAfterSunrise, bs = "cr")+
      #  s(month, bs = "cc")+
      #  s(AgeAtSamplingYears, bs = "cr")+
        s(Seq_depth_nospikein, bs = "cr")+
        Storage+
        Seq_run+
        s(IndividID, bs = "re")+
        s(GroupAtSampling, bs = "re")+
        HoursTilFrozen+
        Sex+
        SocialStatus+
        Condition_weight+
        HoursSinceForagingStart+
        HourlyTemp+
        Temp_max+
        Temp_min+
        sum_rainfall_month,
      correlation = corARMA(form = ~ 1|Year, p = 3),
      data=metadata_all,
      family = gaussian)
    
    summary<- summary(m_taxa$gam) 

    
    fixed_terms<-data.frame(extract_fixed(m_taxa$gam))[,c(2,5:7)] ## need to load the function here
    row.names(fixed_terms)<-extract_fixed(m_taxa$gam)$term
    names(fixed_terms)<-c("Effect", "P.val", "lower_CI", "upper_CI")
    
    smoothed_terms<-data.frame(summary$s.table)
    names(smoothed_terms)<-c("lower_CI", "upper_CI", "Effect", "P.val")
    smoothed_terms$lower_CI <-0
    smoothed_terms$upper_CI <-0
    
    results_df<-rbind(fixed_terms, smoothed_terms)
    results_df$R.sq<-summary$r.sq
    results_df$Effect_abs<-abs( results_df$Effect)
    results_df$Variable<-row.names(results_df)
    results_df$Taxa<-phenotypes[i]
    
    
    #### model with only methods
    
    m_taxa_methods <- mgcv::gamm(metadata_all[, phenotypes[i]]~
        s(Seq_depth_nospikein, bs = "cr")+
        Storage+
        Seq_run+
        HoursTilFrozen,
      correlation = corARMA(form = ~ 1|Year, p = 3),
      data=metadata_all,
      family = gaussian)
    
    results_df$R.sq_methods<-summary(m_taxa_methods$gam)$r.sq
    
    ### model with all terms
    
        m_taxa_full <- mgcv::gamm(metadata_all[, phenotypes[i]]~
        s(HoursAfterSunrise, bs = "cr")+
        s(month, bs = "cc")+
        s(AgeAtSamplingYears, bs = "cr")+
        s(Seq_depth_nospikein, bs = "cr")+
        Storage+
        Seq_run+
        s(IndividID, bs = "re")+
        s(GroupAtSampling, bs = "re")+
        HoursTilFrozen+
        Sex+
        SocialStatus+
        Condition_weight+
        HoursSinceForagingStart+
        HourlyTemp+
        Temp_max+
        Temp_min+
        sum_rainfall_month,
      correlation = corARMA(form = ~ 1|Year, p = 3),
      data=metadata_all,
      family = gaussian)
    
    results_df$R.sq_full<-summary(m_taxa_full$gam)$r.sq
    

    full_model_summary_list[[i]]<-results_df
    
  }, error=function(e){})
}

names(full_model_summary_list)<-phenotypes

#saveRDS(full_model_summary_list, "R_OUTPUT/full_models.RDS")

```

```{r}

full_model_summary<-do.call(rbind, full_model_summary_list)

full_model_summary$Significance<-full_model_summary$P.val < 0.01
full_model_summary<-subset(full_model_summary, Variable !="Intercept")

full_model_summary$Variable<- factor(full_model_summary$Variable, levels = c(
  "s(HoursAfterSunrise)", "s(month)" , "s(AgeAtSamplingYears)",   
  "s(IndividID)", "s(GroupAtSampling)",
  "sum_rainfall_month",  "Temp_max", "Temp_min", "HourlyTemp",
  "SexM",   "Condition_weight",  "SocialStatusSubordinate",    "HoursSinceForagingStart", 
  "StorageFROZEN","s(Seq_depth_nospikein)", "Seq_runRUN2", "Seq_runRUN3","Seq_runRUN4", "HoursTilFrozen"))


```

-   split by storage to assess robustness

```{r}

full_model_storage_i<-list()
full_model_storage_k<-list()


for (k in 1:length(ps_list)){
  
  core_genus_storage<- ps_list[[k]]
 # core_genus_storage<- ps_list[[1]]
  
  for (i in 1:length(phenotypes)){
  
  tryCatch({ #catch errors
    
  #  #print(i)
    
    metadata<-data.frame(sample_data(core_genus_storage))
   # str(metadata)
    asvtable<-data.frame(t(otu_table(core_genus_storage)))
    asvtable_log<-  log10(asvtable+1)
    dim(asvtable_log)
    names(asvtable_log)<-taxanames
    metadata_all<-merge(metadata, asvtable_log, by = 0)
    row.names(metadata_all)<-metadata_all$feature.id
    
     #scale variables
    metadata_all$Seq_depth_nospikein<-as.numeric(scale(metadata_all$Seq_depth_nospikein))
    metadata_all$HoursTilFrozen<-as.numeric(scale(metadata_all$HoursTilFrozen))
    
    metadata_all$TotalAbundance<-log10(metadata_all$TotalAbundance+1) #log10 transform load
    metadata_all$Observed<-scale(metadata_all$Observed) # scale observed
    metadata_all$Shannon<-scale(metadata_all$Shannon) # scale shannon
    
    metadata_all$Condition_weight<-as.numeric(scale(metadata_all$Condition_weight))
    metadata_all$Temp_min<-as.numeric(scale(metadata_all$Temp_min))
    metadata_all$Temp_max<-as.numeric(scale(metadata_all$Temp_max))
    metadata_all$HourlyTemp<-as.numeric(scale(metadata_all$HourlyTemp))
    metadata_all$HoursSinceForagingStart<-as.numeric(scale(log10(metadata_all$HoursSinceForagingStart+0.1)))
    metadata_all$sum_rainfall_month<-as.numeric(scale(log10(metadata_all$sum_rainfall_month+0.1)))
    metadata_all$SocialStatus<-factor(metadata_all$SocialStatus)

    
    # fit gamm
    m_taxa <- mgcv::gamm(metadata_all[, phenotypes[i]]~
       # s(HoursAfterSunrise, bs = "cr")+
                           #s(month, bs = "cc")+
                           #s(AgeAtSamplingYears, bs = "cr")+
                           s(Seq_depth_nospikein, bs = "cr")+
                           #Storage+
                           Seq_run+
                           s(IndividID, bs = "re")+
                           s(GroupAtSampling, bs = "re")+
                           HoursTilFrozen+
                           Sex+
                           SocialStatus+
                           Condition_weight+
                           HoursSinceForagingStart+
                           HourlyTemp+
                           Temp_max+
                           Temp_min+
                           sum_rainfall_month,
                         correlation = corARMA(form = ~ 1|Year, p = 3),
                         data=metadata_all,
                         family = gaussian)
    
    summary<- summary(m_taxa$gam) 

    
    fixed_terms<-data.frame(extract_fixed(m_taxa$gam))[,c(2,5:7)] ## need to load the function here
    row.names(fixed_terms)<-extract_fixed(m_taxa$gam)$term
    names(fixed_terms)<-c("Effect", "P.val", "lower_CI", "upper_CI")
    
    smoothed_terms<-data.frame(summary$s.table)
    names(smoothed_terms)<-c("lower_CI", "upper_CI", "Effect", "P.val")
    smoothed_terms$lower_CI <-0
    smoothed_terms$upper_CI <-0
    
    results_df<-rbind(fixed_terms, smoothed_terms)
    results_df$R.sq<-summary$r.sq
    results_df$Effect_abs<-abs( results_df$Effect)
    results_df$Variable<-row.names(results_df)
    results_df$Taxa<-phenotypes[i]
    results_df$Type <- names(ps_list[k])
    
    full_model_storage_i[[i]]<-results_df
    
  }, error=function(e){})
    
  }
    
    full_model_storage_df<-do.call(rbind, full_model_storage_i)
    
    full_model_storage_k[[k]] <-full_model_storage_df
    
}

full_model_storage_k_df<-do.call(rbind, full_model_storage_k)

#saveRDS(full_model_storage_k_df, "R_OUTPUT/full_models_by_storage.RDS")

#table(full_model_storage_k_df$Taxa)

```

-   Add whether phenotype is robust

```{r}

full_model_summary_frozen<-subset(full_model_storage_k_df, Type == "FROZEN")
full_model_summary_freezedried<-subset(full_model_storage_k_df, Type == "FREEZEDRIED")

frozen_df<-full_model_summary_frozen[,c("Taxa", "Variable", "P.val")]
names(frozen_df)[3]<-"P.val_frozen"
freezedried_df<-full_model_summary_freezedried[,c("Taxa", "Variable", "P.val")]
names(freezedried_df)[3]<-"P.val_freezedried"


full_model_summary1<- merge(full_model_summary, frozen_df, by = c("Taxa", "Variable"))
full_model_summary2<- merge(full_model_summary1, freezedried_df, by = c("Taxa", "Variable"))

full_model_summary2$Sig_robust<- full_model_summary2$P.val_frozen < 0.05 & full_model_summary2$P.val_freezedried < 0.05

#full_model_summary2

```

#### Merge all model stats

```{r}

#full_model_df$R.sq_partial<-full_model_df$R.sq - full_model_df$R.sq_methods
full_model_summary2$r2_fixed <-full_model_summary2$R.sq - full_model_summary2$R.sq_methods
full_model_summary2$r2_smooth <-full_model_summary2$R.sq_full - full_model_summary2$R.sq


### order taxa 
full_model_summary2$Taxa <- factor(full_model_summary2$Taxa, levels = c(rev(top_genera[-17]), rev(diversity_names)))
unique(full_model_summary2$Taxa)
#full_model_summary2$Taxa <- factor(full_model_summary2$Taxa, levels = c(rev(top_genera[-16])))

### add third significance group

full_model_summary2$Sig_robust <-ifelse(full_model_summary2$Sig_robust == FALSE & full_model_summary2$P.val < 0.05, "Overall significant", as.character(full_model_summary2$Sig_robust))
full_model_df<-full_model_summary2

## base models 
base_model_summary2$Taxa <- factor(base_model_summary2$Taxa, levels = c(rev(top_genera[-17]), rev(diversity_names)))

### add third significance group

base_model_summary2$Sig_robust <-ifelse(base_model_summary2$Sig_robust == FALSE & base_model_summary2$P.val < 0.05, "Overall significant", as.character(base_model_summary2$Sig_robust))

base_model_summary_df<-base_model_summary2



```

#### Plotting estimates

```{r fig.height=6, fig.width=13, message=FALSE, warning=FALSE}


fixed_df<-subset(full_model_df, 
              #   Variable != "s(HoursAfterSunrise)" & 
              #     Variable != "s(month)" & 
              #     Variable != "s(AgeAtSamplingYears)"& 
                   Variable != "s(IndividID)"& 
                  Variable != "s(GroupAtSampling)" &
                  Variable != "s(Seq_depth_nospikein)"& 
                   Variable != "Storage" &
                   Variable != "Seq_runRUN2" &
                   Variable != "Seq_runRUN3" &
                   Variable != "Seq_runRUN4" &
                   Variable != "StorageFROZEN" &
                   Variable != "HoursTilFrozen")

fixed_df$Variable<-factor(fixed_df$Variable)
unique(fixed_df$Variable)


levels(fixed_df$Variable) <- c("Rainfall", "Max. temp.", "Min. temp.", 
                              "Sampling temp.", 
                               "Sex (Male)", "Body condition", "Social status (Sub.)",  "Time foraging")


################


smooth_df<-subset(base_model_summary_df, Variable == "s(HoursAfterSunrise)" |  Variable == "s(month)" |  Variable == "s(AgeAtSamplingYears)")

smooth_df$Variable<-factor(smooth_df$Variable, levels = c("s(HoursAfterSunrise)", "s(month)", "s(AgeAtSamplingYears)"))

levels(smooth_df$Variable) <-c("Hours after sunrise", "Month", "Age (years)")



P1<- ggplot(smooth_df, aes(x = Effect, y = Taxa, fill = Sig_robust))+
  theme_bw(base_size = 12)+
  #geom_errorbarh(aes(xmin=upper_CI, xmax=lower_CI, height = 0), alpha = 0.5, col = "black", size = 1)+
  geom_point( pch = 21, alpha = 0.9, col = "black", size = 2)+
  facet_wrap(~Variable, nrow = 1, labeller = label_wrap_gen(width=10))+
  scale_size(range = c(1.5,4))+
  scale_fill_manual(values = c("darkgray" , "deepskyblue", "blue"))+
  scale_color_manual(values = c("darkgray", "deepskyblue", "blue"))+
 theme(legend.position = "none") +
  theme(axis.title = element_blank())+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme( # remove the vertical grid lines
    panel.grid.major.x = element_blank() ,
    panel.grid.minor.x = element_blank() ,
    # explicitly set the horizontal lines (or they will disappear too)
    panel.grid.major.y = element_line( size=.05, color="lightgrey" )) + 
  theme(  plot.margin = margin(0.2, 0.1,0.5,0, "cm"))




P2<-ggplot(fixed_df, aes(x = Effect, y = Taxa, fill = Sig_robust))+
  theme_bw(base_size = 12)+
  geom_vline(xintercept = 0, alpha = 0.3, size = 0.3)+
  geom_errorbarh(aes(xmax=upper_CI+0.08, xmin=lower_CI-0.08, height = 0, col = Sig_robust), size = 1)+
  geom_point(size = 2, pch = 21, alpha = 0.9, col = "black")+
  facet_wrap(~Variable, nrow = 1, labeller = label_wrap_gen(width=10), scales = "free_x")+
  scale_size(range = c(1.5,3))+
  scale_fill_manual(values = c("darkgray","pink" , "red"))+
  scale_color_manual(values = c("darkgray","pink" , "red"))+
 theme(legend.position = "none")+
 theme(axis.title.y=element_blank(),
  axis.text.y=element_blank(),
   axis.ticks.y=element_blank())+
  theme(axis.title = element_blank())+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme( # remove the vertical grid lines
    panel.grid.major.x = element_blank() ,
    panel.grid.minor.x = element_blank() ,
    # explicitly set the horizontal lines (or they will disappear too)
    panel.grid.major.y = element_line( size=.05, color="lightgrey" ))


## Plot R2 contributions of methods and biological variables

r2_partial_df<-distinct(full_model_df, Taxa, .keep_all = T)
head(r2_partial_df)
R2_df<-gather(r2_partial_df, Type, R2, c("R.sq_methods","r2_fixed", "r2_smooth"), factor_key=TRUE)

R2_df$R2<-ifelse(R2_df$R2 <0, 0, R2_df$R2)

P3<-ggplot(R2_df, aes(x = R2, y = Taxa, fill = Type))+
  geom_bar(stat = "identity")+
  scale_fill_manual(values = c("grey", "red", "dodgerblue"))+
  theme_bw()+
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank())+
  theme(legend.position = "none")+
  theme(  plot.margin = margin(1.6, 0.2,0.4,0.1, "cm"))+
 # geom_vline(xintercept = 0.1, linetype = "dashed", size = 1)+
   theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
         axis.title.x = element_blank(),
         panel.grid = element_blank())


ggarrange(P1,P2,P3, ncol = 3, widths = c(3.5,5,1.5))


#ggsave("FIGURES/Figure5.pdf", width = 13, height = 6, units = "in")

```

\newpage

# Figure 6: Do diurnal oscillations decay with age?

-   Split age into categories (young, old, and very old).
-   Fit a GAMM to every genus, with time of day smooths calculated per age category.
-   Extract confidence intervals.

```{r}


taxanames<-as.character(taxa_names(core_genus))
confint_list_age<-list()
#model_stats_list<-list()


for (i in 1:length(taxanames)){
  
  tryCatch({ #catch errors
    
    #print(i)
    #taxa1<-prune_taxa("Lachnoclostridium", core_genus)
    taxa1<-prune_taxa(taxanames[i], core_genus)
    sample_data(taxa1)$taxa_i_abundance<-sample_sums(taxa1)
    metadata<-data.frame(sample_data(taxa1))
    metadata$Seq_depth_nospikein<-as.numeric(scale(metadata$Seq_depth_nospikein))
    
    m_taxa <- mgcv::gamm(log10(taxa_i_abundance+1)~
        s(HoursAfterSunrise, bs = "cr", by = AgeCat)+
        s(month, bs = "cc")+
        # s(AgeAtSamplingYears, bs = "cr")+
        s(Seq_depth_nospikein, bs = "cr")+
        Storage+
        Seq_run+
        HoursTilFrozen+
        s(IndividID, bs = "re")+
        s(GroupAtSampling, bs = "re"),
      correlation = corARMA(form = ~ 1|Year, p = 3),
      data=metadata,
      family = gaussian)
    
    confint_df_time<-confint(m_taxa$gam, parm = "s(HoursAfterSunrise)",  type = "confidence")
    confint_df_time$Taxa<-taxanames[i]
    names(confint_df_time)[3]<-"Time"
    confint_list_age[[i]]<-confint_df_time
    
  }, error=function(e){})
}

names(confint_list_age)<-taxanames

```

```{r}

confint_summary_age<-do.call(rbind, confint_list_age)

unique(confint_summary_age$Taxa)

unique(confint_summary_age$AgeCat)
confint_summary_age$AgeCat<-factor(confint_summary_age$AgeCat, levels = c("YOUNG", "MIDDLEAGED", "OLD"))
levels(confint_summary_age$AgeCat)<-c("Young", "Adult", "Old")

abundance_df<-data.frame(taxa_sums(core_genus))
names(abundance_df)<-"Abundance"
abundance_df$Genus<-row.names(abundance_df)

confint_summary_age$Abundance<-vlookup(confint_summary_age$Taxa, abundance_df, lookup_column = "Genus", result_column = "Abundance")
confint_summary_age<-arrange(confint_summary_age, desc(Abundance))
confint_summary_age$Taxa<-factor(confint_summary_age$Taxa, levels = rev(top_genera[-17]))
```

#### Building figure

```{r fig5abc, message = FALSE, warning = FALSE, fig.width = 10}

sample_data(scaled_core)$AgeCat<-factor(sample_data(scaled_core)$AgeCat, levels = c("YOUNG", "MIDDLEAGED", "OLD"))
levels(sample_data(scaled_core)$AgeCat)<-c("Young", "Adult", "Old")

age_histogram<-ggplot(sample_data(scaled_core), aes(x = HoursAfterSunrise))+
  geom_histogram(fill = "grey", col = "black", binwidth = 0.5)+
  facet_wrap(~AgeCat, scales = "free_y")+
  theme_minimal(base_size = 14)+
  theme( strip.background = element_blank(),
    strip.text.x = element_blank())+
  theme(plot.margin = unit(c(0.2,7.4, 0,0.3), "cm"))+
  theme(axis.title.x = element_blank(), 
        axis.text.x = element_blank())+
  ylab("Frequency")+
   theme( # remove the vertical grid lines
    panel.grid.major.x = element_blank() ,
    panel.grid.minor.x = element_blank() ,
    # explicitly set the horizontal lines (or they will disappear too)
    panel.grid.major.y = element_line( size=.05, color="lightgrey" ))


####################


confint_summary_age_sig<-subset(confint_summary_age, Taxa != "Enterococcus" & Taxa != "Lactococcus" & Taxa != "Eubacterium")

unique(confint_summary_age$Taxa)

Fig6abc<-ggplot(confint_summary_age_sig, aes(x = Time, y = est, col = Taxa))+
  facet_wrap(~AgeCat)+
   geom_line( aes(col = Taxa, size = Abundance))+
  scale_size(range = c(1, 2))+
  ylim(-1.6,1.1)+
  scale_color_manual(values = rev(mypaly[c(-4, -10, -16, -17)]), labels = top_genera[c(-4, -10, -16, -17)])+
  theme_bw(base_size = 14)+
  ylab("Model estimate (log10)")+
  xlab("Hours after sunrise")+
  theme(plot.title = element_text(size=10))+
  guides(size=FALSE)+
  #guides(color = guide_legend(override.aes = list(size = 3))) +
  theme(strip.background = element_blank(),
        strip.text.x = element_blank())+
 # theme(legend.position = "none")+
  geom_hline(yintercept = 0, linetype = "dashed") +
  guides(color = guide_legend(override.aes = list( size = 4, col = mypaly[c(-4, -10, -16, -17)], 
     labels = top_genera[c(-4, -10, -16, -17)])))+
   labs(color='Genus')+
  theme(plot.margin = unit(c(0,0.2, 0.2,0.3), "cm"))

Fig6top<-ggarrange(age_histogram, Fig6abc, ncol = 1, heights = c(0.5,2))



```

### A focus on five genera

-   Look more closely at Clostridium, Bacteroides, Paeniclostridium, Cellulomonas, and Raoultibacter

```{r age-interaction-focus, fig.height=8, fig.width=13, message=FALSE, warning=FALSE}
clostridium_age<-subset(confint_summary_age, Taxa == "Clostridium sensu stricto 1")

clostridium<-ggplot(clostridium_age, aes(x = Time, y = est, col = AgeCat))+
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = AgeCat), linetype=2, alpha=0.1)+
  coord_cartesian(ylim = c(-2,1.5))+
  scale_color_manual(values = c( "blue", "darkgrey", "red"))+
  scale_fill_manual(values = c( "blue", "darkgrey", "red"))+
  theme_bw(base_size = 14)+
  ylab("Model estimate (log10)")+
  xlab("")+
  # theme(legend.position = "none")+
  ggtitle("d) Clostridium")+
  theme(plot.title = element_text(size=12))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(col='Age category')+
    guides(fill = FALSE)


bacteroides_age<-subset(confint_summary_age, Taxa == "Bacteroides")

bacteroides<-ggplot(bacteroides_age, aes(x = Time, y = est, col = AgeCat))+
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = AgeCat), linetype=2, alpha=0.1)+
  coord_cartesian(ylim = c(-2,1.5))+
  scale_color_manual(values = c( "blue", "darkgrey", "red"))+
  scale_fill_manual(values = c( "blue", "darkgrey", "red"))+
  theme_bw(base_size = 14)+
  theme(axis.title.y = element_blank(), axis.text.y = element_blank())+
  xlab("")+
 # theme(legend.position = "none")+
  ggtitle("e) Bacteroides")+
  theme(plot.title = element_text(size=12)) +
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(col='Age category')+
    guides(fill = FALSE)

blautia_age<-subset(confint_summary_age, Taxa == "Blautia")

blautia<-ggplot(blautia_age, aes(x = Time, y = est, col = AgeCat))+
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = AgeCat), linetype=2, alpha=0.1)+
  coord_cartesian(ylim = c(-2,1.5))+
  scale_color_manual(values = c( "blue", "darkgrey", "red"))+
  scale_fill_manual(values = c( "blue", "darkgrey", "red"))+
  theme_bw(base_size = 14)+
  theme(axis.title.y = element_blank(), axis.text.y = element_blank())+
  xlab("Hours after sunrise")+
  # theme(legend.position = "none")+
  ggtitle("f) Blautia")+
  theme(plot.title = element_text(size=12))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(col='Age category')+
    guides(fill = FALSE)

ruminococcus_age<-subset(confint_summary_age, Taxa == "[Ruminococcus] torques group")

ruminococcus<-ggplot(ruminococcus_age, aes(x = Time, y = est, col = AgeCat))+
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = AgeCat), linetype=2, alpha=0.1)+
  coord_cartesian(ylim = c(-2,1.5))+
  scale_color_manual(values = c( "blue", "darkgrey", "red"))+
  scale_fill_manual(values = c( "blue", "darkgrey", "red"))+
  theme_bw(base_size = 14)+
  theme(axis.title.y = element_blank(), axis.text.y = element_blank())+
  xlab("Hours after sunrise")+
  # theme(legend.position = "none")+
  ggtitle("f) Ruminococcus")+
  theme(plot.title = element_text(size=12))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(col='Age category')+
    guides(fill = FALSE)

paeniclostridium_age<-subset(confint_summary_age, Taxa == "Paeniclostridium")

paeniclostridium<-ggplot(paeniclostridium_age, aes(x = Time, y = est, col = AgeCat))+
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = AgeCat), linetype=2, alpha=0.1)+
  coord_cartesian(ylim = c(-2,1.5))+
  scale_color_manual(values = c( "blue", "darkgrey", "red"))+
  scale_fill_manual(values = c( "blue", "darkgrey", "red"))+
  theme_bw(base_size = 14)+
  theme(axis.title.y = element_blank(), axis.text.y = element_blank())+
  xlab("")+
  # theme(legend.position = "none")+
  ggtitle("f) Paeniclostridium")+
  theme(plot.title = element_text(size=12))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(col='Age category')+
    guides(fill = FALSE)



cellulomonas_age<-subset(confint_summary_age, Taxa == "Cellulomonas")

cellulomonas<-ggplot(cellulomonas_age, aes(x = Time, y = est, col = AgeCat))+
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = AgeCat), linetype=2, alpha=0.1)+
  coord_cartesian(ylim = c(-2,1.5))+
  scale_color_manual(values = c( "blue", "darkgrey", "red"))+
  scale_fill_manual(values = c( "blue", "darkgrey", "red"))+
  theme_bw(base_size = 14)+
  theme(axis.title.y = element_blank(), axis.text.y = element_blank())+
  xlab("")+
  # theme(legend.position = "none")+
  ggtitle("g) Cellulomonas")+
  theme(plot.title = element_text(size=12))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(col='Age category')+
    guides(fill = FALSE)

raoultibacter_age<-subset(confint_summary_age, Taxa == "Raoultibacter")

raoultibacter<-ggplot(raoultibacter_age, aes(x = Time, y = est, col = AgeCat))+
  geom_line(size = 1.5)+
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = AgeCat), linetype=2, alpha=0.1)+
  coord_cartesian(ylim = c(-2,1.5))+
  scale_color_manual(values = c( "blue", "darkgrey", "red"))+
  scale_fill_manual(values = c( "blue", "darkgrey", "red"))+
  theme_bw(base_size = 14)+
  theme(axis.title.y = element_blank(), axis.text.y = element_blank())+
  xlab("")+
  # theme(legend.position = "none")+
  ggtitle("h) Raoultibacter")+
  theme(plot.title = element_text(size=12))+
  geom_hline(yintercept = 0, linetype = "dashed")+
  labs(col='Age category')+
    guides(fill = FALSE)


Fig6def<-ggarrange(clostridium, bacteroides,paeniclostridium, cellulomonas, raoultibacter, ncol = 5,  common.legend = TRUE, legend = "right", widths = c(1.2,1,1,1,1))+
  theme(plot.margin = unit(c(0.2,3.5,0.2,0.2), "cm"))


ggarrange(Fig6top, Fig6def, ncol = 1, heights = c(1.5,1))
#ggsave("FIGURES/Figure6.pdf", width = 13, height = 8, units = "in")

```

# Technical replication

-   Import data for 16 technical replicates
-   Calculate load and diversity measures
-   Run linear models with sample ID as random effect
-   Estimate proportion of variation attributed to sample ID using the icc() function

```{r fig.height=8, fig.width=10}

meerkat_replicates<-readRDS("technical_replicate_data_phyloseq.RDS")

meerkat_replicates

### add total abundance of reads to dataframe
sample_data(meerkat_replicates)$TotalAbundance<-sample_sums(meerkat_replicates)

# estimate alpha diversity
sample_data(meerkat_replicates)$Observed<-phyloseq::estimate_richness(meerkat_replicates, measures="Observed")
sample_data(meerkat_replicates)$Observed<-sample_data(meerkat_replicates)$Observed$Observed

sample_data(meerkat_replicates)$Shannon<-phyloseq::estimate_richness(meerkat_replicates, measures="Shannon")
sample_data(meerkat_replicates)$Shannon<-sample_data(meerkat_replicates)$Shannon$Shannon

####### beta diversity 

# weighted unifrac

ord_wunifrac_mds<-ordinate(meerkat_replicates, "MDS", "wunifrac")
wunifrac_dist<-distance(meerkat_replicates, method = "wunifrac")

sample_data(meerkat_replicates)$MDS1_WU<-data.frame(ord_wunifrac_mds$vectors)$Axis.1
sample_data(meerkat_replicates)$MDS2_WU<-data.frame(ord_wunifrac_mds$vectors)$Axis.2


# unweighted

ord_uwunifrac_mds<-ordinate(meerkat_replicates, "MDS", "unifrac")
unifrac_dist<-distance(meerkat_replicates, method = "unifrac")


sample_data(meerkat_replicates)$MDS1_UU<-data.frame(ord_uwunifrac_mds$vectors)$Axis.1
sample_data(meerkat_replicates)$MDS2_UU<-data.frame(ord_uwunifrac_mds$vectors)$Axis.2

## jaccard

ord_jaccard_mds<-ordinate(meerkat_replicates, "MDS", "jaccard")
jaccard_dist<-distance(meerkat_replicates, method = "jaccard")


sample_data(meerkat_replicates)$MDS1_JC<-data.frame(ord_jaccard_mds$vectors)$Axis.1
sample_data(meerkat_replicates)$MDS2_JC<-data.frame(ord_jaccard_mds$vectors)$Axis.2

## morisita


ord_morisita_mds<-ordinate(meerkat_replicates, "MDS", "morisita")
morisita_dist<-distance(meerkat_replicates, method = "morisita")


sample_data(meerkat_replicates)$MDS1_Mor<-data.frame(ord_morisita_mds$vectors)$Axis.1
sample_data(meerkat_replicates)$MDS2_Mor<-data.frame(ord_morisita_mds$vectors)$Axis.2


replicates_df<-data.frame(sample_data(meerkat_replicates))


load_model_re<-lmer(TotalAbundance~1+(1|SampleNo), data = replicates_df)
observed_model_re<-lmer(Observed~1 + (1|SampleNo), data = replicates_df)
shannon_model_re<-lmer(Shannon~1 + (1|SampleNo), data = replicates_df)
MDS1_WU_re<-lmer(MDS1_WU~1 + (1|SampleNo), data = replicates_df)
MDS1_UU_re<-lmer(MDS1_UU~1 + (1|SampleNo), data = replicates_df)
MDS1_JC_re<-lmer(MDS2_JC~1 + (1|SampleNo), data = replicates_df)
MDS1_Mor_re<-lmer(MDS1_Mor~1 + (1|SampleNo), data = replicates_df)

## calculate contribution of sample on all diversity measures


performance::icc(load_model_re)
performance::icc(observed_model_re)
performance::icc(shannon_model_re)
performance::icc(MDS1_WU_re)
performance::icc(MDS1_UU_re)
performance::icc(MDS1_JC_re)
performance::icc(MDS1_Mor_re)


adonis_weighted<-vegan::adonis(wunifrac_dist ~ SampleNo, data = replicates_df)
adonis_unweighted<-vegan::adonis(unifrac_dist ~ SampleNo, data = replicates_df)
adonis_jaccard<-vegan::adonis(jaccard_dist ~ SampleNo, data = replicates_df)
adonis_morisita<-vegan::adonis(morisita_dist ~ SampleNo, data = replicates_df)

adonis_weighted
adonis_unweighted
adonis_jaccard
adonis_morisita

####################### plot  #########################
####################### plot  #########################
####################### plot  #########################

### colours

mypal =  	pal_lancet("lanonc", alpha = 1)(7)
mypal5<- pal_rickandmorty("schwifty")(8)
mypal6<-brewer.pal(8,"BrBG")[-1]
mypalx<-c(mypal, mypal5, mypal6)


rep_load<-ggplot(sample_data(meerkat_replicates), aes(x = fct_reorder(SampleNo, TotalAbundance, .fun = "median"), 
                                        y = TotalAbundance, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_jitter( width = 0, height = 0.01, size = 4, aes(fill = SampleNo), col = "black")+
  ylab("Bacterial load estimate")+
  xlab("Sample")+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  theme(axis.text.x = element_blank())+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))

#### alpha diversity

rep_observed<-ggplot(sample_data(meerkat_replicates), aes(x = fct_reorder(SampleNo, Observed, .fun = "median"), 
                                                          y = Observed, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_jitter(width = 0, height = 0.01, size = 4, aes(fill = SampleNo), col = "black")+
  ylab("Observed richness")+
  xlab("Sample")+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  theme(axis.text.x = element_blank())+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))

rep_shannon<-ggplot(sample_data(meerkat_replicates), aes(x = fct_reorder(SampleNo, Shannon, .fun = "median"), 
                                                          y = Shannon, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_jitter( width = 0, height = 0.01, size = 4, aes(fill = SampleNo), col = "black")+
  ylab("Shannon index")+
  xlab("Sample")+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  theme(axis.text.x = element_blank())+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))

### axis 

rep_MDS1_WU<-ggplot(sample_data(meerkat_replicates), aes(x = fct_reorder(SampleNo, MDS1_WU, .fun = "median"), 
                                                         y =  MDS1_WU, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_jitter( width = 0, height = 0.01, size = 4, aes(fill = SampleNo), col = "black")+
  ylab("Axis 1 Weighted Unifrac")+
  xlab("Sample")+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  theme(axis.text.x = element_blank())+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))



rep_MDS1_UU<-ggplot(sample_data(meerkat_replicates), aes(x = fct_reorder(SampleNo, MDS1_UU, .fun = "median"), 
                                                          y =  MDS1_UU, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_jitter( width = 0, height = 0.01, size = 4, aes(fill = SampleNo), col = "black")+
  ylab("Axis 1 Unweighted Unifrac")+
  xlab("Sample")+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  theme(axis.text.x = element_blank())+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))


rep_MDS1_JC<-ggplot(sample_data(meerkat_replicates), aes(x = fct_reorder(SampleNo, MDS1_JC, .fun = "median"), 
                                                         y =  MDS1_JC, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_jitter( width = 0, height = 0.01, size = 4, aes(fill = SampleNo), col = "black")+
  ylab("Axis 1 Jaccard")+
  xlab("Sample")+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  theme(axis.text.x = element_blank())+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))



rep_MDS1_Mor<-ggplot(sample_data(meerkat_replicates), aes(x = fct_reorder(SampleNo, MDS1_Mor, .fun = "median"), 
                                                         y =  MDS1_Mor, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_jitter( width = 0, height = 0.01, size = 4, aes(fill = SampleNo), col = "black")+
  ylab("Axis 1 Morisita")+
  xlab("Sample")+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  theme(axis.text.x = element_blank())+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))

ggarrange(rep_load, rep_observed, 
            rep_MDS1_JC, 
          rep_MDS1_Mor, rep_MDS1_UU, 
          rep_MDS1_WU,
         # rep_weighted, rep_unweighted, rep_jaccard,
          nrow = 2, ncol = 3, 
          common.legend = T, legend = "right", 
          labels = c("a)", "b)", "c)", "d)", "e)", "f)"), 
          font.label = list(size = 20))



###### ordination plots


rep_weighted<-ggplot(sample_data(meerkat_replicates), aes(x = MDS1_WU, y = MDS2_WU, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_point(aes(fill = SampleNo), size = 4)+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))+
  xlab("Axis 1 WU (36%)")+
  ylab("Axis 2 WU (13%)")


# unweighted


rep_unweighted<-ggplot(sample_data(meerkat_replicates), aes(x = MDS1_UU, y = MDS2_UU, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_point(aes(fill = SampleNo), size = 4)+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))+
  xlab("Axis 1 UU (19%)")+
  ylab("Axis 2 UU (11%)")


## jaccard

ord_jaccard_mds$values$Relative_eig[1:4]

rep_jaccard<-ggplot(sample_data(meerkat_replicates), aes(x = MDS1_JC, y = MDS2_JC, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_point(aes(fill = SampleNo), size = 4)+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))+
  xlab("Axis 1 JC (14%)")+
  ylab("Axis 2 JC (9%)")

## morisita

ord_morisita_mds$values$Relative_eig[1:4]

rep_morisita<-ggplot(sample_data(meerkat_replicates), aes(x = MDS1_Mor, y = MDS2_Mor, group = SampleNo, shape = Storage))+
  geom_line(aes(col = SampleNo), size = 2)+
  geom_point(aes(fill = SampleNo), size = 4)+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)+
  scale_shape_manual(values = c(21,24))+
  guides(fill = guide_legend(override.aes=list(shape=21)))+
  xlab("Axis 1 Morisita (24%)")+
  ylab("Axis 2 Morista (19%)")

ggarrange(rep_jaccard, rep_morisita, rep_unweighted, rep_weighted, ncol = 2, nrow = 2, common.legend = T, legend = "right")




```

# Pilot study

-   18 samples from 9 captive meerkats
-   Samples were either frozen or freezedried

```{r pilot-study, fig.width= 6, fig.height = 4}
## pilot study 

pilot<- readRDS("pilot_study_data_phyloseq.RDS")

pilot


pilot_rarefied<-rarefy_even_depth(pilot, sample.size = min(sample_sums(pilot)), rngsee = 100, replace = TRUE, trimOTUs=TRUE,verbose=TRUE)

core_microbiome<-speedyseq::tax_glom(pilot_rarefied, taxrank = "Family")
core_microbiome<-microbiome::core(core_microbiome, detection = 500, prevalence = 0.6)
speedyseq::plot_bar(core_microbiome, fill = "Family")+scale_fill_manual(values = mypaly)+facet_wrap(~Storage, scales = "free")


### beta diversity 

ord_wunifrac_mds<-ordinate(pilot_rarefied, "MDS", "wunifrac")
ord_unifrac_mds<-ordinate(pilot_rarefied, "MDS", "unifrac")

dist_wunifrac<-distance(pilot_rarefied, method="wunifrac")
dist_unifrac<-distance(pilot_rarefied, method="unifrac")

sample_data(pilot_rarefied)$MDS1<-data.frame(ord_wunifrac_mds$vectors)$Axis.1
sample_data(pilot_rarefied)$MDS2<-data.frame(ord_wunifrac_mds$vectors)$Axis.2
sample_data(pilot_rarefied)$MDS3<-data.frame(ord_wunifrac_mds$vectors)$Axis.3
sample_data(pilot_rarefied)$MDS1_UU<-data.frame(ord_unifrac_mds$vectors)$Axis.1
sample_data(pilot_rarefied)$MDS2_UU<-data.frame(ord_unifrac_mds$vectors)$Axis.2


ggplot(sample_data(pilot_rarefied), aes(x = MDS1, y = MDS2, group = IndividID, fill = Storage))+
  geom_line(size = 2)+
  geom_point( size = 4, pch = 21)+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)

ggplot(sample_data(pilot_rarefied), aes(x = MDS1_UU, y = MDS2_UU, group = IndividID, fill = Storage))+
  geom_line(size = 2)+
  geom_point( size = 4, pch = 21)+
  theme_bw(base_size = 14)+
  scale_fill_manual(values = mypalx)+
  scale_color_manual(values = mypalx)


metadata<-data.frame(sample_data(pilot_rarefied))

adonis2(dist_wunifrac~
    Storage+
    IndividID, 
    data = metadata,
     by = "margin")


adonis2(dist_unifrac~
          Storage+
          IndividID, 
        data = metadata,
        by = "margin")

```

# How long did this script take?

-   Time it look to run this script

```{r}
end_time<-Sys.time()

round((end_time - start_time), 3)

```

# Acknowledgements

Thanks to those who made their data and code open access for myself and others to copy and build upon. Also thanks to the many blog authors for providing explanations of code and stats that helped to build this analysis. Special thanks to Gavin Simpson for his blogs on GAMs, gratia and vegan (https://fromthebottomoftheheap.net/), which I found especially useful throughout this analysis.

# Do me a favour
If you are using this data and/or code please pay it forward and put your own data and code online in a format that is easily reproducible. People who use the line 'data/code available on reasonable request' are takers - don't be that person. I am a researcher that supports equal opportunities for everyone, and stand in solidarity with those negatively affected by political systems and power structures that promote sexist, racist, colonial, and homophobic attitudes.