20230208_Halden_analysis.Rmd

---
title: "Microbiome convergence and deterministic community assembly along successional biocrust gradients on potash salt heaps"
author: "Juliette Ohan"
date: "7/28/2022"
output:
  word_document: default
  html_document: default
editor_options: 
  markdown: 
    wrap: 72
---

This is an R document to analyze the chemical values and biological
sequences from Biocrusts on Salt Heaps in Germany aka the "Halden."

Measurements: Abiotics - EC, Chl a, pH, DOC, TDN (site + stage,
3x rep) output:

Biotics - DNA sequencing of 16S rRNA gene (only site + stage, 3x rep)
output: ASV table with taxonomy

qPCR of 16S rRNA bacterial and archaeal genes (site + stage) output:

Co-variates: 2 sites = OD, WT 3 stages = heap, initial, biocrust (3x
replicates)

# READING DATA

```{r}

library(tidyverse)
library(ggplot2)
library(dplyr)

#getwd()
#setwd

#import data from present working directory, if it is a file on your work computer.. THIS IS ALREADY SCALED TO 14347 READS

#Get ASV file

asv.data <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/20220404_with_archaea/abundance_with_archea.csv", 
            header = TRUE, 
           sep = ",",
           dec = ".")

#change first column to row names

asv.data <- data.frame(asv.data, row.names =1)


#get mapping file

salt_map <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/all/salt_map.csv", 
            header = TRUE, 
           sep = ",",
           dec = ".")

#reorder levels

salt_map$type <- factor(salt_map$type, levels = c("heap","initial", "biocrust"))
salt_map$label <- factor(salt_map$label, levels = c ("OD_heap", "OD_initial", "OD_biocrust",
                                                     "WT_heap", "WT_initial", "WT_biocrust"))
#subset by type

salt_map_OD <- subset(salt_map, site == "OD", select = c("NAME", "label", "Sample_num", "seqID", "site", "type"))

salt_map_WT <- subset(salt_map, site == "WT", select = c("NAME", "label","Sample_num", "seqID", "site", "type"))

#match and replace seqID with Names

asv.data2 <- asv.data %>%
      rename_at(as.vector(na.omit(salt_map$seqID[match(names(asv.data), salt_map$seqID)])), 
               ~as.vector(na.omit(salt_map$NAME[match(names(asv.data), salt_map$seqID)])))


```

# PRE-PROCESS

```{r}

#merge biological replicates and remove asvs that are not present in all replicates, so you only have six treatments

#make separate df for each treatment 

OD_crust <- asv.data[, c("Juliette11", "Juliette12", "Juliette13")]
OD_inter <- asv.data[, c("Juliette22", "Juliette15", "Juliette16")]
OD_heap <- asv.data[, c("Juliette17", "Juliette18", "Juliette19")]
WT_crust <- asv.data[, c("Juliette1", "Juliette2", "Juliette3")]
WT_inter <- asv.data[, c("Juliette4", "Juliette5", "Juliette6")]
WT_heap <- asv.data[, c("Juliette10", "Juliette8")]

#removing all rows with ALL zeros 

WT_crust <- WT_crust[apply(WT_crust[,-1], 1, function(x) !all(x==0)),]
WT_inter <- WT_inter[apply(WT_inter[,-1], 1, function(x) !all(x==0)),]
OD_crust <- OD_crust[apply(OD_crust[,-1], 1, function(x) !all(x==0)),]
OD_inter <- OD_inter[apply(OD_inter[,-1], 1, function(x) !all(x==0)),]
OD_heap  <- OD_heap [apply(OD_heap [,-1], 1, function(x) !all(x==0)),]

#different function for WT_heap since it only has two replicates

WT_heap <- WT_heap[rowSums(WT_heap[])>0,]

```

# AMPVIS2 HEATMAP 

```{r dev = c("png", "jpg", "pdf")}

#install.packages("ampvis2")

library(ampvis2)

#load in tables 

data <- amp_load(
  otutable = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/20220404_with_archaea/abundance_with_archea_ampvis.csv ",
  metadata = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/all/OUT_300721/halden_mapping_ampvis.txt",
  taxonomy = NULL,
  fasta = NULL,
  tree = NULL,
  pruneSingletons = TRUE
)

#reorder

data$metadata$type <- factor(data$metadata$type, levels = c("Heap","Initial", "Biocrust"))



#make a heatmap 

order_by_y_vec = paste(data$tax$Phylum, data$tax$Family, sep = "; ")

amp_heatmap(
  data = data,
  tax_aggregate = "Family",
  tax_add = "Phylum",
  tax_show = 50,
  showRemainingTaxa = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  order_y_by = rev(sort(unique(order_by_y_vec))), #call previous variable made 

)


amp_heatmap(
  data = data,
  tax_aggregate = "Phylum",
 # tax_add = "Phylum",
  tax_show = 20,
  showRemainingTaxa = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  #color_vector = c( "whitesmoke", "royalblue4")
  #color_vector = c( "#440154FF", "#481567FF", "#482677FF", "#453781FF","#404788FF","#39568CFF","#33638DFF","#2D708EFF","#287D8EFF","#238A8DFF","#1F968BFF","#20A387FF","#29AF7FFF","#3CBB75FF", "#55C667FF", "#73D055FF", "#95D840FF", "#B8DE29FF", "#DCE319FF", "#FDE725FF")
)


order_by_y_vec = paste(data$tax$Phylum, data$tax$Genus, sep = "; ")

#Genus 
   
  amp_heatmap(
  data = data,
  tax_aggregate = "Genus",
  tax_add = "Phylum",
  tax_show = 40,
  showRemainingTaxa = TRUE,
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  order_y_by = rev(unique(order_by_y_vec)), #call previous variable made 
  plot_colorscale = "sqrt",
  #color_vector = c("royalblue4", "whitesmoke", "darkred")
)


order_by_y_vec = paste(data$tax$Phylum, data$tax$Class, sep = "; ")

sort(order_by_y_vec)


amp_heatmap(
  data = data,
  tax_aggregate = "Class",
  tax_add = "Phylum",
  tax_show = 25,
  showRemainingTaxa = TRUE,
  order_y_by = rev(sort(unique(order_by_y_vec))), #call previous variable made 
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt"
)


```

