# Leaf-shelters facilitate the colonization of arthropods and enhance microbial diversity on plants

# This heatmap script will use the relative abundance of bacterial taxa to show visualy the relative abundance of different taxa between control and leaf shelter

# Cleaning the script
rm(list = ls(all = TRUE))
graphics.off()

# Library
library(tidyverse)
library(readxl)
library(ggtext)
library(RColorBrewer)

# Trial 2 - Relative abundance of bacterial taxa ####
# Directory/Data
getwd()
setwd("/Users/daniloferreiraborgesdossantos/Library/CloudStorage/OneDrive-ThePennsylvaniaStateUniversity/Documents/Sync/Research/PhD/PennState/Papers/LeafShelter/Danilo_edits/Analyzes/LefSe")
dir()

# This file will provide the plant ID and the specific OTUs associated to each sample
metadata <- read_excel("leaf.shelter.metadata.T2.xlsx", na="NA") %>%
 dplyr::select(plant_id, treatment) %>%
 drop_na(treatment)

# This is the abundance table that will be used to calculate the relative abundance of OTUs
otu_counts <- read_table("otu_table_T2.txt")%>%
 dplyr::select(otu, starts_with("lf"))%>%
 pivot_longer(-otu, names_to="plant_id", values_to = "count")

# Taxonomy will assign the species name for the specific OTU present in each sample
taxonomy <- read_tsv("taxonomy.T2.tsv")%>%
 dplyr::select("otu", Taxon)%>%
 rename_all(tolower)%>%
 mutate(taxon = str_replace_all(taxon, "\\(\\d+\\)", ""),
        taxon = str_replace(taxon, ";$", ""))%>%
 separate(taxon,
          into=c("kingdom", "phylum", "class", "order", "family", "genus", 
                 "species"),
          sep=";")

# Calculating the OTU relative abundance
otu_rel_abund <- inner_join(metadata, otu_counts, by="plant_id") %>%
 inner_join(., taxonomy, by="otu") %>%
 group_by(plant_id)%>%
 mutate(rel_abund = 100*count / sum(count))%>%
 ungroup() %>%
 dplyr::select(-count) %>%
 pivot_longer(
  c("kingdom", "phylum", "class", "order", "family", "genus", "species"),
  names_to="level",
  values_to="taxon")

## Family ####
# This part will calculate the relative abundance for bacterial families
taxon_rel_abund <- otu_rel_abund %>%
 filter(level=="family") %>%
 group_by(treatment, plant_id, taxon) %>%
 dplyr::summarize(rel_abund = sum(rel_abund), .groups="drop") %>%
 group_by(treatment, taxon) %>%
 dplyr::summarise(mean_rel_abund = mean(rel_abund), .groups = "drop")%>%
 dplyr::mutate(taxon = str_replace(taxon,
                                   "(.*)_unclassified", "Unclassified<br>*\\1*"),
               taxon = str_replace(taxon,
                                   "^([^<]*)$", "*\\1*"),
               taxon = str_replace_all(taxon,
                                       "_", " "))
# This command summarizes max and mean of the relative abundance for families
taxon_pool <- taxon_rel_abund %>%
 group_by(taxon) %>%
 summarize(pool = max(mean_rel_abund) < .5,
           mean = mean(mean_rel_abund),
           .groups="drop")

