7 Beta diversity

load("data/data.Rdata")

7.1 Diversity calculations

beta_div <- list()

n=0
for (strategy in strategies){
  n=n+1
  bin_counts_norm_filt_strategy <- bin_counts_norm_filt %>%
    filter(overall_strategy == {{strategy}}) %>%
    select(secondary_cluster,sample,read_count) %>%
    pivot_wider(names_from = "sample", values_from = "read_count", values_fn = sum) %>%
    column_to_rownames(var="secondary_cluster")
    
   cluster_tree_strategy <- keep.tip(cluster_tree, tip=rownames(bin_counts_norm_filt_strategy))
  
    beta_q0n <- bin_counts_norm_filt_strategy %>%
      hillpair(., q=0, metric="C") %>% as.dist()
    
    beta_q1n <- bin_counts_norm_filt_strategy %>%
      hillpair(., q=1, metric="C")
    
    beta_q1p <- bin_counts_norm_filt_strategy %>%
      hillpair(., q=1, metric="C", tree = cluster_tree_strategy)
  
  beta_div_strategy <- list(q0n=beta_q0n,q1n=beta_q1n,q1p=beta_q1p)
  beta_div[[n]] <- beta_div_strategy
}

names(beta_div) <- strategies
save(beta_div,file="results/beta_diversity.Rdata")

7.2 Procrustes analysis

7.2.1 Richness

load("results/beta_diversity.Rdata")

#Generate NMDS ordinations
nmds_q0n <- list()
for(strategy in strategies){
  set.seed(1)
  nmds_q0n[[strategy]] <- beta_div[[strategy]]$q0n %>%
  metaMDS(.,trymax = 999, k=2, trace=0)
}

#Procrustes analysis 
combinations <- expand.grid(seq(1, 10), seq(1, 10))
protest_q0n <- tibble(
  v1 = character(),
  v2 = character(),
  scale = numeric(),
  sig = numeric()
)

for(i in c(1:100)){
  x<-combinations[i,1]
  y<-combinations[i,2]
  protest <- protest(nmds_q0n[[x]], nmds_q0n[[y]])
  protest_q0n <- bind_rows(protest_q0n,tibble(v1=names(nmds_q0n)[x],v2=names(nmds_q0n)[y],scale=protest$scale,sig=protest$signif))
}

colfunc<-colorRampPalette(c("#44ce1b","#bbdb44","#f7e379","#f2a134","#e51f1f"))
protest_q0n %>%
  ggplot(aes(x=v1,y=v2,fill=sig)) +
    geom_tile(colour = "white") +
    scale_fill_gradientn(colours=colfunc(10)) +
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x=NULL,y=NULL)

7.2.2 Neutral

load("results/beta_diversity.Rdata")

#Generate NMDS ordinations
nmds_q1n <- list()
for(strategy in strategies){
  set.seed(1)
  nmds_q1n[[strategy]] <- beta_div[[strategy]]$q1n %>%
  metaMDS(.,trymax = 999, k=2, trace=0)
}

#Procrustes analysis 
combinations <- expand.grid(seq(1, 10), seq(1, 10))
protest_q1n <- tibble(
  v1 = character(),
  v2 = character(),
  scale = numeric(),
  sig= numeric()
)

for(i in c(1:100)){
  x<-combinations[i,1]
  y<-combinations[i,2]
  protest <- protest(nmds_q1n[[x]], nmds_q1n[[y]])
  protest_q1n <- bind_rows(protest_q1n,tibble(v1=names(nmds_q1n)[x],v2=names(nmds_q1n)[y],scale=protest$scale,sig=protest$signif))
}

colfunc<-colorRampPalette(c("#44ce1b","#bbdb44","#f7e379","#f2a134","#e51f1f"))
protest_q1n %>%
  ggplot(aes(x=v1,y=v2,fill=sig)) +
    geom_tile(colour = "white") +
    scale_fill_gradientn(colours=colfunc(10)) +
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x=NULL,y=NULL)

7.2.3 Phylogenetic

load("results/beta_diversity.Rdata")

#Generate NMDS ordinations
nmds_q1p <- list()
for(strategy in strategies){
  set.seed(1)
  nmds_q1p[[strategy]] <- beta_div[[strategy]]$q1p %>%
  metaMDS(.,trymax = 999, k=2, trace=0)
}

#Procrustes analysis 
combinations <- expand.grid(seq(1, 10), seq(1, 10))
protest_q1p <- tibble(
  v1 = character(),
  v2 = character(),
  scale = numeric(),
  sig= numeric()
)

for(i in c(1:100)){
  x<-combinations[i,1]
  y<-combinations[i,2]
  protest <- protest(nmds_q1p[[x]], nmds_q1p[[y]])
  protest_q1p <- bind_rows(protest_q1p,tibble(v1=names(nmds_q1p)[x],v2=names(nmds_q1p)[y],scale=protest$scale,sig=protest$signif))
}

colfunc<-colorRampPalette(c("#44ce1b","#bbdb44","#f7e379","#f2a134","#e51f1f"))
protest_q1p %>%
  ggplot(aes(x=v1,y=v2,fill=sig)) +
    geom_tile(colour = "white") +
    scale_fill_gradientn(colours=colfunc(10)) +
    theme_minimal() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x=NULL,y=NULL)