# AMPVIS2 - CORE VENN DIAGRAMS 

```{r dev = c("png", "jpg", "pdf")}

library(ampvis2)

#venn diagram of core otus

venn_all <- amp_venn(
  data,
  group_by = "type",
  cut_a = 0.001,
  cut_f = 65,
  text_size = 5,
  normalise = TRUE,
  detailed_output = TRUE
)

#between sites

venn_site <- amp_venn(
  data,
  group_by = "site",
  cut_a = .001,
  cut_f = 65,
  text_size = 5,
  normalise = TRUE,
  detailed_output = TRUE
)

#subset only OD 

data_OD <- amp_subset_samples(data, site %in% c("OD"))

venn_OD <- amp_venn(
  data_OD,
  group_by = "type",
  cut_a = 0.001,
  cut_f = 65,
  text_size = 5,
  normalise = TRUE,
  detailed_output = TRUE
)

#subset only WT

data_WT <- amp_subset_samples(data, site %in% c("WT"))

venn_WT <- amp_venn(
  data_WT,
  group_by = "type",
  cut_a = 0.001,
  cut_f = 65,
  text_size = 5,
  normalise = TRUE,
  detailed_output = TRUE
)

#subset crust

data_biocrust <- amp_subset_samples(data, type %in% c("Biocrust"))

venn_biocrust <- amp_venn(
  data_biocrust,
  group_by = "site",
  cut_a = 0.001,
  cut_f = 65,
  text_size = 5,
  normalise = TRUE,
  detailed_output = TRUE
)

#save ASVs to an object 

df_biocrust = venn_biocrust$Otutable

#save df to .csv on local

write.csv(x = df_biocrust, file = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/20220404_with_archaea/df_biocrust.csv", row.names = FALSE)

data_inter <- amp_subset_samples(data, type %in% c("Initial"))

venn_inter <- amp_venn(
 data_inter,
  group_by = "site",
  cut_a = 0.001,
   cut_f = 65,
  text_size = 5,
  normalise = TRUE,
  detailed_output = TRUE
)

#save ASVs to an object 

df_inter = venn_inter$Otutable

#save df to .csv on local

write.csv(x = df_inter, file = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/20220404_with_archaea/df_inter.csv", row.names = FALSE)

#subset heap

data_heap <- amp_subset_samples(data, type %in% c("Heap"))

 venn_heap <- amp_venn(
  data_heap,
  group_by = "site",
  cut_a = 0.001,
  cut_f = 65,
  text_size = 5,
  normalise = TRUE,
  detailed_output = TRUE
)

#save ASVs to an object 
 
df_heap = venn_heap$Otutable
 
#save df to .csv on local

write.csv(x = df_heap, file = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/20220404_with_archaea/df_heap.csv", row.names = FALSE)
```

# AMPVIS - CORE HEATMAPS 

