#Analysing the diet of Alaskan Red knot chicks using meta-barcoding (2019-2024)
#Script by Roos Winters
#In this script, the asv-file is transformed/prepped for data  analysis. 

#####Start Script####
remove(list=ls()) # clear everything in memory

#Load required libraries
library(tidyverse)
library(cowplot)
library(lme4)
library(qdapTools)
library(writexl)
library(readxl)
library(devtools)
library(vegan) 
library(openxlsx)
library(ggpubr)

#Set working directory 
setwd("~/Alaskan chick diet - paper/Final_figures/" )


#load ASV file 
seq_data <- read_xlsx("./Roos_data_and_357_358.xlsx")

all_samples <- read_xlsx("./Mapping_file_all.xlsx") %>% 
  rename(sample = Seq_sample_name) %>% 
  select(sample)

#####1. add missing samples to the dataset #####
samples_in_sequencedata <- seq_data %>% 
  pivot_longer(cols = -c(`#OTU_ID`, taxonomy),names_to = "sample", values_to = "reads") %>% 
  select(sample) %>% 
  unique() 

missing_samples_vec <- setdiff(all_samples$sample, samples_in_sequencedata$sample)

#add them as 0-columns 
seq_data <- seq_data %>%
  bind_cols(as_tibble(setNames(
    replicate(length(missing_samples_vec), 0, simplify = FALSE),
    missing_samples_vec
  ))) %>%
  relocate(all_of(missing_samples_vec), .before = c(taxonomy, identity))

#####1. Negative controls#####
#Check the negative extracts for possible contamination.  
#Do not forget to add the NTC's

neg_controls <- read_excel("./Negative_control_list.xlsx") %>% 
  pull(Negative_Controls)

# Filter + calculate sums for those controls
neg_control <- seq_data %>% 
  select(`#OTU_ID`, all_of(neg_controls), taxonomy) %>% 
  rowwise() %>% 
  filter(sum(c_across(-c(`#OTU_ID`, taxonomy))) > 0) %>%   # keep only rows with reads
  mutate(sum_control_reads = sum(c_across(-c(`#OTU_ID`, taxonomy)))) %>% 
  ungroup()

neg_control$`#OTU_ID` #Gives a list of all contaminants 

total_reads <- seq_data %>% 
  rowwise() %>% 
  mutate(sum_reads = sum(c_across(-c("#OTU_ID", "taxonomy", "identity")))) %>% 
  select(c("#OTU_ID", taxonomy, sum_reads))

reads_comparison <- left_join(total_reads, neg_control, by="#OTU_ID")  %>%
  select(c("#OTU_ID", taxonomy.x, sum_reads, sum_control_reads)) %>%
  rename(taxonomy = taxonomy.x) %>%
  filter(!is.na(sum_control_reads)) %>% 
  mutate(difference = sum_reads-sum_control_reads)

contaminants <- reads_comparison %>% 
  filter(difference == 0) %>% 
  pull("#OTU_ID")

seq_data <- seq_data %>% filter(!`#OTU_ID` %in% contaminants)

#Many OTU_ID's that occur in the negative controls do not occur in any other sample (difference = 0). These can be removed from the dataset without further investigation. This leaves 28 / 89 OTU's 

reads_comparison_2 <- reads_comparison %>% 
  filter(!difference == 0)

#Let's transport an overview of the negative controls to excel 
#Nice if it's one document: use workbooks
#remove double ## to save

##wb <- createWorkbook()
##addWorksheet(wb, "Negative_controls")
##writeData(wb, "Negative_controls", neg_control)
##addWorksheet(wb, "Reads_comparison")
##writeData(wb, "Reads_comparison", reads_comparison)
##addWorksheet(wb, "Reads_comparison_2")
##writeData(wb, "Reads_comparison_2", reads_comparison_2)

# save workbook
##saveWorkbook(wb, "./overview_negative_controls_all.xlsx", overwrite = TRUE)

