--- 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") ```