---
title: "CAZyme analysis"
author: "Rebecca Gellman"
date: "2023-09-12"
output: html_document
---
This document contains the code used to produce figures related to CAZyme annotation in Figures 2 and 3.

Load relevant packages
```{r}
library("tidyverse")
library("readxl")
library("cowplot")
library("stringr")
library("RColorBrewer")
library("MetBrewer")
library("svglite")
library("dendsort")
library("seriation")
library("dendextend")
library("ComplexHeatmap")
library("readxl")
library("circlize")
```


Read in files and format 
Load in .csv files, name the columns, add a "Species" column so that we know what isolate each gene came from. Repeat with all genomes
```{r}
pc2379 <- read_csv("~/2379_Pcopri.csv", c("Gene_ID", "CAzyme")) %>% 
  mutate(isolate = "P. copri 2379") %>% 
  mutate(origin = "Hadza") %>% 
  separate(CAzyme, into = c("CAzyme", "CAzyme02", "CAzyme03"), sep = "-") %>% 
  separate(CAzyme, into = c("CAzyme_01", "subfamily_01"), sep = "_") %>% 
  separate(CAzyme02, into = c("CAzyme_02", "subfamily_02"), sep = "_") %>% 
  select(-X3)

#repeat with all genomes. adding another as example
pc2383 <- read_csv("~/2383_Pcopri.csv", c("Gene_ID", "CAzyme")) %>% 
  mutate(isolate = "P. copri 2383") %>% 
  mutate(origin = "Hadza") %>% 
  separate(CAzyme, into = c("CAzyme", "CAzyme02", "CAzyme03"), sep = "-") %>% 
  separate(CAzyme, into = c("CAzyme_01", "subfamily_01"), sep = "_") %>% 
  separate(CAzyme02, into = c("CAzyme_02", "subfamily_02"), sep = "_")%>%
  select(-X3)

#combine into a single df
pc_combo_df <- rbind(pc2379,pc2383,pc2446,pc2477,pc2489,pc2497,pc2632,pc_dsm18205_western, pc_yf2_western) %>%  #combine hadza bugs
  mutate(type = str_sub(CAzyme_01, 1,2)) %>% 
  filter(type != "CB") %>% 
  filter(type != "GT") %>% 
  filter(type != "CE") #remove CAZyme types that aren't GH or PL

pc_pruned <-  pc_combo_df %>% 
  select(CAzyme_01, isolate, origin, type) 

bp_combo_all <- pc_pruned %>% 
  mutate(isolate2 = isolate) %>% 
  separate(isolate2, into = c("genus", "species", "strain"), sep = " ") %>% 
  left_join(new_ids, by = c("isolate"= "isolate"))

#write.csv(bp_combo_all, "~/all_CAzymes_01.csv")
```
This will generate a df with all CAZymes from all isolates of interest

Load data for making visualizations
```{r}
bp_combo_all <- bp_combo_all %>% 
  select(CAzyme_01, isolate, origin, type, genus, species, strain) 

carb_types <- read_excel("/Users/rebeccagellman/Desktop/Stanford M&I/Sonnenburg Lab/results/CAzyme_list_update.xlsx") %>% 
  mutate(CAzyme_01 = CAzyme)

genome_sizes <- read_xlsx("/Users/rebeccagellman/Desktop/Stanford M&I/Sonnenburg Lab/results/genome_summary.xlsx") 

new_ids <- read_xlsx("/Users/rebeccagellman/Desktop/Stanford M&I/Sonnenburg Lab/results/strain_ids.xlsx")
```