#Check all data in excel, determine which asv are contaminants (e.g. present in negative extraction control and samples from same extraction batch)
contamination <- c("asv.17", "asv.30", "asv.180", "asv.251","asv.453", "asv.473", "asv.1471", "asv.1589", "asv.3038")

# Select all rows in seq_data where the OTU is a contaminant
contaminant_rows <- seq_data %>% 
  filter(`#OTU_ID` %in% contamination)

contaminant_rows_nonzero <- contaminant_rows %>%
  select(where(~ any(. != 0)))   

contaminant_rows_long <- contaminant_rows_nonzero %>%
  pivot_longer(cols = -c(`#OTU_ID`, taxonomy),names_to = "sample", values_to = "reads") %>%
  filter(reads > 0)

cols17 <- c("NIOZ401.167.174",
            "NIOZ401.173.180",
            "NIOZ402.009.018",
            "NIOZ402.011.020")

seq_data[seq_data$"#OTU_ID" == "asv.17", cols17] <- 0

cols30 <- c("NIOZ402.085.094", "NIOZ402.091.100", "NIOZ402.103.112", "NIOZ402.107.116")

seq_data[seq_data$"#OTU_ID" == "asv.30", cols30] <- pmax(seq_data[seq_data$"#OTU_ID" == "asv.30", cols30] - 510, 0)

seq_data <- seq_data %>%
  filter(!(`#OTU_ID` %in% c("asv.251", "asv.180", "asv.453", 
                            "asv.473", "asv.1471", "asv.1589", "asv.3038")))

rm(all_samples, contaminant_rows, contaminant_rows_long, contaminant_rows_nonzero, neg_control, reads_comparison, reads_comparison_2, samples_in_sequencedata, total_reads, cols17, cols30, contamination, contaminants, missing_samples_vec)

#####2. Positive controls#####
#check what groups appear in the positive control
pos_controls <- c("NIOZ398.177.178",
                  "NIOZ399.177.180", 
                  "NIOZ400.177.182", 
                  "NIOZ401.177.184", 
                  "NIOZ402.177.186", 
                  "NIOZ403.135.146")

pos_control <- seq_data %>% select(c('#OTU_ID', "NIOZ398.177.178",
                                     "NIOZ399.177.180", 
                                     "NIOZ400.177.182", 
                                     "NIOZ401.177.184", 
                                     "NIOZ402.177.186", 
                                     "NIOZ403.135.146",, taxonomy)) %>% 
  rowwise() %>% 
  mutate(sum_reads = sum(c_across(-c("#OTU_ID", "taxonomy")))) %>% 
  filter(sum_reads >10) 

#####3. Identity cutoff #### 
seq_data$identity[is.na(seq_data$identity)] <- "0"
seq_data$identity <- as.numeric(seq_data$identity)

seq_data <- seq_data %>% 
  rowwise() %>% 
  mutate(sum_reads = sum(c_across(-c("#OTU_ID", "taxonomy", "identity")))) %>%
  mutate(identity = round(identity)) 

check <- seq_data %>% filter(taxonomy == "Eukaryota;Metazoa;Mollusca;Gastropoda;Littorinimorpha;Littorinidae;Littorina;Littorina littorea")

