--- title: "HeadspaceAnalysisV2TICnorm" author: "SML" date: "3/26/2020" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` This analysis is for the GCMS headspace data collected at Supra. The peak areas are from confirmed compounds with standards ran as well as some putitativly identified peaks. The peak areas were normalized to the TIC for each sample, therefore the values being analyzed in this dataset are RELATIVE ABUNDANCES and NOT absolute quantitative values. I removed the peak area for 1-octanol in this analysis, it was added as an internal standard but there was not enough sensitvity without the SPME fiber to use this as an internal standard. Goals of the analysis: - Create some exploratory graphs to look for outliers -> PCA plots, colored by different factors - full dataset - 15C samples - 25C samples - Comparison of strains WITHIN each temperature: - subset into 15C and 25C - do an ANOVA (for each compound, comparing between 11 strains + blank) - Tukey HSD post hoc analysis for each compound within each temperature. - correlation analysis for each temperature group. ```{r message=FALSE} library(readxl) library(tigerstats) library(tidyverse) library(ggplot2) library(ggpubr) library(visreg) library(car) library(agricolae) library(devtools) #library(ggbiplot) library(xlsx) library(corrplot) library(Hmisc) library(cowplot) library(ggcorrplot) library(RColorBrewer) library(gplots) library(ggfortify) library(ggthemes) ``` Read in the data: ```{r} df<-read_xlsx("data/revised peak areas 20200314.xlsx", sheet= "TIC_Norm_areas") %>% data.frame() df <- filter(df, strain != "curve") df <- filter(df, strain != "Model_Wine") df <- filter(df, strain != "Blank") df <- filter(df, strain != "QC") df$aucl <- as.numeric(df$aucl) full <- df hot <- filter(df, temperature == "25") cold <- filter(df, temperature =="15") ``` _______________________________________________________________________________________________________________________________ PCA Plots for whole dataset. - no blanks or QC's - no loadings - saved as jpegs ######PCA for Full df colored by strain and temp#################### ```{r} #?theme() #change default of ggplot2 to center the title. #theme_update(plot.title = element_text(hjust = 0.5)) #do the PCA full_pcaformat <- full[6:ncol(full)] %>% data.frame() row.names(full_pcaformat) <- full$sample_ID full_pca <- prcomp(full_pcaformat, center = TRUE, scale. = TRUE) #summary(cold_pca) #Colored by strain #strain <- full$strain all_species_pca <- autoplot(full_pca, data = full, colour = 'species', scale = T, size = 1.25, loadings = FALSE, loadings.label = FALSE,loadings.label.repel=T, loadings.label.size = 3,loadings.colour = 'black', loadings.label.colour = 'black') + theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"), legend.key.size = unit(.025, "cm"), legend.key.width = unit(.25,"cm"))+ scale_color_colorblind(labels = c(expression(italic("S. cerevisiae")), expression(italic("S. uvarum")))) #ggsave("output/pca_plots/pca_full_Strainv2.jpeg") #Colored by Temperature #temperature <- full$temperature all_temp_pca <- autoplot(full_pca, data = full, colour = 'temperature', scale = T, size = 1.25, loadings = FALSE, loadings.label = FALSE,loadings.label.repel=T, loadings.label.size = 3,loadings.colour = 'black', loadings.label.colour = 'black') + theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"), legend.key.size = unit(.025, "cm"), legend.key.width = unit(.25,"cm")) + scale_color_colorblind(labels = c(expression("15C"), expression("25C"))) #ggsave("output/pca_plots/pca_full_temp.jpeg") ``` ______________________________________ PCA plots of 15C (cold) dataset: ```{r} cold_pcaformat <- cold[6:ncol(cold)] %>% data.frame() row.names(cold_pcaformat) <- cold$sample_ID cold_pca <- prcomp(cold_pcaformat, center = TRUE, scale. = TRUE) #summary(cold_pca) cold$strain #PCA plot of cold temp fermentations colored by strain: strain <- cold$strain cold_strain_pca <- autoplot(cold_pca, data = cold, colour = 'strain', scale = TRUE, size = 1.25, loadings = FALSE, loadings.label = FALSE,loadings.label.repel=T, loadings.label.size = 3,loadings.colour = 'black', loadings.label.colour = 'black') + theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"), legend.key.size = unit(.025, "cm"), legend.key.width = unit(.25,"cm"), legend.text=element_text(size=6.5))+ scale_color_calc() #ggsave("output/pca_plots/pca_15C_strain.jpeg") #make a PCA plot of the cold dataset colored by cerevisiae or uvarum: species <- cold$species cold_species_pca <- autoplot(cold_pca, data = cold, colour = 'species', scale = TRUE, size = 1.25, loadings = FALSE, loadings.label = FALSE,loadings.label.repel=T, loadings.label.size = 3,loadings.colour = 'black', loadings.label.colour = 'black') + theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"), legend.key.size = unit(.025, "cm"), legend.key.width = unit(.25,"cm")) + scale_color_colorblind(labels = c(expression(italic("S. cerevisiae")), expression(italic("S. uvarum")))) #ggsave("output/pca_plots/pca_15C_species.jpeg") ``` labels = c(expression(italic("S. cerevisiae")), expression(italic("S. uvarum"))) ___________________________________ PCA plots of 25C (hot) dataset: ```{r} hot_pcaformat <- hot[6:ncol(hot)] %>% data.frame() row.names(hot_pcaformat) <- hot$sample_ID hot_pca <- prcomp(hot_pcaformat, center = TRUE, scale. = TRUE) #summary(hot_pca) #Hot df colored by strain strain <- hot$strain hot_strain_pca <- autoplot(hot_pca, data = hot, colour = 'strain', scale = TRUE, size = 1.25, loadings = FALSE, loadings.label = FALSE,loadings.label.repel=T, loadings.label.size = 3,loadings.colour = 'black', loadings.label.colour = 'black') + theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"), legend.key.size = unit(.025, "cm"), legend.key.width = unit(.25,"cm"), legend.text=element_text(size=6.5)) + scale_color_calc() #ggsave("output/pca_plots/pca_25C_strain.jpeg") #make PCA plot of hot fermentations colored by strain type species <- hot$type hot_species_pca <-autoplot(hot_pca, data = hot, colour = 'species', scale = TRUE, size = 1.25, loadings = FALSE, loadings.label = FALSE,loadings.label.repel=T, loadings.label.size = 3,loadings.colour = 'black', loadings.label.colour = 'black') + theme(panel.background = element_blank(), panel.border=element_rect(fill=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background=element_blank(), axis.text.x=element_text(colour="black"), axis.text.y=element_text(colour="black"), axis.ticks=element_line(colour="black"), plot.margin=unit(c(1,1,1,1),"line"), legend.key.size = unit(.025, "cm"), legend.key.width = unit(.25,"cm")) + scale_color_colorblind(labels = c(expression(italic("S. cerevisiae")), expression(italic("S. uvarum")))) #ggsave("output/pca_plots/pca_25C_plot_species.jpeg") ``` Arrange figures into a grid and made into a single figure. ```{r} pca_grid<-plot_grid(all_temp_pca, all_species_pca, cold_strain_pca, cold_species_pca, hot_strain_pca, hot_species_pca, ncol = 2, labels =c('A', 'B', 'C', 'D', 'E', 'F'), label_size = 10, align = "hv" ) #ggsave("output/pca_plots/pca_gridV4.jpeg") ``` Analyze the differences between each strain WITHIN a temperature using an ANOVA. - Test for homogeneity of variance - levene's test - ANOVA - Post Hoc analysis - Tukey's HSD _____________________________________________________ 25C fermentations ANOVA and Post Hoc analysis: Levene's test ```{r} filtered_df <- hot #this does the test on each column of filtered_df starting at column 4 levene.results <- apply(filtered_df[,6:ncol(filtered_df)],2,function(x) {leveneTest(x~filtered_df$strain)}) #this makes a table of pvalues levene.plvals<-lapply(levene.results,function(x) {c(x$`Pr(>F)`)}) %>% as.data.frame() #Make an excel file for all the hot temp anova related results. #write.xlsx(levene.plvals, file="output/Spreadsheets/25C_output_res.xlsx", sheetName="LevenesTest", row.names=TRUE) ``` ANOVA ```{r} #define df I am analyzing x<-hot #now create an empty df anova.out <- data.frame(ss_trt = rep(NA, 33), ss_error = rep(NA, 33), f_val = rep(NA, 33), p_val = rep(NA, 33), stringsAsFactors = FALSE) row.names(anova.out) <- colnames(x[,6:ncol(x)]) #Create the for loop (compound ~ yeast strain) within in a single trt for (i in 6:ncol(x)) { lm.res<-lm(x[[i]]~ x$strain) anova.res<-anova(lm.res) anova.out$ss_trt[i-5] <- anova.res$`Sum Sq`[2] anova.out$ss_error[i-5] <- anova.res$`Sum Sq`[1] anova.out$f_val[i-5] <- anova.res$`F value`[1] anova.out$p_val[i-5] <- anova.res$`Pr(>F)`[1] } anova.out #Write the results to an excel file inside a folder for this analysis. #write.xlsx(anova.out, file="output/Spreadsheets/25C_output_res.xlsx", sheetName="ANOVA_results", append = TRUE, row.names=TRUE) ``` For the post hoc analysis I will use a Tukey's HSD test and generate letter codes to indicate significance. Strains that share the same letter are not significantly different. ```{r} #Making a loop that will output the results into a df. x <- hot #Make and empty df to fill hsd_df.out <-data.frame(X2_phenylethanol = rep(NA, 11), benzalcohol = rep(NA, 11), dihydroxy_hexamethoxyflavone_putative = rep(NA, 11), hexanoic_acid = rep(NA, 11), octanoic_acid = rep(NA, 11), X2_phenylethyl_acetate = rep(NA, 11), methionol = rep(NA, 11), X3_methylbutanoic_acid = rep(NA, 11), X2.methylbutanoic_acid = rep(NA, 11), ethyl_decanoate = rep(NA, 11), acetic_acid_putative = rep(NA, 11), ethyl_octanoate = rep(NA, 11), hexanol = rep(NA, 11), hexyl_acetate = rep(NA, 11), ethyl_hexanoate = rep(NA, 11), X3methyl_butanol = rep(NA, 11), X2methyl_butanol = rep(NA, 11), X2methylbutyl_acetate = rep(NA, 11), X2methylpropanol = rep(NA, 11), ethyl_3_methyl_butanoate = rep(NA, 11), ethyl_2_methyl_butanoate = rep(NA, 11), ethyl_butanoate = rep(NA, 11), ethyl_2_methyl_propanoate = rep(NA, 11), ethyl_propanoate = rep(NA, 11), ethyl_acetate = rep(NA, 11), ethyl_formate_putative = rep(NA, 11), acetaldehyde_putative = rep(NA, 11), ethananolamine_putative = rep(NA, 11), ethanol= rep(NA, 11), glycerol= rep(NA, 11), glucose= rep(NA, 11), fructose= rep(NA, 11), aucl = rep(NA, 11), stringsAsFactors = FALSE) row.names(hsd_df.out) <- sort(unique(x$strain)) for (i in 6:ncol(x)) { model <- aov(x[[i]]~strain, data = x) hsd_res <- HSD.test(model, "strain", group = TRUE) groups <- data.frame(hsd_res$groups) groups <- groups[order(row.names(groups)),] hsd_df.out[,i-5] <- groups$groups } hsd_df.out #Save the output into the anova results file: #write.xlsx(hsd_df.out, file="output/Spreadsheets/25C_output_res.xlsx", sheetName="HSD_PostHoc_results", append = TRUE, row.names=TRUE) ``` _____________________________________________________ 15C fermentations ANOVA and post hoc analysis: Levene's test ```{r} filtered_df <- cold #this does the test on each column of filtered_df starting at column 4 levene.results <- apply(filtered_df[,6:ncol(filtered_df)],2,function(x) {leveneTest(x~filtered_df$strain)}) #this makes a table of pvalues levene.plvals<-lapply(levene.results,function(x) {c(x$`Pr(>F)`)}) %>% as.data.frame() levene.plvals #Write an excel file that will keep all of our anova results from cold fermentation analysis: #write.xlsx(levene.plvals, file="output/Spreadsheets/15C_output_res.xlsx", sheetName="LevenesTest", row.names=TRUE) ``` ANOVA for 15C (cold) fermentations ```{r} #define df I am analyzing as x so that I dont have to change all the code for cold fermentaitons x<-cold #now create an empty df anova.out <- data.frame(ss_trt = rep(NA, 33), ss_error = rep(NA, 33), f_val = rep(NA, 33), p_val = rep(NA, 33), stringsAsFactors = FALSE) row.names(anova.out) <- colnames(x[,6:ncol(x)]) #Create the for loop (compound ~ yeast strain) within in a single trt for (i in 6:ncol(x)) { lm.res<-lm(x[[i]]~ x$strain) anova.res<-anova(lm.res) anova.out$ss_trt[i-5] <- anova.res$`Sum Sq`[2] anova.out$ss_error[i-5] <- anova.res$`Sum Sq`[1] anova.out$f_val[i-5] <- anova.res$`F value`[1] anova.out$p_val[i-5] <- anova.res$`Pr(>F)`[1] } anova.out #Write the results to an excel file inside a folder for this analysis. #write.xlsx(anova.out, file="output/Spreadsheets/15C_output_res.xlsx", sheetName="ANOVA_results", append = TRUE, row.names=TRUE) ``` For the post hoc analysis I will use a Tukey's HSD test and generate letter codes to indicate significance. Strains that share the same letter are not significantly different. ```{r} #Making a loop that will output the results into a df. x <- cold #Make and empty df to fill hsd_df.out <-data.frame(X2_phenylethanol = rep(NA, 11), benzalcohol = rep(NA, 11), dihydroxy_hexamethoxyflavone_putative = rep(NA, 11), hexanoic_acid = rep(NA, 11), octanoic_acid = rep(NA, 11), X2_phenylethyl_acetate = rep(NA, 11), methionol = rep(NA, 11), X3_methylbutanoic_acid = rep(NA, 11), X2.methylbutanoic_acid = rep(NA, 11), ethyl_decanoate = rep(NA, 11), acetic_acid_putative = rep(NA, 11), ethyl_octanoate = rep(NA, 11), hexanol = rep(NA, 11), hexyl_acetate = rep(NA, 11), ethyl_hexanoate = rep(NA, 11), X3methyl_butanol = rep(NA, 11), X2methyl_butanol = rep(NA, 11), X2methylbutyl_acetate = rep(NA, 11), X2methylpropanol = rep(NA, 11), ethyl_3_methyl_butanoate = rep(NA, 11), ethyl_2_methyl_butanoate = rep(NA, 11), ethyl_butanoate = rep(NA, 11), ethyl_2_methyl_propanoate = rep(NA, 11), ethyl_propanoate = rep(NA, 11), ethyl_acetate = rep(NA, 11), ethyl_formate_putative = rep(NA, 11), acetaldehyde_putative = rep(NA, 11), ethananolamine_putative = rep(NA, 11), ethanol= rep(NA, 11), glycerol= rep(NA, 11), glucose= rep(NA, 11), fructose= rep(NA, 11), aucl = rep(NA, 11), stringsAsFactors = FALSE) row.names(hsd_df.out) <- sort(unique(x$strain)) for (i in 6:ncol(x)) { model <- aov(x[[i]]~strain, data = x) hsd_res <- HSD.test(model, "strain", group = TRUE) groups <- data.frame(hsd_res$groups) groups <- groups[order(row.names(groups)),] hsd_df.out[,i-5] <- groups$groups } hsd_df.out #write output to the excel file #write.xlsx(hsd_df.out, file="output/Spreadsheets/15C_output_res.xlsx", sheetName="HSD_PostHoc_results", append = TRUE, row.names=TRUE) ``` END OF ANOVA AND POST HOC ANALYSIS FOR BOTH 25C AND 15C FERMENTATIONS ___________________________________________________________________________________________ ___________________________________________________________________________________________ Heatmaps ```{r} library(RColorBrewer) library(ComplexHeatmap) ``` ```{r} #15C all compounds averaged cold_averages<-read_xlsx("data/15C.xlsx", sheet= "averages_for_correlation_hm") %>% data.frame() rnames<- cold_averages$strain cold_averages <- cold_averages[2:ncol(cold_averages)] row.names(cold_averages) <- rnames cold_averages <- as.matrix(cold_averages) #25C all compounds averaged hot_averages<-read_xlsx("data/25C.xlsx", sheet= "averages_for_correlation_hm") %>% data.frame() rnames<- hot_averages$strain hot_averages <- hot_averages[2:ncol(hot_averages)] row.names(hot_averages) <- rnames hot_averages <- as.matrix(hot_averages) #Sclaed by column,since we need to absorb the variation between column. trans<- t(as.matrix(hot_averages)) trans<- t(as.matrix(cold_averages)) x <- trans rm <- rowMeans(x) x <- sweep(x, 1, rm) sx <- apply(x, 1, sd) x <- sweep(x, 1, sx, "/") Heatmap(x, name = "x", col = brewer.pal(11,"RdYlBu"), cluster_rows = FALSE, cluster_columns = FALSE, row_split = rep(c("A", "B", "C", "D", "E" ), c(4,5,7,8,4)), cluster_row_slices = FALSE, cluster_column_slices = FALSE) ```