# Plot
inner_join(taxon_rel_abund, taxon_pool, by="taxon") %>%
 mutate(taxon = if_else(pool, "Other", taxon)) %>%
 group_by(treatment, taxon) %>%
 summarize(mean_rel_abund = sum(mean_rel_abund),
           mean = min(mean),
           .groups="drop") %>%
 mutate(taxon = factor(taxon),
        taxon = fct_reorder(taxon, mean, .desc=FALSE)) %>%
 ggplot(aes(x=treatment, y=taxon, fill=mean_rel_abund)) +
 geom_tile() +
 geom_text(aes(label = format(round(mean_rel_abund, 1), nsmall = 1))) +
 theme_classic() +
 scale_fill_gradient(name="Mean<br>Relative<br>Abun (%)",
                     low = "#FFFFFF",
                     high = "#FF0000",
                     expand = c(0,0),
                     limits=c(0,NA)) +
 theme(axis.text.x = element_markdown(),
       legend.key.size = unit(20, "pt"),
       legend.title = element_markdown(size = 15),
       legend.title.align = .5,
       axis.line = element_blank(),
       axis.ticks = element_blank(),
       legend.text = element_text(colour="black", size=15),
       axis.text = element_text(size = 15, color = "black")) +
 labs(x=NULL,
      y=NULL) +
 scale_x_discrete(expand=c(0, 0)) +
 coord_fixed(ratio = .3) +
 scale_y_discrete(labels = c("Other", "Acidobacteriaceae", "Neisseriaceae", "Caulobacteracea", "Pseudomonadaceae", "Armatimonadaceae", "Comamonadacea", "Chitinophagaceae", "Acetobacteraceae", "Sphingobacteriaceae", "Sphingomonadaceae", "Methylobacteriaceae"))

ggsave("MICROBIOME.HEAT.MAP.T2.family.tiff", width=8, height=10, dpi = 300)

## Genus ####
# This part will calculate the relative abundance for bacterial genus
taxon_rel_abund <- otu_rel_abund %>%
 filter(level=="genus") %>%
 group_by(treatment, plant_id, taxon) %>%
 summarize(rel_abund = sum(rel_abund), .groups="drop") %>%
 group_by(treatment, taxon) %>%
 summarise(mean_rel_abund = mean(rel_abund), .groups = "drop")%>%
 mutate(taxon = str_replace(taxon,
                            "(.*)_unclassified", "Unclassified<br>*\\1*"),
        taxon = str_replace(taxon,
                            "^([^<]*)$", "*\\1*"),
        taxon = str_replace_all(taxon,
                                "_", " "))
# This command summarizes max and mean of the relative abundance for genus
taxon_pool <- taxon_rel_abund %>%
 group_by(taxon) %>%
 summarize(pool = max(mean_rel_abund) < .5,
           mean = mean(mean_rel_abund),
           .groups="drop")

# Plot
inner_join(taxon_rel_abund, taxon_pool, by="taxon") %>%
 mutate(taxon = if_else(pool, "Other", taxon)) %>%
 group_by(treatment, taxon) %>%
 summarize(mean_rel_abund = sum(mean_rel_abund),
           mean = min(mean),
           .groups="drop") %>%
 mutate(taxon = factor(taxon),
        taxon = fct_reorder(taxon, mean, .desc=FALSE)) %>%
 ggplot(aes(x=treatment, y=taxon, fill=mean_rel_abund)) +
 geom_tile() +
 geom_text(aes(label = format(round(mean_rel_abund, 1), nsmall = 1))) +
 theme_classic() +
 scale_fill_gradient(name="Mean<br>Relative<br>Abun (%)",
                     low = "#FFFFFF",
                     high = "#FF0000",
                     expand = c(0,0),
                     limits=c(0,NA)) +
 theme(axis.text.x = element_markdown(),
       legend.key.size = unit(20, "pt"),
       legend.title = element_markdown(size = 15),
       legend.title.align = .5,
       axis.line = element_blank(),
       axis.ticks = element_blank(),
       legend.text = element_text(colour="black", size=15),
       axis.text = element_text(size = 15, color = "black"),
       axis.text.y = element_text(size = 15, face="italic")) +
 labs(x=NULL,
      y=NULL) +
 scale_x_discrete(expand=c(0, 0)) +
 coord_fixed(ratio = .3) +
 scale_y_discrete(labels = c("Other", "Granulicella", "Silvimonas", "Pseudomonas", "Mucilaginibacter", "Variovax", "Sphingomonas", "Methylobacterium", "Methylobacteriaceae"))