cutoff <- seq_data %>% 
  select(c("#OTU_ID", taxonomy, sum_reads, identity)) %>% 
  separate(taxonomy, into = c("Domain", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"), sep = ";", fill = "right") %>% 
  mutate(across(c(Domain, Kingdom, Phylum, Class, Order, Family, Genus, Species),~ ifelse(grepl("Unassigned|unknown|unkown", ., ignore.case = TRUE), NA, .)))

data_summary <- cutoff %>%
  pivot_longer(cols = Domain:Species, names_to = "Taxonomic_Rank", values_to = "Taxon") %>%
  filter(!is.na(Taxon)) %>%
  group_by(identity, Taxonomic_Rank, Taxon) %>%
  summarise(Read_Counts = sum(sum_reads, na.rm = TRUE),OTU_Counts = n_distinct(`#OTU_ID`, na.rm = TRUE)) %>%
  mutate(Read_Percentage = Read_Counts / sum(Read_Counts) * 100, OTU_Percentage = OTU_Counts / sum(OTU_Counts) * 10) %>%
  ungroup()

# ensure numeric so 100 > 98 > ... (not character sorting)
data_summary <- data_summary %>% mutate(identity = as.numeric(identity))

unassigned <- data_summary %>%
  select(-Taxon, -Read_Percentage, -OTU_Percentage) %>%
  group_by(identity, Taxonomic_Rank) %>%
  summarise(Read_Counts = sum(Read_Counts),
            ASV_Counts  = sum(OTU_Counts),
            .groups = "drop")

cumulative_count <- unassigned %>%
  group_by(Taxonomic_Rank) %>%
  arrange(desc(identity), .by_group = TRUE) %>%      # HIGH -> LOW for cumsum
  mutate(Cum_Read_Counts = cumsum(Read_Counts),
         Cum_ASV_Counts  = cumsum(ASV_Counts)) %>%
  ungroup()

#Make your own legend 
p1 <- cumulative_count %>% 
  ggplot(aes(x = identity, y = Cum_Read_Counts, color = Taxonomic_Rank, group = Taxonomic_Rank)) +
  geom_line(linewidth = 1, position = "identity") +
  labs(x = "BLASTn Identity Cut-off (%)", y = "Count") +
  xlim(90,100)+
  geom_vline(xintercept=97, colour= "black", linetype="dashed")+
  scale_color_manual(name = 'Taxonomic Groups', breaks = c("Domain", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"),
                     values = c("Domain" = "#b15928", "Kingdom" = "#a6cee3","Phylum" = "#e31a1c", "Class" = "#6a3d9a", "Order" = "#fdae6b", "Family" = "hotpink2", "Genus" = "#1f78b4","Species" = "#33a02c" )) +
  theme_cowplot() + background_grid() +
  theme(plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position='none', axis.title.x = element_blank()) 

p1

p2<-cumulative_count %>% 
  ggplot(aes(x = identity, y = Cum_ASV_Counts, color = Taxonomic_Rank, group = Taxonomic_Rank)) +
  geom_line(linewidth = 1, position = "identity") +
  labs(x = "BLASTn Identity Cut-off (%)", y = "Count") +
  xlim(90,100)+
  geom_vline(xintercept=97, colour= "black", linetype="dashed")+
  scale_color_manual(name = 'Taxonomic Groups', breaks = c("Domain", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"), values = c("Domain" = "#b15928", "Kingdom" = "#a6cee3","Phylum" = "#e31a1c", "Class" = "#6a3d9a", "Order" = "#fdae6b", "Family" = "hotpink2", "Genus" = "#1f78b4","Species" = "#33a02c" ))+
  theme_cowplot() + background_grid() +
  theme(plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position = 'none', legend.title = element_blank(), axis.title.x = element_blank())
p2


#plot_grid(p1, p2, ncol = 2, nrow = 1, align="h")

#proportions of all groups
tot_proportions <- cumulative_count %>% 
  group_by(Taxonomic_Rank) %>% 
  summarise(total_cum_reads = sum(Read_Counts), 
            total_cum_asv = sum(ASV_Counts)) 

proportions <- left_join(cumulative_count, tot_proportions, by="Taxonomic_Rank") %>% 
  group_by(identity, Taxonomic_Rank) %>% 
  mutate(prop_reads = (Cum_Read_Counts/total_cum_reads)*100, 
         prop_asv = (Cum_ASV_Counts/total_cum_asv)*100)

p3 <- proportions %>% 
  ggplot(aes(x = identity, y = prop_reads, color = Taxonomic_Rank, group = Taxonomic_Rank)) +
  geom_line(linewidth = 1, position = "identity") +
  labs(x = "BLASTn Identity Cut-off (%)", y = "Percentage") +
  xlim(90,100)+
  ylim(0,100) +
  geom_vline(xintercept=97, colour= "black", linetype="dashed")+
  scale_color_manual(name = 'Taxonomic Groups', breaks = c("Domain", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"),values = c("Domain" = "#b15928", "Kingdom" = "#a6cee3","Phylum" = "#e31a1c", "Class" = "#6a3d9a", "Order" = "#fdae6b", "Family" = "hotpink2", "Genus" = "#1f78b4","Species" = "#33a02c" )) +
  theme_cowplot() + background_grid() +
  theme(plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position='bottom', axis.title.x = element_blank())

p3

p4 <- proportions %>% 
  ggplot(aes(x = identity, y = prop_asv, color = Taxonomic_Rank, group = Taxonomic_Rank)) +
  geom_line(linewidth = 1, position = "identity") +
  labs(x = "BLASTn Identity Cut-off (%)", y = "Percentage") +
  xlim(90,100)+
  ylim(0,100)+
  geom_vline(xintercept=97, colour= "black", linetype="dashed")+
  scale_color_manual(name = 'Taxonomic Groups', breaks = c("Domain", "Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species"),values = c("Domain" = "#b15928", "Kingdom" = "#a6cee3","Phylum" = "#e31a1c", "Class" = "#6a3d9a", "Order" = "#fdae6b", "Family" = "hotpink2", "Genus" = "#1f78b4","Species" = "#33a02c" )) +
  theme_cowplot() + background_grid() +
  theme( plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position = 'none', legend.title = element_blank(), axis.title.x = element_blank())

#plot_grid(p3, p4,ncol = 2, nrow = 1, align="h")

title_left <- ggdraw() + draw_label("Annotated reads", fontface = 'bold', size = 14)
title_right <- ggdraw() + draw_label("Annotated ASV's", fontface = 'bold', size = 14)

legend_plot <- ggplot(proportions, aes(x = identity, y = prop_asv, color = Taxonomic_Rank)) +
  geom_blank() +   # invisible, but creates legend
  scale_color_manual(
    name = 'Taxonomic Groups',
    breaks = c("Domain","Kingdom","Phylum","Class","Order","Family","Genus","Species"),
    values = c("Domain" = "#8dd3c7",
               "Kingdom" = "#ffffb3", 
               "Phylum" = "#bebada",
               "Class"= "#fb8072", 
               "Order"   = "#7EC0EE", 
               "Family" = "#fdb462",
               "Genus"   = "#b3de69",
               "Species" = "#fccde5")) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.margin = margin(t = -5, b = -5),  
        legend.box.margin = margin(0, 0, 0, 0))

x_axis_title <- cowplot::ggdraw() + draw_label("BLASTn Identity Cut-off (%)", fontface = 'bold', size = 12, hjust = 0.5)

top_titles <- cowplot::plot_grid(title_left, title_right, ncol = 2)

plots <- cowplot::plot_grid(p1, p3, p2, p4, ncol = 4, align = "vh")

middle <- cowplot::plot_grid(top_titles, plots, x_axis_title,
  ncol = 1, rel_heights = c(0.15, 1, 0.08))

final <- cowplot::plot_grid(middle,legend_plot, ncol = 1, rel_heights = c(1, 0.08))

final

seq_data <- seq_data %>% 
  filter(identity >= 97.0)

#####3. Aggregate data #####
names(seq_data)
head(row.names(seq_data))

unique(seq_data$taxonomy)

ASV<-seq_data[ ,2:633]

TaxTab<- seq_data[635] #select the taxonomy for each ASV-number 
names(TaxTab)
head(TaxTab)

TaxTab <- TaxTab %>%
  separate(taxonomy, into = c("Domain", "Kingdom", "Phylum", "Class","Order", "Family", "Genus", "Species"), sep = ";") 

TaxTab$Domain[TaxTab$Domain == "Unassigned"] <- NA
TaxTab$Species <- gsub(".*(unknown_).*", "\\1", TaxTab$Species, ignore.case = TRUE)
TaxTab$Species[TaxTab$Species == "unknown_"] <- NA
TaxTab$Genus <- gsub(".*(unknown_).*", "\\1", TaxTab$Genus, ignore.case = TRUE)
TaxTab$Genus[TaxTab$Genus == "unknown_"] <- NA
TaxTab$Family <- gsub(".*(unknown_).*", "\\1", TaxTab$Family, ignore.case = TRUE)
TaxTab$Family[TaxTab$Family == "unknown_"] <- NA

TaxTab$L8<-str_c(TaxTab$Domain,';',TaxTab$Kingdom,';',TaxTab$Phylum,';',TaxTab$Class,';',TaxTab$Order,';',TaxTab$Family,';',TaxTab$Genus, ';',TaxTab$Species) #This code creates a new column called L1_2_3_4_5_6_7 by combining (or concatenating) several existing columns. 

TaxTab$L7<-str_c(TaxTab$Domain,';',TaxTab$Kingdom,';',TaxTab$Phylum,';',TaxTab$Class,';',TaxTab$Order,';',TaxTab$Family,';',TaxTab$Genus)

TaxTab$L6<-str_c(TaxTab$Domain,';',TaxTab$Kingdom,';',TaxTab$Phylum,';',TaxTab$Class,';',TaxTab$Order,';',TaxTab$Family) #Same as above, but then without L7

TaxTab$L5<-str_c(TaxTab$Domain,';',TaxTab$Kingdom,';',TaxTab$Phylum,';',TaxTab$Class,';',TaxTab$Order)

TaxTab$L4<-str_c(TaxTab$Domain,';',TaxTab$Kingdom,';',TaxTab$Phylum,';',TaxTab$Class)

TaxTab$L3<-str_c(TaxTab$Domain,';',TaxTab$Kingdom,';',TaxTab$Phylum)

TaxTab$L2<-str_c(TaxTab$Domain,';',TaxTab$Kingdom)

dim(TaxTab)

#Merge with same/asv data, pick how specific you want the taxonomy
#For me, tax6 seems the best at the moment? 
TaxL2aggr<-aggregate(ASV, by =list(TaxTab$L2), FUN=sum)
TaxL3aggr<-aggregate(ASV, by =list(TaxTab$L3), FUN=sum)
TaxL4aggr<-aggregate(ASV, by =list(TaxTab$L4), FUN=sum)
TaxL5aggr<-aggregate(ASV, by =list(TaxTab$L5), FUN=sum)
TaxL6aggr<-aggregate(ASV, by =list(TaxTab$L6), FUN=sum)
TaxL7aggr<-aggregate(ASV, by =list(TaxTab$L7), FUN=sum)
TaxL8aggr<-aggregate(ASV, by =list(TaxTab$L8), FUN=sum)

#I wont use anything less specific than family-level, so use TAXL6 - L8. Lets first start with L6 

curculionidae_sp <- TaxL8aggr %>% 
  separate(Group.1, into = c("Domain", "Kingdom", "Phylum", "Class","Order", "Family", "Genus", "Species"), sep = ";") %>% 
  filter(Family == "Curculionidae") 

Family_tax <- TaxL6aggr %>% separate(Group.1, into = c("Domain", "Kingdom", "Phylum", "Class","Order", "Family"), sep = ";") 

Family_tax <- Family_tax %>%
  filter(Domain == "Eukaryota") %>% 
  filter(Kingdom == "Metazoa") %>% 
  filter(!Phylum == "Chordata")

unique(Family_tax$Family)

#library(writexl)
#write_xlsx(Family_tax, "./taxonomy_family_level.xlsx")

#####2. Loading data & first impression#####
rm(ASV, cumulative_count, curculionidae_sp, cutoff, data_summary, legend_plot, p1, p2, p3, p4, proportions, seq_data, TaxL2aggr, TaxL3aggr, TaxL4aggr, TaxL5aggr, TaxL6aggr, TaxL7aggr, TaxL8aggr, TaxTab, title_left, title_right, tot_proportions, unassigned, x_axis_title, final, middle, plots, top_titles, pos_control)


#load ASV file, remove controls (negative and positive), calculate the sum reads for each family, and filter out those families that have very little reads (<10)

data <- Family_tax %>% 
  select(!all_of(neg_controls)) %>%
  select(!all_of(pos_controls)) %>%
  select(-starts_with("NIOZ358.")) %>%
  rowwise() %>% 
  mutate(sum_reads = sum(c_across(-c("Domain", "Kingdom", "Phylum", "Class", "Order", "Family")))) %>% 
  filter(!sum_reads < 10)

unique(data$Family)

#Plot the total number of reads to see which families occur the most
data %>% 
  arrange(desc(sum_reads)) %>%
  ggplot(aes(x=Family, y=sum_reads, fill=Family)) + 
  geom_bar(stat="identity") +
  theme_cowplot() + background_grid() +
  theme(plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position='none', axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1))

#Transform the data into a long-format 
fam_data_long <- data %>% 
  select(!sum_reads) %>%
  mutate(taxonomy = str_c(Order, '-', Family), .before=NIOZ398.003.004) %>%
  select(!c(Domain, Kingdom, Phylum, Class, Order, Family)) %>%
  pivot_longer(!taxonomy, names_to = "Sample", values_to = "count") 

#calculate the total reads per sample as well as the proportion of the sample each family represents 
fam_data_long <- fam_data_long %>% 
  group_by(Sample) %>%
  mutate(sum_reads_sample = sum(count)) %>% 
  filter(sum_reads_sample >= 10) %>%
  ungroup() %>% 
  mutate(RRA = (count/sum_reads_sample)*100)

#Plot of the diet (per sample) with all families. 
fam_data_long %>% 
  ggplot(aes(x=Sample, y=RRA, fill=taxonomy)) + 
  geom_bar(stat="identity") +
  labs(x = "Sample",
       y = "% Diet") +
  theme_cowplot() + background_grid() +
  theme(plot.title = element_text(hjust = 0.5), text = element_text(size = 18), legend.position='bottom') 

data_97 <- unique(fam_data_long$taxonomy) 

#Families that represent => 75% of the diet in a sample
#abundant_groups <- fam_data_long %>% filter(RRA >= 75)

#unique(abundant_groups$taxonomy) #29 different families 

#Change NA's in RRA to 0
fam_data_long$RRA[fam_data_long$RRA == "NaN"] <- 0

#####3. top 25 families #####
rel_reads <- fam_data_long %>% 
  select(c(taxonomy, RRA)) %>%
  group_by(taxonomy) %>%
  summarise(Relative_read_abundance = mean(RRA)) 

families_25 <- rel_reads %>% 
  filter(Relative_read_abundance >= 0.5) #families that represent over 0.5% of the total diet 

#Make a list of the top 25 families, and one with the families that make up over 75% of a sample. Then merge these lists together to create a list with all families, using unique to remove families that occur in both lists. 
list1 <- families_25$taxonomy
#list2 <- unique(abundant_groups$taxonomy) 
#list3 = c(list1, list2) 
families <- unique(list1)
rm(list1)

#select only the families with high abundance
fam_data_long_2 <- fam_data_long %>% 
  filter(taxonomy %in% families) %>% 
  filter(!taxonomy == "Diptera-Calliphoridae") %>% 
  select(!sum_reads_sample) 

#Select all the families that do not have a high abundance, and merge them under "Other"
others <- fam_data_long %>% 
  filter(!taxonomy %in% families | taxonomy == "Diptera-Calliphoridae") %>%  
  group_by(Sample) %>% 
  summarise(count = sum(count), .groups = "drop") %>% 
  mutate(taxonomy = "Other") %>% 
  left_join(
    fam_data_long %>% 
      group_by(Sample) %>% 
      summarise(sum_reads_sample = sum(count), .groups = "drop"),
    by = "Sample"
  ) %>% 
  mutate(RRA = (count / sum_reads_sample) * 100) %>% 
  select(-sum_reads_sample)

unique(others$taxonomy)  
#merge the datasets together 
fam_long <- rbind(fam_data_long_2, others)

unique(fam_long$taxonomy) 

#library(writexl)
#write_xlsx(fam_long, "./filtered_family_data_long_new.xlsx")