```{r dev = c("png", "jpg", "pdf")}

library(ampvis2)
library(dplyr)
library(readr)

#keep only rows that are "Core"

df_heap2 <- df_heap[df_heap$Shared == "Core", ]

#remove last column (using dplyr)

df_heap2 <- select(df_heap2, -Shared)

#rename row headers (using dplyr)

df_heap2 <- df_heap2 %>%
  rename(
    OTU2 = OTU,
  )

asv.data <- read_csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/20220404_with_archaea/asv.table.csv", col_names = TRUE)

#rename row headers (using dplyr)

asv.data <- asv.data %>%
  rename(
    OTU = ...1,
  )

#Match up values from asv.data to df_heap2

df_heap3 =
  asv.data %>%
  filter(`OTU` %in% df_heap2$OTU2) %>%
  select(c(`OTU`, names(df_heap2) %>% grep("Juliette", ., value = T))) %>%
  right_join(df_heap2[,1:8], by = c("OTU" = "OTU2"))

data_heap <- amp_load(
  otutable = df_heap3,
   metadata = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/all/OUT_300721/halden_mapping_ampvis.txt",
  taxonomy = NULL,
  fasta = NULL,
  tree = NULL,
  pruneSingletons = TRUE
)

#order by phylum

order_by_y_vec = paste(data_heap$tax$Phylum, data_heap$tax$Class, sep = "; ")


plot_heap <- amp_heatmap(
  data = data_heap,
  tax_aggregate = "Class",
  tax_add = "Phylum",
  tax_show = 20,
  showRemainingTaxa = TRUE,
  order_y_by = rev(unique(order_by_y_vec)), #call previous variable made 
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  color_vector = c( "whitesmoke", "royalblue4")
)


#----------------------------------------------------------------------------------------------------------------------------------

#INITIAL CORE HEATMAP

#keep only rows that are "Core"

df_inter2 <- df_inter[df_inter$Shared == "Core", ]

#remove last column (using dplyr)

df_inter2 <- select(df_inter2, -Shared)

#rename row headers (using dplyr)

df_inter2 <- df_inter2 %>%
  rename(
    OTU2 = OTU,
  )

# Match up values from asv.data to df_inter2

df_inter3 =
  asv.data %>%
  filter(`OTU` %in% df_inter2$OTU2) %>%
  select(c(`OTU`, names(df_inter2) %>% grep("Juliette", ., value = T))) %>%
  right_join(df_inter2[,1:8], by = c("OTU" = "OTU2"))


data_inter <- amp_load(
  otutable = df_inter3,
   metadata = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/all/OUT_300721/halden_mapping_ampvis.txt",
  taxonomy = NULL,
  fasta = NULL,
  tree = NULL,
  pruneSingletons = TRUE
)

#order by phylum 

order_by_y_vec = paste(data_inter$tax$Phylum, data_inter$tax$Class, sep = "; ")

plot_inter <- amp_heatmap(
  data = data_inter,
  tax_aggregate = "Class",
  tax_add = "Phylum",
  tax_show = 20,
  showRemainingTaxa = TRUE,
  order_y_by = rev(sort(unique(order_by_y_vec))),
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  color_vector = c( "whitesmoke", "royalblue4")
)

#----------------------------------------------------------------------------------------------------------------------------------
#BIOCRUST CORE HEATMAP

#keep only rows that are "Core"

df_biocrust2 <- df_biocrust[df_biocrust$Shared == "Core", ]

#remove last column (using dplyr)

df_biocrust2 <- select(df_biocrust2, -Shared)

#rename row headers (using dplyr)

df_biocrust2 <- df_biocrust2 %>%
  rename(
    OTU2 = OTU,
  )

# Match up values from asv.data to df_biocrust2

df_biocrust3 =
  asv.data %>%
  filter(`OTU` %in% df_biocrust2$OTU2) %>%
  select(c(`OTU`, names(df_biocrust2) %>% grep("Juliette", ., value = T))) %>%
  right_join(df_biocrust2[,1:8], by = c("OTU" = "OTU2"))


data_biocrust <- amp_load(
  otutable = df_biocrust3,
   metadata = "D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/all/OUT_300721/halden_mapping_ampvis.txt",
  taxonomy = NULL,
  fasta = NULL,
  tree = NULL,
  pruneSingletons = TRUE
)

order_by_y_vec = paste(data_biocrust$tax$Phylum, data_biocrust$tax$Class, sep = "; ")


plot_biocrust <- amp_heatmap(
  data = data_biocrust,
  tax_aggregate = "Class",
  tax_add = "Phylum",
  tax_show = 20,
  showRemainingTaxa = TRUE,
  order_y_by = rev(sort(unique(order_by_y_vec))),
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  color_vector = c( "whitesmoke", "royalblue4")
)


#add them all up into one figure

library(ggpubr)
core_all <- ggarrange(plot_heap, plot_inter, plot_biocrust,
                         labels = c("A.",  "B.", "C."),
                         label.y = 0.05,
                         label.x = 0.05,
                         font.label = list(size = 18),
                         nrow = 1,
                         ncol = 3,
                         #widths = c(1,0.05,1,0.05,1, 0.05,1,0.05,1 ),
                         common.legend = TRUE, legend = "right"
                      )
core_all



#----------------------------------------------------------------------------------------------------------------------------------
# HEAP CORE HEATMAP - GENUS LEVEL CLASSIFICATION

#order by phylum

order_by_y_vec = paste(data_heap$tax$Phylum, data_heap$tax$Genus, sep = "; ") 

g1 <- amp_heatmap(
  data = data_heap,
  tax_aggregate = "Genus",
  tax_add = "Phylum",
  tax_show = 20,
  showRemainingTaxa = TRUE,
  #order_y_by = rev(unique(order_by_y_vec)), #call previous variable made 
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  color_vector = c( "whitesmoke", "royalblue4")
)

#----------------------------------------------------------------------------------------------------------------------------------
# INTER CORE HEATMAP - GENUS LEVEL CLASSIFICATION

#order by phylum

order_by_y_vec = paste(data_inter$tax$Phylum, data_inter$tax$Genus, sep = "; ")

g2 <- amp_heatmap(
  data = data_inter,
  tax_aggregate = "Genus",
  tax_add = "Phylum",
  tax_show = 20,
  showRemainingTaxa = TRUE,
  order_y_by = rev(sort(unique(order_by_y_vec))),
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  color_vector = c( "whitesmoke", "royalblue4")
)

#----------------------------------------------------------------------------------------------------------------------------------
# BIOCRUST CORE HEATMAP - GENUS LEVEL CLASSIFICATION

#order by phylum
order_by_y_vec = paste(data_biocrust$tax$Phylum, data_biocrust$tax$Genus, sep = "; ")


g3 <-amp_heatmap(
  data = data_biocrust,
  tax_aggregate = "Genus",
  tax_add = "Phylum",
  tax_show = 20,
  showRemainingTaxa = TRUE,
  order_y_by = rev(sort(unique(order_by_y_vec))),
  tax_class = TRUE,
  plot_values = TRUE, #this shows the numerical value of reads
  group_by = "site",
  facet_by = "type",
  plot_colorscale = "sqrt",
  color_vector = c( "whitesmoke", "royalblue4")
)

core_all_genus <- ggarrange(g1, g2, g3,
                         labels = c("A.",  "B.", "C."),
                         label.y = 0.05,
                         label.x = 0.05,
                         font.label = list(size = 18),
                         nrow = 1,
                         ncol = 3,
                         #widths = c(1,0.05,1,0.05,1, 0.05,1,0.05,1 ),
                         common.legend = TRUE, legend = "right"
                      )
core_all_genus


```

# ABIOTIC VALUES BAR PLOT