ggsave("MICROBIOME.T2.GENUS.PDF", width=8, height=10, dpi = 500)

# Trial 2 - Relative abundance of bacterial taxa 
# Cleaning the script
rm(list = ls(all = TRUE))
graphics.off()

# Library
library(tidyverse)
library(readxl)
library(ggtext)
library(RColorBrewer)

# Directory/Data
getwd()
setwd("/Users/daniloferreiraborgesdossantos/Library/CloudStorage/OneDrive-ThePennsylvaniaStateUniversity/Documents/Sync/Research/PhD/PennState/Papers/LeafShelter/Danilo_edits/Analyzes/Updated_data")
dir()

# This file will provide the plant ID and the specific OTUs associated to each sample
metadata <- read.csv("metadata.T1.csv", na="NA") %>%
 dplyr::select(plant_id, treatment) %>%
 drop_na(treatment)

# This is the abundance table that will be used to calculate the relative abundance of OTUs
otu_counts <- read_table("otu.count.T1.txt") %>%
 dplyr::select(otu, starts_with("lf")) %>%
 pivot_longer(-otu, names_to="plant_id", values_to = "count")

# Taxonomy will assign the species name for the specific OTU present in each sample
taxonomy <- read_tsv("taxonomy.lefse.T1.tsv")%>%
 dplyr::select(otu, Taxon) %>%
 rename_all(tolower) %>%
 mutate(taxon = str_replace_all(taxon, "\\(\\d+\\)", ""),
        taxon = str_replace(taxon, ";$", "")) %>%
 separate(taxon,
          into=c("kingdom", "phylum", "class", "order", "family", "genus", 
                 "species"),
          sep=" ")

# Calculating the OTU relative abundance
otu_rel_abund <- inner_join(metadata, otu_counts, by="plant_id")%>%
 inner_join(., taxonomy, by="otu")%>%
 dplyr::group_by(plant_id)%>%
 dplyr::mutate(rel_abund = 100*count / sum(count))%>%
 ungroup()%>%
 dplyr::select(-count)%>%
 pivot_longer(
  c("kingdom", "phylum", "class", "order", "family", "genus", "species"),
  names_to="level",
  values_to="taxon")

## Family ####
# This part will calculate the relative abundance for bacterial families
taxon_rel_abund <- otu_rel_abund %>%
 filter(level=="family") %>%
 group_by(treatment, plant_id, taxon) %>%
 dplyr::summarize(rel_abund = sum(rel_abund), .groups="drop") %>%
 group_by(treatment, taxon) %>%
 dplyr::summarise(mean_rel_abund = mean(rel_abund), .groups = "drop")%>%
 dplyr::mutate(taxon = str_replace(taxon,
                                   "(.*)_unclassified", "Unclassified<br>*\\1*"),
               taxon = str_replace(taxon,
                                   "^([^<]*)$", "*\\1*"),
               taxon = str_replace_all(taxon,
                                       "_", " "))
# This command summarizes max and mean of the relative abundance for families
taxon_pool <- taxon_rel_abund %>%
 group_by(taxon) %>%
 summarize(pool = max(mean_rel_abund) < .5,
           mean = mean(mean_rel_abund),
           .groups="drop")