Fig. 3 CAZyme heatmap. 
Load in relevant dfs
```{r}
bp_combo_all <- bp_combo_all %>% 
  select(CAzyme_01, isolate, origin, type, genus, species, strain)

bug_data <- bp_combo_all %>% 
  left_join(new_ids, by = "isolate") %>% 
  select(-CAzyme_01, -type) 
bug_data <- bug_data[!duplicated(bug_data$isolate_id), ]

carb_types <- read_excel("/Users/rebeccagellman/Desktop/Stanford M&I/Sonnenburg Lab/results/CAzyme_list_update.xlsx") %>% 
  mutate(CAzyme_01 = CAzyme) %>%
  select(CAzyme_01, Broad_Substrate)
carb_types$Broad_Substrate[is.na(carb_types$Broad_Substrate)] <- "Other"
carb_types <- carb_types[!duplicated(carb_types$CAzyme_01), ]
```

Make a matrix to go into complex heatmap
```{r}
rep.zero <- function(x) ifelse(is.na(x), 0, x) #function to turn NA into 0

tmp <- bp_combo_all %>% 
  select(-isolate_id) %>% 
  left_join(new_ids, by = "isolate")

tmp1 <- select(tmp,CAzyme_01) %>% 
  group_by(CAzyme_01) %>% 
  tally(name = "common") %>% 
  ungroup() 


t1 <- tmp %>% 
  left_join(tmp1, by = c('CAzyme_01' = 'CAzyme_01')) %>% 
  #arrange(-common) %>% 
  group_by(isolate_id, CAzyme_01) %>% 
  summarise(n = n()) %>% 
  select(CAzyme_01, isolate_id, n) %>%
  pivot_wider(names_from = isolate_id,
              values_from = n) %>% 
  mutate_at(vars(!starts_with("C")), rep.zero)


t3 <- t1 %>% 
  mutate_at(vars(!starts_with("C")), rep.zero) 

CAzyme <- t3$CAzyme_01

t4 <- select(t3, -CAzyme_01) #%>%
t4 <- as.matrix(t4)
rownames(t4) <- CAzyme
# select(-CAzyme_01)

substrate <- as_tibble(t3) %>%   
  left_join(carb_types, by = "CAzyme_01") %>% 
  select(-CAzyme)
substrate$Broad_Substrate[is.na(substrate$Broad_Substrate)] <- "Other"
substrate <- substrate[!duplicated(substrate$CAzyme_01), ]

mat <- t4   #this is the matrix for making the heatmap

```

Make heatmap using ComplexHeatmap
```{r}
#manual color palettes
col_fun = colorRamp2(c(0, 10, 36), c("#ffffff", "#707fdb", "#151b43"))
bread <- list(Substrate = c("Animal" = "#ef9d5f", "Other" = "#707fdb", "Plant" = "#151b43"))
genera <- colorRamp2(c("B.", "P."), c("#c9421e", "#efd05f"))
  list(genus = c("B." = "#c9421e", "P." = "#efd05f", origin = c("Hadza" = "#707fdb", "reference" = "#151b43")))
origins <- colorRamp2(c("Hadza", "reference"), c("black", "gray"))
  list(origin = c("Hadza" = "orange", "reference" = "blue"))


#define annotations
b_annotation = HeatmapAnnotation(origin = bug_data$origin,
                                 genus = bug_data$genus,
                                 col = list(origin = c("Hadza" = "black", "reference"= "grey"), genus = c("B."="#c9421e", "P."= "#efd05f")),
                                 show_annotation_name = FALSE, show_legend = FALSE,
                                 annotation_legend_param = list(genus = list(title = "Genus", labels = c("Bacteroides", "Prevotella")), 
                                                                                              origin = list(title = "Origin", labels = c("Hadza", "Reference"))))



lgd = Legend(col_fun = col_fun, title = "foo", at = c(0, 20, 40), 
             legend_height = unit(4, "cm"))

hta = ComplexHeatmap::Heatmap(mat, 
             col = col_fun,
             show_row_dend = TRUE,
             show_row_names = TRUE,
             row_names_gp = gpar(fontsize= 9),
             column_names_gp = gpar(fontsize= 9),
             show_column_dend = TRUE,
             show_heatmap_legend = FALSE,
             use_raster = TRUE,
             bottom_annotation = b_annotation,
             raster_by_magick = TRUE)

ComplexHeatmap::draw(hta)
dev.set(dev.next())
dev.off()
```