```{r}

library(ggplot2)

#load in abiotic dataset 

abiotic_wide1 <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Abiotic_Data/20221123_abiotic_wide.csv", 
                  header = TRUE, 
           sep = ",",
           dec = ".")

#reorder

abiotic_wide1$stage <- factor(abiotic_wide1$stage, levels =c("Heap","Initial", "Biocrust"))

#barplot for pH

p1 <- 
  ggplot(data = abiotic_wide_BSC, aes(x = stage, y = pH_avg, fill = site)) + 
  geom_bar(stat="identity", color = "black", position=position_dodge()) + 
  # keep only upper error bars
  geom_errorbar(aes(ymin = pH_avg, ymax = pH_avg+pH_sd), width = 0.2, 
  position=position_dodge(.9)) +
  theme_classic() +
 # facet_wrap(~ type, scales = "free", ncol = 5, ) +
  theme(
    plot.title = element_text(hjust = 0.5),
    #remove legend
  legend.position = "none") + 
  labs(title = "pH", x = "", y = "") +
    scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black")) +
   # add spacing 
  theme(panel.spacing = unit(2, "lines")) +
  #change text size
  theme(plot.title = element_text(size = 20)) +
  theme(axis.title = element_text(size = 20))  +
  #change axis labels
  theme(axis.text.y = element_text(size = 20)) +
  theme(axis.text.x = element_text(size = 15)) +
  #only change legend title
  theme(legend.title = element_text(size = 20)) +
  #only change legend text
  theme(legend.text = element_text(size = 20)) 
  

#barplot for EC

p2 <- 
  ggplot(data = abiotic_wide_BSC, aes(x = stage, y = EC_avg, fill = site)) + 
  geom_bar(stat="identity", color = "black", position=position_dodge()) + 
  # keep only upper error bars
  geom_errorbar(aes(ymin = EC_avg, ymax = EC_avg+EC_sd), width = 0.2, 
  position=position_dodge(.9)) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5),
        #remove legend
  legend.position = "none") +  
  
  labs(title = "EC", x = "", y = "mS/cm") +
    scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black"))+
    # add spacing 
  theme(panel.spacing = unit(2, "lines")) +
  #change text size
  theme(plot.title = element_text(size = 20)) +
  theme(axis.title = element_text(size = 20))  +
  #change axis labels
  theme(axis.text.y = element_text(size = 20)) +
  theme(axis.text.x = element_text(size = 15)) +
  #only change legend title
  theme(legend.title = element_text(size = 20)) +
  #only change legend text
  theme(legend.text = element_text(size = 20)) 



#barplot for Chlorophyll a

p3 <- ggplot(data = abiotic_wide_BSC, aes(x = stage, y = Chla_avg, fill = site)) + 
  geom_bar(stat="identity", color = "black", position=position_dodge()) + 
  # keep only upper error bars
  geom_errorbar(aes(ymin = Chla_avg, ymax = Chla_avg+Chla_sd), width = 0.2, 
  position=position_dodge(.9)) +
  theme_classic() +
 # facet_wrap(~ type, scales = "free", ncol = 5, ) +
  theme(
    plot.title = element_text(hjust = 0.5),
        #remove legend
  legend.position = "none") +  
  labs(title = "Chl a", x = "", y = "μg/g dry soil") +
    scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black"))+
  # add spacing 
  theme(panel.spacing = unit(2, "lines")) +
  #change text size
  theme(plot.title = element_text(size = 20)) +
  theme(axis.title = element_text(size = 20))  +
  #change axis labels
  theme(axis.text.y = element_text(size = 20)) +
  theme(axis.text.x = element_text(size = 15)) +
  #only change legend title
  theme(legend.title = element_text(size = 20)) +
  #only change legend text
  theme(legend.text = element_text(size = 20)) 


#barplot for DOC

p4 <- ggplot(data = abiotic_wide_BSC, aes(x = stage, y = DOC_avg, fill = site)) + 
  geom_bar(stat="identity", color = "black", position=position_dodge()) + 
  # keep only upper error bars
  geom_errorbar(aes(ymin = DOC_avg, ymax = DOC_avg+DOC_sd), width = 0.2, 
  position=position_dodge(.9)) +
  theme_classic() +
 # facet_wrap(~ type, scales = "free", ncol = 5, ) +
  theme(
    plot.title = element_text(hjust = 0.5),
        #remove legend
  legend.position = "none") +  
  labs(title = "DOC", x = "", y = "μg/g dry soil") +
    scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black"))+
  # add spacing 
  theme(panel.spacing = unit(2, "lines")) +
  #change text size
  theme(plot.title = element_text(size = 20)) +
  theme(axis.title = element_text(size = 20))  +
  #change axis labels
  theme(axis.text.y = element_text(size = 20)) +
  theme(axis.text.x = element_text(size = 15)) +
  #only change legend title
  theme(legend.title = element_text(size = 20)) +
  #only change legend text
  theme(legend.text = element_text(size = 20)) 


#barplot for TDN

p5 <- ggplot(data = abiotic_wide_BSC, aes(x = stage, y = TNb_avg, fill = site)) + 
  geom_bar(stat="identity", color = "black", position=position_dodge()) + 
  # keep only upper error bars
  geom_errorbar(aes(ymin = TNb_avg, ymax = TNb_avg+TNb_sd), width = 0.2, 
  position=position_dodge(.9)) +
  theme_classic() +
 # facet_wrap(~ type, scales = "free", ncol = 5, ) +
  theme(
    plot.title = element_text(hjust = 0.5),
    #remove legend
  legend.position = "none") +   
  labs(title = "TDN", x = "", y = "μg/g dry soil") +
    scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black"))+
  # add spacing 
  theme(panel.spacing = unit(2, "lines")) +
  #change text size
  theme(plot.title = element_text(size = 20)) +
  theme(axis.title = element_text(size = 20))  +
  #change axis labels
  theme(axis.text.y = element_text(size = 20)) +
  theme(axis.text.x = element_text(size = 15)) +
  #only change legend title
  theme(legend.title = element_text(size = 20)) +
  #only change legend text
  theme(legend.text = element_text(size = 20)) 

#combine plots in 1 row

library(ggpubr)

abiotic_BSC <- ggarrange(p1, NULL, p2, NULL, p3, NULL, p4, NULL, p5, #CCA_plot,
          labels = c("A.", "", "B.", "","C.", "","D.", "","E."),
          label.y = 0.1,
          label.x = 0.1,
          font.label = list(size = 18),
             nrow = 1, widths = c(1,0.05,1,0.05,1,0.05,1,0.05,1 ),
          common.legend = TRUE, legend = "right")
abiotic_BSC

#combine plots in 2 rows 

abiotic_BSC <- ggarrange(p1, p2, NULL, p3,  p4,  p5, #CCA_plot,
                         labels = c("B.",  "C.", "", "D.", "E.", "F."),
                         label.y = 0.1,
                         label.x = 0.1,
                         font.label = list(size = 18),
                         nrow = 2, 
                         ncol = 3,
                         #widths = c(1,0.05,1,0.05,1, 0.05,1,0.05,1 ),
                         common.legend = TRUE, legend = "right")
abiotic_BSC

#TOTAL CARBON

ggplot(abiotic_wide1, aes(x = stage, y = pH, fill = site)) +
 geom_boxplot() + facet_wrap(~type) + 
  scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black")) + 
    #change labels
 labs(title = "pH", x = "Stage", y = "NA") +
 #theme_minimal()+
   theme_bw()+
  #change theme
  theme(
    #add borders
  panel.border = element_rect(color = "grey85", fill = NA),
    #element_rect(fill = "black"),
  #legend
  legend.position = "right",
  #center title
  plot.title = element_text(hjust = 0.5),
  #change all text size
  text = element_text(size = 20)
         ) 

```