# Plot
inner_join(taxon_rel_abund, taxon_pool, by="taxon") %>%
 mutate(taxon = if_else(pool, "Other", taxon)) %>%
 group_by(treatment, taxon) %>%
 summarize(mean_rel_abund = sum(mean_rel_abund),
           mean = min(mean),
           .groups="drop") %>%
 mutate(taxon = factor(taxon),
        taxon = fct_reorder(taxon, mean, .desc=FALSE)) %>%
 ggplot(aes(x=treatment, y=taxon, fill=mean_rel_abund)) +
 geom_tile() +
 geom_text(aes(label = format(round(mean_rel_abund, 1), nsmall = 1))) +
 theme_classic() +
 scale_fill_gradient(name="Mean<br>Relative<br>Abun (%)",
                     low = "#FFFFFF",
                     high = "#FF0000",
                     expand = c(0,0),
                     limits=c(0,NA)) +
 theme(axis.text.x = element_markdown(),
       legend.key.size = unit(20, "pt"),
       legend.title = element_markdown(size = 15),
       legend.title.align = .5,
       axis.line = element_blank(),
       axis.ticks = element_blank(),
       legend.text = element_text(colour="black", size=15),
       axis.text = element_text(size = 15, color = "black")) +
 labs(x=NULL,
      y=NULL) +
 scale_x_discrete(expand=c(0, 0)) +
 coord_fixed(ratio = .3) +
 scale_y_discrete(labels = c("Other", "Neisseriaceae", "Rhizobiaceae", "Comamonadaceae", "Acidobacteriaceae", "Sphingobacteriaceae", "Pseudomonadaceae", "Microbacteriaceae", "Oxalobacteraceae", "Methylocystaceae"))

ggsave("MICROBIOME.HEATMAP.T1.Family.tiff", width=8, height=10, dpi = 350)

## Genus ####
# This part will calculate the relative abundance for bacterial genus
taxon_rel_abund <- otu_rel_abund %>%
 filter(level=="genus") %>%
 group_by(treatment, plant_id, taxon) %>%
 dplyr::summarize(rel_abund = sum(rel_abund), .groups="drop") %>%
 group_by(treatment, taxon) %>%
 dplyr::summarise(mean_rel_abund = mean(rel_abund), .groups = "drop")%>%
 dplyr::mutate(taxon = str_replace(taxon,
                                   "(.*)_unclassified", "Unclassified<br>*\\1*"),
               taxon = str_replace(taxon,
                                   "^([^<]*)$", "*\\1*"),
               taxon = str_replace_all(taxon,
                                       "_", " "))

# This command summarizes max and mean of the relative abundance for genus
taxon_pool <- taxon_rel_abund %>%
 group_by(taxon) %>%
 summarize(pool = max(mean_rel_abund) < .5,
           mean = mean(mean_rel_abund),
           .groups="drop")

# Plot
inner_join(taxon_rel_abund, taxon_pool, by="taxon") %>%
 mutate(taxon = if_else(pool, "Other", taxon)) %>%
 group_by(treatment, taxon) %>%
 summarize(mean_rel_abund = sum(mean_rel_abund),
           mean = min(mean),
           .groups="drop") %>%
 mutate(taxon = factor(taxon),
        taxon = fct_reorder(taxon, mean, .desc=FALSE)) %>%
 ggplot(aes(x=treatment, y=taxon, fill=mean_rel_abund)) +
 geom_tile() +
 geom_text(aes(label = format(round(mean_rel_abund, 1), nsmall = 1))) +
 theme_classic() +
 scale_fill_gradient(name="Mean<br>Relative<br>Abun (%)",
                     low = "#FFFFFF",
                     high = "#FF0000",
                     expand = c(0,0),
                     limits=c(0,NA)) +
 theme(axis.text.x = element_markdown(),
       legend.key.size = unit(20, "pt"),
       legend.title = element_markdown(size = 15),
       legend.title.align = .5,
       axis.text.y = element_text(face = "italic"),
       axis.line = element_blank(),
       axis.ticks = element_blank(),
       legend.text = element_text(colour="black", size=15),
       axis.text = element_text(size = 15, color = "black")) +
 labs(x=NULL,
      y=NULL) +
 scale_x_discrete(expand=c(0, 0)) +
 coord_fixed(ratio = .3) +
 scale_y_discrete(labels = c("Other", "Granulicella", "Silvimonas", "Agrobacterium", "Variovax", "Acidobacteriaceae", "Sphingobacteriaceae", "Pseudomonas", "Microbacteriaceae", "Massilia", "Methylocystaceae"
 ))

ggsave("MICROBIOME.HEATMAP.T1.genus.PDF", width=8, height=10, dpi = 500)