Fig. 3 CAZymes by size
```{r}
cazymes_by_size <- bp_combo_all %>% 
  group_by(isolate) %>% 
  summarise(cazy_count = n()) %>% 
  #mutate(isolate  = isolate) %>% 
  inner_join(genome_sizes, by = "isolate") %>% 
  drop_na(isolate) %>% 
  mutate(cazy_by_mb = cazy_count/genome_size_mb) %>% 
  mutate(cazy_by_bp = cazy_count/genome_size_bp) %>% 
  mutate(cazy_by_gene  = cazy_count/num_genes)

fig3b <- ggplot(by_type, aes(x = isolate_id, y = cazy_by_mb)) +geom_col(aes(fill= factor(Broad_Substrate, levels = c("Plant", "Other", "Animal")))) +theme_classic()+ theme(text = element_text(family = "Arial", size = 16), axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), axis.title.x = element_blank(), axis.text.x = element_text(face= "italic", angle = 60, vjust = .95, hjust = 1)) + scale_fill_manual(values = general_carbs, name = "Substrate Category") + ylab("CAZymes by Mb")+ theme(legend.position = "none")
print(fig3b)
###THIS IS 3C
fig3c <- ggplot(by_type, aes(x = isolate_id, y = percent)) +geom_col(aes(fill= factor(Broad_Substrate, levels = c("Plant", "Other", "Animal")))) +theme_classic()+ theme(text = element_text(family = "Arial", size = 16), axis.ticks.x = element_blank(), axis.ticks.y = element_blank(), axis.title.x = element_blank(), axis.text.x = element_text(face= "italic", angle = 60, vjust = .95, hjust = 1)) + scale_fill_manual(values = general_carbs) + ylab("Percent CAZymes") + labs(fill = "Substrate Category") + theme(legend.position = "none")
```

Fig. 3 Mucin CAZymes
```{r}
mucin <- cazymes_by_type %>%  
  filter(mucin == "Mucin") %>%   #filter only mucin cazymes
  #filter(origin == "Hadza") %>% 
  group_by(isolate, CAzyme) %>% 
  # summarise(cazy_counts = n()) %>% 
  # ungroup() %>% 
  left_join(new_ids, by = "isolate")

mucin_no_gh2 <- mucin %>% 
  filter(CAzyme == "GH18" | CAzyme == "GH20") %>% #only GH18 and GH20
  inner_join(genome_sizes, by = "isolate") %>% 
  left_join(new_ids, by = "isolate") %>% 
  mutate(cazy_by_mb = cazy_counts/genome_size_mb)

mucin_all <- mucin %>% 
  filter(CAzyme == "GH18" | CAzyme == "GH20") %>% #only GH18 and GH20
  inner_join(genome_sizes, by = "isolate") %>% 
  left_join(new_ids, by = "isolate") %>% 
  mutate(cazy_by_mb = cazy_counts/genome_size_mb)

fig3d <- ggplot(mucin_no_gh2, aes(x = factor(isolate_id), y = cazy_counts, fill = CAzyme)) + geom_bar(stat = "identity") + theme_classic() + theme(text = element_text(family = "Arial", size = 16)) + theme(axis.title.x = element_blank()) + theme(axis.text.x = element_text(angle = 90, vjust = .95, hjust = 1, face = "italic")) + scale_fill_manual(values=animal_carbs) + ylab("Number of Genes") + theme(legend.position = "none")

figs3 <- ggplot(mucin_all, aes(x = factor(isolate_id), y = cazy_counts, fill = CAzyme)) + geom_bar(stat = "identity") + theme_classic() + theme(text = element_text(family = "Arial", size = 16)) + theme(axis.title.x = element_blank()) + theme(axis.text.x = element_text(angle = 90, vjust = .95, hjust = 1, face = "italic")) + scale_fill_manual(values=animal_carbs) + ylab("Number of Genes") + theme(legend.position = "none")
```