# STATISTICAL ANALYSIS - ABIOTIC AND QPCR DATA

```{r}

library("dplyr")
library("ggpubr")
library("MASS")
library(car)

qPCR_wide <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Raw_Data/qPCR/qPCR_wide.csv", 
                  header = TRUE, 
           sep = ",",
           dec = ".")

abiotic_wide <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Abiotic_Data/20221122_abiotic_wide.csv", 
                  header = TRUE, 
           sep = ",",
           dec = ".")

#LINEAR MODELS
  
#chlorophyll

log_Chl_a <- log(abiotic_wide_BSC$Chla)
shapiro.test(log_Chl_a)
chla_test <- lm(log_Chl_a ~ site*stage, data = abiotic_wide_BSC)
summary(chla_test)
plot(lm(log_Chl_a ~ site*stage, data = abiotic_wide_BSC))
  

#pairwise comparisons  

lsmeans(chla_test, pairwise~site*stage, adjust="tukey")

lsmeans(chla_test, pairwise~site*stage, adjust="none")
marginal = lsmeans(chla_test, ~ site*stage)
multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")
  

#qPCR: 16S bacteria -----------

#transforming 16S bacteria

log_bac_16S <- log(qPCR_wide$bac_16S)
shapiro.test(log_bac_16S) 

#linear model

test1 <- lm(qPCR_wide$bac_16S ~ site*type, data = qPCR_wide)
summary(test1) 
plot(test1) 

test2 <- lm(log_bac_16S ~ site*type, data = qPCR_wide)
summary(test2) 
plot(test2)

#generalized linear model

test2g <- glm(log_bac_16S ~ site*type, data = qPCR_wide)
test2g
plot(test2g)

#calculate residuals

res1 <- residuals(test1) 
res2 <- residuals(test2)

#shapiro test for normality for residuals 
shapiro.test(res1) 
shapiro.test(res2) 

AIC(test1, test2) 

#ANOVA

test_a1 <- aov(log_bac_16S ~ site*type, data = qPCR_wide)

summary(test_a1)

#pairwise comparisons 

lsmeans(test2, pairwise~site*type, adjust="none")

marginal = lsmeans(test2, ~ site*type)

multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")
    
#qPCR: 16S bacteria per g DRY SOIL

#transforming 16S bacteria

log_bac_16S_ds <- log(qPCR_wide$bac_16S_ds)
shapiro.test(log_bac_16S_ds)

#linear model

test1 <- lm(qPCR_wide$bac_16S_ds ~ site*type, data = qPCR_wide)
summary(test1) 
plot(test1)

test2 <- lm(log_bac_16S_ds ~ site*type, data = qPCR_wide)
summary(test2)
plot(test2) 

#generalized linear model

test2g <- glm(log_bac_16S_ds ~ site*type, data = qPCR_wide)
test2g
plot(test2g)

#calculate residuals

res1 <- residuals(test1) 
res2 <- residuals(test2)

#shapiro test for normality for residuals

shapiro.test(res1)
shapiro.test(res2)


AIC(test1, test2) 

#ANOVA 

test_a1 <- aov(log_bac_16S_ds ~ site*type, data = qPCR_wide)
summary(test_a1)

#pairwise comparisons 

lsmeans(test2, pairwise~site*type, adjust="none")
marginal = lsmeans(test2, ~ site*type)
multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")  
  
#qPCR: 16S archaea per g DRY SOIL  

#remove 0 values   

arc_16S_data_ds <- qPCR_wide[-9,]

#transform 

log_arc_16S_ds <- log(arc_16S_data_ds$arc_16S_ds)
shapiro.test(log_arc_16S_ds)
ggqqplot(log_arc_16S_ds) 

test3 <- lm(arc_16S_data_ds$arc_16S_ds ~ site*type, data = arc_16S_data_ds)
summary(test3) 
res3 <- residuals(test3)
shapiro.test(res3)

#try a different transformation

spreadLevelPlot(test3) 

test4 <- lm(log_arc_16S_ds ~ site*type, data = arc_16S_data_ds)
summary(test4) 
res4 <- residuals(test4)
shapiro.test(res4)

#plot

plot(test3) 
plot(test4) 

AIC(test3, test4) 

#pairwise comparisons

lsmeans(test4, pairwise~site*type, adjust="none")
marginal = lsmeans(test4, ~ site*type)
multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")

  
  
  
#qPCR Archaeal 16S

#remove 0 values  

arc_16S_data <- qPCR_wide[-9,]
#transform 

log_arc_16S <- log(arc_16S_data$arc_16S)
shapiro.test(log_arc_16S) 
ggqqplot(log_arc_16S)

test3 <- lm(qPCR_wide$arc_16S ~ site*type, data = qPCR_wide)

summary(test3)
res3 <- residuals(test3)

shapiro.test(res3)

#try a different transformation

spreadLevelPlot(test3)

test4 <- lm(log_arc_16S ~ site*type, data = arc_16S_data)
summary(test4) 
res4 <- residuals(test4)
shapiro.test(res4) 

#plot
plot(test3)
plot(test4)

AIC(test3, test4)

#pairwise comparisons

lsmeans(test4, pairwise~site*type, adjust="none")
marginal = lsmeans(test4, ~ site*type)
multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")

#DOC

shapiro.test(abiotic_wide_BSC$DOC)
log_DOC <- log(abiotic_wide_BSC$DOC)
shapiro.test(log_DOC)
test5 <- (lm(log_DOC ~ site*stage, data = abiotic_wide_BSC))

#test residuals

res5 <- residuals(test5)
shapiro.test(res5)
plot(test5)
summary(test5)

#pairwise comparisons 

lsmeans(test5, pairwise~site*stage, adjust="none")
marginal = lsmeans(test5, ~ site*stage)
multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")

#pH
shapiro.test(abiotic_wide_BSC$pH)
plot(lm(abiotic_wide_BSC$pH ~ site*stage, data = abiotic_wide_BSC))
ph <- lm(abiotic_wide_BSC$pH ~ site*stage, data = abiotic_wide_BSC)
summary(ph)

#pairwise comparisons 

lsmeans(ph, pairwise~site*stage, adjust="none")
marginal = lsmeans(ph, ~ site*stage)
multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")

#TDN

log_TDN <- log(abiotic_wide_BSC$TNb)
shapiro.test(log_TDN)
test6 <- (lm(log_TDN ~ site*stage, data = abiotic_wide_BSC))

#test residuals

res6 <- residuals(test6)
shapiro.test(res6)
plot(test6)
summary(test6)

#pairwise comparisons 

lsmeans(test6, pairwise~site*stage, adjust="none")
marginal = lsmeans(test6, ~ site*stage)
multcomp::cld(marginal,
              alpha=0.05,
              Letters=letters, adjust="none")


#EC

ec_set <- lm(EC ~ site*stage, data = abiotic_wide_BSC)
shapiro.test(abiotic_wide_BSC$EC) 
summary(ec_set)

#test residuals

resec_set <- residuals(ec_set)
shapiro.test(resec_set) 
plot(resec_set)

#log xform

log_EC <- log(abiotic_wide_BSC$EC)
test7 <- (lm(log_EC ~ site*stage, data = abiotic_wide_BSC))
res7 <- residuals(test7)
shapiro.test(res7)
plot(test7)


ec_set <- lm(formula = EC ~ site+stage, data = abiotic_wide_BSC)

#create vector of VIF values

vif_values <- vif(ec_set)

#create horizontal bar chart to display each VIF value

barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "steelblue")

#add vertical line at 5

abline(v = 5, lwd = 3, lty = 2)

ggqqplot(abiotic_wide_BSC$EC)

box_cox <- boxcox(log_EC ~ site*stage, data = abiotic_wide_BSC)

(lambda <- box_cox$x[which.max(box_cox$y)])  

new_model <- lm(((log_EC^lambda-1)/lambda) ~ site*stage, data = abiotic_wide_BSC)
qqnorm(new_model$residuals)
plot(new_model)
res8 <- residuals(new_model)
shapiro.test(res8) 6

#check for outliers

ggplot(abiotic_wide_BSC) +
  aes(x = , y = EC) +
  geom_boxplot(fill = "#0c4c8a") +
  theme_minimal()

#try GLM

test3g <- glm(log_EC ~ site*stage, data = abiotic_wide_BSC)
test3g 
res9 <- residuals(test3g)
shapiro.test(res9)

#try a different boxcox

box_cox2 <- log_EC^coef(powerTransform(log_EC))
test8 <- lm(box_cox2~site*stage, data = abiotic_wide_BSC)
res10 <- residuals(test8)
shapiro.test(res10)

#try a different transformation

spreadLevelPlot(lm(EC ~ site*stage, data = abiotic_wide_BSC)) 
pow_EC <- (abiotic_wide_BSC$EC^-0.7790933)
test11 <- (lm(pow_EC ~ site*stage, data = abiotic_wide_BSC))
res11 <- residuals(test11)
shapiro.test(res11)
plot(test11)
summary(test11) 

#non-constant error variance test (Breusch-Pagan test)
ncvTest(test8)

#now do non-parametric tests for this
kruskal.test(EC ~ stage, data = abiotic_wide_BSC) 

#ABIOTIC PARAMETERS BY SITE AND TYPE (site: WT/OD, )
wilcox.test(EC ~ site, data = abiotic_wide_BSC) 

pairwise.wilcox.test(abiotic_wide_BSC$EC, abiotic_wide_BSC$stage,
                     p.adjust.method = "BH")

```

# qPCR BAR PLOT -

```{r dev = c("png", "jpg", "pdf")}

library(ggplot2)
library(scales)

#load data

qPCR_long <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Raw_Data/qPCR/qPCR_long.csv", 
                  header = TRUE, 
           sep = ",",
           dec = ".")
str(qPCR_long)

#add pseudocount to whole matrix before log transformation

qPCR_long[,4] <- qPCR_long[,4] + 1

#reorder genes

qPCR_long$gene <- factor(qPCR_long$gene, levels = c("bac16S" ,"arc16S", "ds.bac16S" ,"ds.arc16S"))

#reorder stages "type"

qPCR_long$type <- factor(qPCR_long$type, levels = c("Heap" ,"Initial" ,"Biocrust"))

qPCR_long2 <- subset(qPCR_long,
                     qPCR_long$gene != "ds.bac16S" & 
                     qPCR_long$gene != "ds.arc16S"
                    # qPCR_long$gene != "bac16S" & 
                    # qPCR_long$gene != "arc16S" &
                     )

#rename headers

qPCR_long2$gene <- factor(qPCR_long2$gene, levels = c("bac16S","arc16S" ),
                     labels = c ("Bacterial 16S", "Archaeal 16S"))

#plot copies/ng DNA for all bac, arc

q1 <- ggplot(qPCR_long2, aes(x=type, y=copies, fill = site)) +
  geom_bar(stat="identity", color="black", 
           position=position_dodge()) +
  scale_y_log10(labels = scales::scientific) + #log scale, scientific notation
  geom_errorbar(aes(ymin = copies, ymax = copies+sd), width = 0.2, #keep only upper error bars
  position=position_dodge(.9)) +
  labs(x = "", y = "copies ng-1 DNA") + 
  scale_x_discrete(labels = new_type) +
  theme_classic() + 
 #increase text size
  theme(text = element_text(size = 20)) +
  facet_wrap(~ gene) +
  #change barplot colors
  scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black")) 

#plot copies/ g dry soil for bac, arc

qPCR_long3 <- subset(qPCR_long,
                    # qPCR_long$gene != "ds.bac16S" & 
                    # qPCR_long$gene != "ds.arc16S" &
                     qPCR_long$gene != "bac16S" & 
                     qPCR_long$gene != "arc16S"
                     )
#rename headers

qPCR_long3$gene <- factor(qPCR_long3$gene, levels = c("ds.bac16S","ds.arc16S" ),
                     labels = c ("Bacterial 16S", "Archaeal 16S"))

 q2 <- ggplot(qPCR_long3, aes(x=type, y=copies, fill = site)) +
  geom_bar(stat="identity", color="black", 
           position=position_dodge()) +
  scale_y_log10(labels = scales::scientific) + #log scale, scientific notation
  geom_errorbar(aes(ymin = copies, ymax = copies+sd), width = 0.2, #keep only upper error bars
  position=position_dodge(.9)) +
  labs(x = "", y = "copies g-1 dry soil") + 
  theme_classic() + 
   #increase text size
  theme(text = element_text(size = 20)) +
   #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10) +
  #scale_y_continuous(labels = scales::scientific) + #scientific notation
  #facet_wrap(~ gene, scales = "free")
  scale_x_discrete(labels = new_type) +
  facet_wrap(~ gene) +
  #change barplot colors
  scale_fill_manual(values = c("OD" = "lightgrey",
                               "WT" = "black")) 
 
#assemble plots together
 
library(ggpubr)
 
qPCR_plots <- ggarrange(q1, q2, 
          labels = c("A.", "B."),
          label.y = 0.07,
          label.x = 0.05,
          nrow = 2,
          font.label = list(size = 18),
         common.legend = TRUE, legend = "right")

qPCR_plots 
 
```

# BETA NTI STEGEN 2013

```{r}

#load package for Beta NTI

library("picante") #version 1.8.2

#read in otu table from phyloseq object 

otu = otu_table(phytted)

## read in the phylogeny

phylo = phy_tree(phytted)

phylo


## make sure the names on the phylogeny are ordered the same as the names in otu table

match.phylo.otu = match.phylo.data(phylo, otu)

str(match.phylo.otu)

## calculate empirical betaMNTD

beta.mntd.weighted = as.matrix(comdistnt(t(match.phylo.otu$data),cophenetic(match.phylo.otu$phy),abundance.weighted=T));
dim(beta.mntd.weighted);
beta.mntd.weighted[1:5,1:5];
write.csv(beta.mntd.weighted,'betaMNTD_weighted.csv',quote=F);

identical(colnames(match.phylo.otu$data),colnames(beta.mntd.weighted)); 
identical(colnames(match.phylo.otu$data),rownames(beta.mntd.weighted)); 

# calculate randomized betaMNTD

beta.reps = 999;

rand.weighted.bMNTD.comp = array(c(-999),dim=c(ncol(match.phylo.otu$data),ncol(match.phylo.otu$data),beta.reps));
dim(rand.weighted.bMNTD.comp);


for (rep in 1:beta.reps) {
  
  rand.weighted.bMNTD.comp[,,rep] = as.matrix(comdistnt(t(match.phylo.otu$data),taxaShuffle(cophenetic(match.phylo.otu$phy)),abundance.weighted=T,exclude.conspecifics = F));
  
  print(c(date(),rep));
  
}

weighted.bNTI = matrix(c(NA),nrow=ncol(match.phylo.otu$data),ncol=ncol(match.phylo.otu$data));
dim(weighted.bNTI);

for (columns in 1:(ncol(match.phylo.otu$data)-1)) {
  for (rows in (columns+1):ncol(match.phylo.otu$data)) {
    
    rand.vals = rand.weighted.bMNTD.comp[rows,columns,];
    weighted.bNTI[rows,columns] = (beta.mntd.weighted[rows,columns] - mean(rand.vals)) / sd(rand.vals);
    rm("rand.vals");
    
  };
};

rownames(weighted.bNTI) = colnames(match.phylo.otu$data);
colnames(weighted.bNTI) = colnames(match.phylo.otu$data);
weighted.bNTI;
write.csv(weighted.bNTI,"weighted_bNTI.csv",quote=F);

pdf("weighted_bNTI_Histogram.pdf")
  hist(weighted.bNTI)
dev.off()

```

# "Raup-Crick Abundance" STEGEN 2013

```{r}

#spXsite = as.data.frame(otu)

spXsite = as.data.frame(t(otu))

raup_crick_abundance = function(spXsite, plot_names_in_col1=FALSE, classic_metric=FALSE, split_ties=TRUE, reps=999, set_all_species_equal=FALSE, as.distance.matrix=TRUE, report_similarity=FALSE){
	
	
	##this section moves plot names in column 1 (if specified as being present) into the row names of the matrix and drops the column of names
	if(plot_names_in_col1){
		row.names(spXsite)<-spXsite[,1]
		spXsite<-spXsite[,-1]
		}
	
	
	## count number of sites and total species richness across all plots (gamma)
	n_sites<-nrow(spXsite)
	gamma<-ncol(spXsite)

	##build a site by site matrix for the results, with the names of the sites in the row and col names:
	results<-matrix(data=NA, nrow=n_sites, ncol=n_sites, dimnames=list(row.names(spXsite), row.names(spXsite)))
	
	##make the spXsite matrix into a new, pres/abs. matrix:
	ceiling(spXsite/max(spXsite))->spXsite.inc
	
	##create an occurrence vector- used to give more weight to widely distributed species in the null model:
	occur<-apply(spXsite.inc, MARGIN=2, FUN=sum)
	
	##create an abundance vector- used to give more weight to abundant species in the second step of the null model:
	abundance<-apply(spXsite, MARGIN=2, FUN=sum)
	
	##make_null:
	##looping over each pairwise community combination:
	
	for(null.one in 1:(nrow(spXsite)-1)){
		for(null.two in (null.one+1):nrow(spXsite)){
			
			null_bray_curtis<-NULL
			for(i in 1:reps){
				
				##two empty null communities of size gamma:
				com1<-rep(0,gamma)
				com2<-rep(0,gamma)
				
				##add observed number of species to com1, weighting by species occurrence frequencies:
				com1[sample(1:gamma, sum(spXsite.inc[null.one,]), replace=FALSE, prob=occur)]<-1
				com1.samp.sp = sample(which(com1>0),(sum(spXsite[null.one,])-sum(com1)),replace=TRUE,prob=abundance[which(com1>0)]);
				com1.samp.sp = cbind(com1.samp.sp,1); # head(com1.samp.sp);
				com1.sp.counts = as.data.frame(tapply(com1.samp.sp[,2],com1.samp.sp[,1],FUN=sum)); colnames(com1.sp.counts) = 'counts'; # head(com1.sp.counts);
				com1.sp.counts$sp = as.numeric(rownames(com1.sp.counts)); # head(com1.sp.counts);
				com1[com1.sp.counts$sp] = com1[com1.sp.counts$sp] + com1.sp.counts$counts; # com1;
				#sum(com1) - sum(spXsite[null.one,]); ## this should be zero if everything work properly
				rm('com1.samp.sp','com1.sp.counts');			
				
				##same for com2:
				com2[sample(1:gamma, sum(spXsite.inc[null.two,]), replace=FALSE, prob=occur)]<-1
				com2.samp.sp = sample(which(com2>0),(sum(spXsite[null.two,])-sum(com2)),replace=TRUE,prob=abundance[which(com2>0)]);
				com2.samp.sp = cbind(com2.samp.sp,1); # head(com2.samp.sp);
				com2.sp.counts = as.data.frame(tapply(com2.samp.sp[,2],com2.samp.sp[,1],FUN=sum)); colnames(com2.sp.counts) = 'counts'; # head(com2.sp.counts);
				com2.sp.counts$sp = as.numeric(rownames(com2.sp.counts)); # head(com2.sp.counts);
				com2[com2.sp.counts$sp] = com2[com2.sp.counts$sp] + com2.sp.counts$counts; # com2;
				# sum(com2) - sum(spXsite[null.two,]); ## this should be zero if everything work properly
				rm('com2.samp.sp','com2.sp.counts');

				null.spXsite = rbind(com1,com2); # null.spXsite;
				
				##calculate null bray curtis
				null_bray_curtis[i] = vegdist(null.spXsite,method='bray');
				
			}; # end reps loop

			## empirically observed bray curtis
			obs.bray = vegdist(spXsite[c(null.one,null.two),],method='bray');

			##how many null observations is the observed value tied with?
			num_exact_matching_in_null = sum(null_bray_curtis==obs.bray);
			
			##how many null values are smaller than the observed *dissimilarity*?
			num_less_than_in_null = sum(null_bray_curtis<obs.bray);
			
			rc = (num_less_than_in_null )/reps; # rc;
			
			if(split_ties){
				
				rc = ((num_less_than_in_null +(num_exact_matching_in_null)/2)/reps)
			};
			
			
			if(!classic_metric){
					
					##our modification of raup crick standardizes the metric to range from -1 to 1 instead of 0 to 1
					
					rc = (rc-.5)*2
			};

			results[null.two,null.one] = round(rc,digits=2); ##store the metric in the results matrix
			
			print(c(null.one,null.two,date()));
			
		}; ## end null.two loop
		
	}; ## end null.one loop
	
	if(as.distance.matrix){ ## return as distance matrix if so desired
		results<-as.dist(results)
	}	
	
return(results)

}; ## end function


#then run the function

raup_crick_abundance(spXsite)

```


# plotting figures for Community Assembly analysis and Beta-NTI

```{r}

#Make a percent stacked barplot

library(ggplot2)

#load dataset

com_ass_df <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/R/betaNTI/com_ass.csv", 
                  header = TRUE, 
           sep = ",",
           dec = ".")

#order the site names

com_ass_df$label <- factor(com_ass_df$label, levels =c("OD heap","OD initial", "OD biocrust", "WT heap", "WT initial", "WT biocrust"))

#order the Community Processes

com_ass_df$process <- factor(com_ass_df$process, levels = c("Undominated", "Homogenizing dispersal","Dispersal limitation", "Homogenizing selection", "Variable selection"))

#then plot

comm_ass_plot <- ggplot(com_ass_df, aes(x = label, y = perc_site_pairs, fill = process)) + 
  geom_bar(stat = "identity", color = "black") +
  #change colors 
  scale_fill_manual(values = c("Undominated" = "snow",
                               "Homogenizing dispersal" = "gray80",
                               "Dispersal limitation" = "gray55",
                               "Homogenizing selection" = "sandybrown",
                               "Variable selection" = "palegreen4")) +
  #axis labels
  xlab ("site") + ylab ("Percent of site pairs") +
  #legend customization
  guides(fill = guide_legend(title = "Assembly Process")
  )

#make the text size good

comm_ass_plot +   
  #increase font size overall
  #theme(text = element_text(size = 20))
  #only change axis titles 
  theme(axis.title = element_text(size = 20))  +
  #change axis labels
  theme(axis.text.y = element_text(size = 20)) +
  theme(axis.text.x = element_text(size = 15)) +
  #only change legend title
  theme(legend.title = element_text(size = 20)) +
  #only change legend text
  theme(legend.text = element_text(size = 20))   


#make a box and whisker plot for beta NTI 

#load dataset

beta_NTI_df <- read.csv("D:/Users/juliette.ohan/Desktop/2_Saltheaps_(Ohan)/Experiments/Scripts_(R_and_Bioinformatics)/R/betaNTI/b_NTI_all.csv", 
                        header = TRUE, 
                        sep = ",",
                        dec = ".")

#order the site names

beta_NTI_df$label <- factor(beta_NTI_df$label, levels =c("OD_heap","OD_initial", "OD_biocrust", "WT_heap", "WT_initial", "WT_biocrust"))


#plot

bNTI_plot <- ggplot(beta_NTI_df, aes(x=label, y=bNTI)) + 
  geom_boxplot() +
  #add axis labels
  xlab ("site") + ylab ("βNTI") +
    #add points
   #geom_jitter(color="black", size=0.4, alpha=0.9) +
  #add reference lines
  geom_hline(yintercept = 2, linetype = 2) +
  geom_hline(yintercept = -2, linetype = 2) 

bNTI_plot +   
  #increase font size overall
  #theme(text = element_text(size = 20))
  #only change axis titles 
  theme(axis.title = element_text(size = 20))  +
  #change axis labels
  theme(axis.text.y = element_text(size = 20)) +
  theme(axis.text.x = element_text(size = 15)) +
  #only change legend title
  theme(legend.title = element_text(size = 20)) +
  #only change legend text
  theme(legend.text = element_text(size = 20))  

```