library(BiocManager) library(shiny) library(dplyr) library(tidyverse) library(ggplot2) library(vegan) library(data.table) library(readr) library(phyloseq) library(broom) library(shinycssloaders) #function to take level_5 file and separate for future use separate_taxa <- function(x, meta) { numSamples <<-ncol(x)-1 #separate taxa column colnames(x)[1] <- "index" level_5_sep <- separate(x,"index", into=c("Kingdom","Phylum","Class","Order","Family"),sep=";") #create a dataset of samples with sample name, abundance, and 1% threshold samples <- x %>% summarize_if(is.numeric,sum,na.rm = TRUE) samples <- pivot_longer(samples,1:ncol(samples),names_to="Sample",values_to="Abundance") samples$threshold<-samples$Abundance*0.01 samples$threshold_rare<-min(samples$Abundance)*0.01 #create a dataset of treatments and how they relate to sample columns treatments <- as.data.frame(unique(meta[,2])) colnames(treatments) <- c("treatment") for (x in 1:nrow(treatments)) { treatments[x,2] <- min(which(meta$treatment==treatments[x,1])+1) treatments[x,3] <- max(which(meta$treatment==treatments[x,1])+1) } colnames(treatments) <- c("treatment","min_col","max_col") #create different datasets based on taxonomy phylum <- cbind(paste(level_5_sep$Kingdom,level_5_sep$Phylum), level_5_sep[,6:ncol(level_5_sep)]) colnames(phylum)[1] <- "Phylum" class <- cbind(paste(level_5_sep$Kingdom,level_5_sep$Phylum, level_5_sep$Class),level_5_sep[,6:ncol(level_5_sep)]) colnames(class)[1] <- "Class" order <- cbind(paste(level_5_sep$Kingdom,level_5_sep$Phylum, level_5_sep$Class,level_5_sep$Order),level_5_sep[,6:ncol(level_5_sep)]) colnames(order)[1] <- "Order" family <- cbind(paste(level_5_sep$Kingdom,level_5_sep$Phylum, level_5_sep$Class,level_5_sep$Order,level_5_sep$Family), level_5_sep[,6:ncol(level_5_sep)]) colnames(family)[1] <- "Family" #create list with separate taxa files, samples file, and treatment file sep_taxa<-list(phylum=phylum,class=class,order=order,family=family,treatments=treatments,samples=samples) return(sep_taxa) } #function to create "other" category and reorder create_other <- function(longdata,rarified) #LongData/Rare,numSamples { #long data with other category otherdata <- longdata otherdata <- otherdata %>% mutate_if(is.factor,as.character) otherrows=nrow(longdata) #create Other category for(i in 1:numSamples) for(j in seq(i,otherrows,numSamples)) if(rarified==TRUE){ if(otherdata$Abundance[j]% mutate_if(is.character,as.factor) otherdata <- as.data.frame(otherdata) #convert to wide format to sort taxa by overall abundance widedata <- otherdata widedata <- widedata %>% pivot_wider(names_from="Sample",values_from="Abundance",values_fill=list(Abundance=0),values_fn=list(Abundance=sum)) #calculate overall abundance and reorder taxa widedata$abundance <- rowSums(widedata[,-1],na.rm=TRUE) widedata[[1]]<-reorder(widedata[[1]],-widedata$abundance) #delete overall abundance and convert back to long widedata <- widedata[,-ncol(widedata)] newdata <- widedata %>% pivot_longer(cols=-1,names_to="Sample",values_to="Abundance") #add treatment data newdata <- merge(newdata,meta_global,by.x="Sample",by.y=1) names(newdata)[ncol(newdata)] <- "treatment" return(newdata) } # makes datasets needed for future analysis create_datasets <- function(sep_taxa,taxa) { #select taxonomic dataset x <- as.data.frame(sep_taxa[[taxa]]) #count the number of samples numSamples=nrow(sep_taxa$samples) #samples in columns ColumnData <- x %>% group_by_at(1) %>% summarize_if(is.numeric, sum, na.rm = TRUE) #reorder column data by overall abundance ColumnData$abundance<-rowSums(ColumnData[,-1],na.rm = TRUE) ColumnData[[1]]<-reorder(ColumnData[[1]],ColumnData$abundance) ColumnData<-ColumnData[,-ncol(ColumnData)] #long data LongData <- ColumnData %>% pivot_longer(cols=-1,names_to="Sample",values_to="Abundance") #long data with other category #long data with other (raw) LongDataOther <- create_other(LongData,FALSE) #create data sets for phyloseq otumat <- data.matrix(subset(ColumnData,select=c(2:ncol(ColumnData)))) rownames(otumat) <- paste0("sp",1:nrow(otumat)) taxmat <- as.matrix(subset(ColumnData,select=c(1))) rownames(taxmat) <- paste0("sp",1:nrow(taxmat)) OTU = otu_table(otumat, taxa_are_rows=TRUE) TAX = tax_table(taxmat) OTU_t <-t(OTU) physeq_samples <- meta_global rownames(physeq_samples)<-meta_global$sample physeq=phyloseq(OTU,TAX,sample_data(as.data.frame(physeq_samples))) #create rarefy datasets physeq2 <- rarefy_even_depth(physeq,rngseed=10,sample.size=min(sample_sums(physeq)),replace=TRUE) OTU_rarefy<-otu_table(physeq2) OTU_rarefy_t<-t(OTU_rarefy) TAX_rarefy <- tax_table(physeq2) #pivot_longer for rarified data ColumnDataRare <- as.data.frame(merge(TAX_rarefy,OTU_rarefy,by="row.names")) #delete row names ColumnDataRare <- ColumnDataRare[-1] #reorder column data by overall abundance ColumnDataRare$abundance<-rowSums(ColumnDataRare[,-1],na.rm = TRUE) ColumnDataRare[[1]]<-reorder(ColumnDataRare[[1]],ColumnDataRare$abundance) ColumnDataRare<-ColumnDataRare[,-ncol(ColumnDataRare)] LongDataRare <- ColumnDataRare %>% pivot_longer(cols=-1,names_to="Sample",values_to="Abundance") #long data with other (rarified) LongDataRareOther <- create_other(LongDataRare,TRUE) #create data set with diversity indices richness<-specnumber(OTU_t) shannon<-diversity(OTU_t,index="shannon") simpson<-diversity(OTU_t,index="simpson") DiversityResults<-cbind(meta_global[,2:ncol(meta_global)],richness,shannon,simpson) #create data set with diversity indices based on rarefy data richness_rarefy<-specnumber(OTU_rarefy_t) shannon_rarefy<-diversity(OTU_rarefy_t,index="shannon") simpson_rarefy<-diversity(OTU_rarefy_t,index="simpson") DiversityResults_rarefy<-cbind(meta_global[,2:ncol(meta_global)],richness_rarefy,shannon_rarefy,simpson_rarefy) #rename column names to be consistent with DiversityResults colnames(DiversityResults_rarefy) <- c("treatment","richness","shannon","simpson") #create list with datasets dataset_list <-list(ColumnData=ColumnData,LongData=LongData,LongDataOther=LongDataOther,OTU=OTU,OTU_t=OTU_t,physeq=physeq, OTU_rarefy=OTU_rarefy,OTU_rarefy_t=OTU_rarefy_t,physeq2=physeq2,ColumnDataRare=ColumnDataRare,LongDataRare=LongDataRare, LongDataRareOther=LongDataRareOther,DiversityResults=DiversityResults,DiversityResults_rarefy=DiversityResults_rarefy,tax=TAX) #return data set list return(dataset_list) } # Function to find core taxa, found in all treatments and samples using ColumnData core_taxa <- function(x) { core <- x %>% filter_if(is.numeric,all_vars(.>0)) core$abundance <- rowSums(core[,-1],na.rm=TRUE) core <- core[order(-core$abundance),] return(core) } # Creates caption/title for core taxa table using metadata, ColumnData, and core make_core_caption <- function(df,core) { tm <- paste(unique(meta_global$treatment),collapse="/") numcore <- nrow(core) total <- nrow(df) caption <- paste("

Core taxa in",tm,"treatments;",numcore,"of",total,"taxa are shared between treatments

") return(toString(caption)) } # Function to find the taxa that are unique to a treatment # with inputs of ColumnData and treatments table and returns it as a list unique_taxa <- function(x,t) { unique_list <- list() for (i in 1:nrow(t)) { unique_sample <- x %>% select(!(t[i,"min_col"]:t[i,"max_col"])) %>% filter_if(is.numeric,all_vars(.==0)) treatment<-as_tibble(t) if(nrow(unique_sample)>0) { unique_list[[i]] <- cbind(treatment[i,1],unique_sample[,1]) #adding the name to the objects in the list is important for creating different tables later name<-treatment[i,1] names(unique_list)[i]<-name } } return(unique_list) } # Create function to graph rarefaction curve with ggplot graph_rare <- function(x,ylab,bytype) { raresample=list() # Write values for each sample from rarecurve function that is # input into a list. Sample names come from samples dataframe for (i in 1:numSamples) { raresample[[i]]<-cbind(samples[i,1],as.data.frame(x[[i]]), attributes(x[[i]])$Subsample) } # Bind the data together from the list raregraph <- bind_rows(raresample) raregraph <-merge(raregraph,meta_global,by.x ="Sample",by.y=1) colnames(raregraph) <- c("Sample","num_taxa","num_samples","Treatment") ggplot(raregraph,aes(x=num_samples,y=num_taxa,group=Sample, color=.data[[bytype]]))+labs(x="Sample Size",y=ylab)+theme_bw()+geom_line(size=1)+guides(fill=guide_legend(title=bytype)) } #creates options for dropdown menus taxachoices <- list("Phylum","Class","Order","Family") raw_or_rare <- list("Raw Data","Rarified Data") abs_or_rel <- list("Absolute abundance","Relative abundance") sample_or_trt <- list("Sample","Treatment") diversitychoices <- list("Richness"="richness","Shannon Diversity"="shannon","Simpson Diversity"="simpson") distancechoices <- list("Jaccard"="jaccard","Bray-Curtis"="bray","Morisita-Horn"="horn") ordinationchoices <- list("NMDS","PCoA") ui <- fluidPage( titlePanel("Bean Beetle Microbiome Analysis"), tabsetPanel(id="tabs", tabPanel(value="tab1",title="Welcome", mainPanel( h2("Welcome to the Bean Beetle Microbiome Analysis App"), p("The Bean Beetle Microbiome Project is a research/teaching collaboration of institutions across the US that is studying the microbiome of", em("Callosobruchus maculatus"), "in the context of course-based undergraduate research experiences (CUREs)."), p("For more information, on the Bean Beetle Microbiome Project, please see", tags$a(href="https://www.beanbeetles.org/microbiome/the-bean-beetle-microbiome-project/", "https://www.beanbeetles.org/microbiome/the-bean-beetle-microbiome-project/")), p("This app is designed to lead students through the community analysis of level-5 (family-level) datasets they have produced in DNA Subway."), p("Before running the app, students should remove any unidentified taxa, chloroplasts, and mitochondria. The level-5 file should be formatted with the first column as the taxa and the subsequent columns as samples with unique sample identifiers. Students should also prepare a metadata file with the first column as samples and the second column as treatments. Both files should be in .csv format."), p("Although this app was designed for community analysis of data from the Bean Beetle Microbiome Project, it could be used for community analysis of any level-5 data."), p("This app was coded by Carolyne Huang, an undergraduate student at Emory University, with assistance from Dr. Chris Beck, Department of Biology, Emory University. The Bean Beetle Microbiome Project is supported by NSF awards DUE-1821533 to Emory University and DUE-1821184 to Morehouse College.") ) ), tabPanel(value = "tab2", title = "Data Upload", sidebarLayout( sidebarPanel( fileInput("file1", "Choose level 5 CSV File", accept = c( "text/csv", "text/comma-separated-values,text/plain", ".csv")), tags$hr(), fileInput("file2", "Choose metadata CSV File", accept = c( "text/csv", "text/comma-separated-values,text/plain", ".csv")), tags$hr(), p("Click Run App after selecting your files and then double-check the format"), actionButton("run", "Run App"), ), mainPanel( h1("Level 5 File"), p("Upload your level 5 file using the interface at the left. The first column should be the combined taxa including kingdom with taxonomic levels separated by semi-colons. The remaining columns should be the abundance data for each sample."), tags$hr(), tableOutput("level_5_contents"), tags$hr(), h1("Metadata file"), p("Upload your metadata file using the interface at the left. The first column should be the samples and the second column should be the treatment for each sample."), tableOutput("metadata_contents") ) ) ), tabPanel(value="tab4",title="Core Taxa", sidebarPanel(selectInput("taxon_core","Select a taxon",taxachoices),width = 2), mainPanel( h2("Core Taxa"), p("Core taxa are those taxa found in all samples."), verticalLayout(htmlOutput("coreCaption"),tableOutput("coreTaxa")), width = 10 ) ), tabPanel(value="tab5", title="Unique Taxa", sidebarPanel(selectInput("taxon_unique","Select a taxon",taxachoices)), mainPanel( h2("Unique Taxa"), p("Unique taxa are those taxa found only in a single treatment."), #this output interface allows for multiple output items textOutput("no_unique_caption"), uiOutput("unique_tables") ) ), tabPanel(value="tab6",title="Rarefaction", sidebarPanel( selectInput("taxon_rarefy","Select a taxonomic level",taxachoices), selectInput("by_type","Graph by",sample_or_trt)), mainPanel( h2("Sample Rarefaction Curves"), h3("Raw Data"), plotOutput("initialRarefaction",width="100%",height="400px")%>% withSpinner(color="#0dc5c1"), h3("Even Rarefaction to minimum number of sequences"), plotOutput("evenRarefaction",width="100%",height="400px")%>% withSpinner(color="#0dc5c1") ) ), tabPanel(value="tab3",title="Taxonomy Bar Graphs",fluid=TRUE, sidebarPanel (selectInput("taxon_bar","Select a taxonomic level",taxachoices), selectInput("rawrare_bar","Select which data to use",raw_or_rare), selectInput("abs_rel","Graph by",abs_or_rel), width = 2), mainPanel(plotOutput("bargraph",width="100%",height = "auto"),width = 10) ), tabPanel(value="tab9",title="Taxonomy Heatmap", sidebarPanel (selectInput("taxon_heatmap","Select a taxonomic level",taxachoices), selectInput("rawrare_heatmap","Select which data to use",raw_or_rare), width = 2), mainPanel(plotOutput("heatmap",height="auto")) ), tabPanel(value="tab7",title="Alpha Diversity", sidebarPanel(selectInput("taxon_alpha_test","Select a taxonomic level",taxachoices), selectInput("rawrare","Select which data to use",raw_or_rare), selectInput("divmeasure","Select diversity measure",diversitychoices)), mainPanel( # h3("Alpha Diversity"), plotOutput("alpha_plots"), tags$hr(), htmlOutput("alpha_caption"), textOutput("anova_caption"), tableOutput("alpha_anova"), textOutput("posthoc_caption"), tableOutput("alpha_posthoc") ) ), tabPanel(value="tab8",title="Beta Diversity", sidebarPanel(selectInput("taxon_beta_test","Select a taxonomic level",taxachoices), selectInput("rawrare_beta","Select which data to use",raw_or_rare), selectInput("distmeasure","Select distance measure",distancechoices), selectInput("ordmethod","Select ordination method",ordinationchoices), selectInput("samptreat","Graph by",sample_or_trt)), mainPanel( h2("Beta Diversity"), plotOutput("ordination_plot"), textOutput("ordination_caption"), tags$hr(), textOutput("permanova_caption"), tableOutput("beta_permanova") ) ) #close tabsetPanel ) #close ui ) server <- function(input, output) { #define level_5 and metadata level_5<-0 metadata<-0 observeEvent(input$tabs,{ if(input$tabs=="tab2"){ #upload level-5 csv file - sidebar for tab2 level_5_upload<-reactive({ inFile1 <- input$file1 if (is.null(inFile1)) return(NULL) read_csv(inFile1$datapath) }) #upload metadata csv file - sidebar for tab2 metadata_upload<-reactive({ inFile2 <- input$file2 if (is.null(inFile2)) return(NULL) read_csv(inFile2$datapath) }) #run uploads and output tables for tab2 observeEvent(input$run, { level_5<-level_5_upload() level_5 <- filter_all(level_5, any_vars(. != 0)) metadata<-metadata_upload() output$level_5_contents <- renderTable({head(level_5)}) output$metadata_contents <- renderTable({metadata}) #make sure that column names in metadata are sample and treatment colnames(metadata)<-c("sample","treatment") sep_taxa <-separate_taxa(level_5,metadata) #create treatments and samples dataset; double arrow makes the lists global samples<<-sep_taxa$samples treatments<<-sep_taxa$treatments meta_global<<-metadata #create datasets at each taxonomic level; double arrow makes the lists global Phylum <<- create_datasets(sep_taxa, "phylum") Class <<- create_datasets(sep_taxa, "class") Order <<- create_datasets(sep_taxa, "order") Family <<- create_datasets(sep_taxa, "family") }) } else if(input$tabs=="tab3"){ if(exists("Phylum")) { #select data and create taxonomy bar graph for tab3 whichBarGraph <- reactive({ taxa <- input$taxon_bar my_list<-get(taxa) if(input$rawrare_bar=="Raw Data"){ my_data <- my_list$LongDataOther } else{ #Rarified Data my_data <- my_list$LongDataRareOther } if(input$abs_rel=="Absolute abundance"){ ggplot(data=my_data,aes(x=Sample,y=Abundance))+ geom_bar(aes(fill=my_data[,2]),position="stack",stat="identity")+ labs(fill=taxa,y="Absolute Abundance")+ facet_grid(.~treatment,space="free_x",scales="free_x")+ theme(legend.position="bottom")+ guides(fill=guide_legend(ncol=2)) } else{ #relative abundance ggplot(data=my_data,aes(x=Sample,y=Abundance))+ geom_bar(aes(fill=my_data[,2]),position="fill",stat="identity")+ labs(fill=taxa,y="Relative Abundance")+ facet_grid(.~treatment,space="free_x",scales="free_x")+ theme(legend.position="bottom")+ guides(fill=guide_legend(ncol=2)) } }) #determine height of taxonomy bar graph for tab3 BarGraphHeight <- reactive({ taxa <- input$taxon_bar my_list<-get(taxa) my_data <- my_list$LongDataOther num_taxa <- length(unique(my_data[,2])) height <-300 + num_taxa*10 return(height) }) #output taxonomy bar graph observe({output$bargraph <- renderPlot({whichBarGraph()},height=BarGraphHeight())}) } } else if(input$tabs=="tab4"){ if(exists("Phylum")) { #creates the core dataset based on taxonomic level selected for tab4 core_data<-reactive({ taxa <- input$taxon_core my_list<-get(taxa) core<-core_taxa(my_list$ColumnData) return(core) }) #creates the caption for the core dataset based on the taxonomic level selected core_data_caption<-reactive({ taxa <- input$taxon_core my_list<-get(taxa) core<-core_taxa(my_list$ColumnData) caption<-make_core_caption(my_list$ColumnData,core) return(caption) }) #returns the caption and core taxa table to the output output$coreCaption<- renderUI({HTML(core_data_caption())}) output$coreTaxa <- renderTable({core_data()},digits=0) } } else if(input$tabs=="tab5"){ if(exists("Phylum")) { #code to create unique taxa output for tab5 #creates the unique taxa dataset based on taxonomic level selected. #It is a list with different objects for each treatment unique_data<-reactive({ taxa <- input$taxon_unique my_list<-get(taxa) unique_taxa_list <- list() unique_taxa_list<-unique_taxa(my_list$ColumnData,treatments) return(unique_taxa_list) }) #creates a list of the number of output tables based on number of treatments #tables are labeled with treatment names output$unique_tables <- renderUI({ table_output_list <- lapply(treatments$treatment, function(i){ tablename <- paste0("table", i) htmlOutput(tablename) }) #this is important, but I'm not sure what it does tagList(table_output_list) }) #renders the tables #the observe command allows a reactive function to be run observe({ unique_data_tables<-unique_data() if(!length(unique_data_tables)) { output$no_unique_caption<-renderText({"There are no unique taxa."}) } for (i in 1:nrow(treatments)) { local({ my_i<-treatments[i,1] #this code creates the caption numunique<-length(unique_data_tables[[my_i]]$treatment) taxa <- input$taxon_unique my_list<-get(taxa) total<-nrow(my_list$ColumnData) caption_title <- paste("Unique taxa in",my_i,"treatment;",numunique,"of",total,"taxa are unique to this treatment") tablename <- paste0("table",my_i) output[[tablename]] <- renderTable( { unique_data_tables[[my_i]] },caption=caption_title,caption.placement = getOption("xtable.caption.placement", "top") ) }) } }) } } else if(input$tabs=="tab6"){ if(exists("Phylum")) { #select data and create initial rarefaction graph for tab6 whichInitialRarefaction <- reactive({ taxa <- input$taxon_rarefy bytype<-input$by_type my_list<-get(taxa) my_data <- my_list$OTU_t graph_rare(rarecurve(my_data,label=FALSE,step=20),taxa,bytype) }) #output initial rarefaction graph output$initialRarefaction <- renderPlot({whichInitialRarefaction()}) #select data and create even rarefaction graph for tab6 whichEvenRarefaction <- reactive({ taxa <- input$taxon_rarefy bytype<-input$by_type my_list<-get(taxa) my_data <- my_list$OTU_rarefy_t graph_rare(rarecurve(my_data,label=FALSE,step=20),taxa,bytype) }) output$evenRarefaction <- renderPlot({whichEvenRarefaction()}) } } else if(input$tabs=="tab7"){ if(exists("Phylum")) { #creates boxplot of alpha diversity depending on taxa, data type, and index whichAlphaPlot <- reactive({ taxa<-input$taxon_alpha_test my_list<-get(taxa) if(input$rawrare=="Raw Data"){ my_data <- my_list$DiversityResults } else{ my_data<- my_list$DiversityResults_rarefy } whichdiv <- input$divmeasure dataMedian <- summarise(group_by(my_data, treatment), MD = round(median(.data[[whichdiv]]),2)) ylabel <-names(diversitychoices)[grep(whichdiv,diversitychoices)] ggplot(my_data,aes(x=treatment,y=.data[[whichdiv]]))+geom_boxplot()+theme_bw()+labs(x="Treatment",y=ylabel)+geom_text(data = dataMedian, aes(treatment, MD, label = MD), position = position_dodge(width = 0.8), size = 3, vjust = -0.5) }) output$alpha_plots <- renderPlot({whichAlphaPlot()}) if(nrow(treatments)==2){ #create caption from results of t-test alpha_caption <- reactive({ taxa <- input$taxon_alpha_test my_list <- get(taxa) if(input$rawrare=="Raw Data"){ my_data <- my_list$DiversityResults } else{ my_data <- my_list$DiversityResults_rarefy } whichdiv <- input$divmeasure result <- t.test(my_data[[whichdiv]]~my_data[,1]) result_t <- round(result$statistic,digits=2) result_df <- round(result$parameter,digits=2) result_p <- round(result$p.value,digits=2) pvalue <- "" if(result_p<0.01){ pvalue <- "p<0.01" } else{ pvalue <- paste("p=",result_p) } #create caption from components of result from t-test caption <- HTML(paste("Welch's Two-sided T-test
t=",result_t,", df=",result_df,",",pvalue)) }) output$alpha_caption<-renderUI({alpha_caption()}) } else{ #create ANOVA table and post-hoc comparisons alpha_anova <- reactive({ taxa <- input$taxon_alpha_test my_list <- get(taxa) if(input$rawrare=="Raw Data"){ my_data <- my_list$DiversityResults } else{ my_data <- my_list$DiversityResults_rarefy } anova_result <- 0 posthoc_result <- 0 anova_tables <- list() whichdiv <- input$divmeasure anova_result <- anova(aov(my_data[[whichdiv]]~my_data[,1])) posthoc_result <- tidy(TukeyHSD(aov(my_data[[whichdiv]]~my_data[,1]))) posthoc_result <- set(posthoc_result,,1,NULL) rownames(anova_result) <- c("Treatment","Residuals") anova_caption <-"ANOVA Table" posthoc_caption <-"Tukey's HSD" anova_tables <- list(anova_result,anova_caption,posthoc_result,posthoc_caption) #return list return(anova_tables) }) #output ANOVA table and caption output$anova_caption <- renderText({alpha_anova()[[2]]}) output$alpha_anova <- renderTable({alpha_anova()[[1]]},rownames=TRUE) #output post-hoc table and caption output$posthoc_caption <- renderText({alpha_anova()[[4]]}) output$alpha_posthoc<-renderTable({alpha_anova()[[3]]}) } } #close tab 7 - alpha diversity tab } else if(input$tabs=="tab8"){ if(exists("Phylum")) { #returns physeq object whichPhyseq <- reactive({ taxa <- input$taxon_beta_test my_list <- get(taxa) if(input$rawrare_beta=="Raw Data"){ my_physeq <- my_list$physeq } else{ my_physeq <- my_list$physeq2 } return(my_physeq) }) #makes ordination data whichOrdinationData <- reactive({ my_ord_data <- ordinate(whichPhyseq(),method=input$ordmethod,distance=input$distmeasure) return(my_ord_data) }) #plots ordination whichOrdinationPlot <- reactive({ if(input$samptreat=="Sample"){ plot_ordination(whichPhyseq(),whichOrdinationData(),color="sample")+theme_bw()+coord_fixed() } else if(input$samptreat=="Treatment"){ plot_ordination(whichPhyseq(),whichOrdinationData(),color="treatment")+stat_ellipse(type="t")+theme_bw()+coord_fixed() } }) #makes ordination caption (NMDS stress) whichOrdinationCaption <- reactive({ orddata <- whichOrdinationData() if(input$ordmethod=="NMDS"){ ordcap <- paste("NMDS results: stress=",round(orddata$stress,digits=4)) } #else{ # ordcap <- "PCoA results" #} }) #plots ordination plot in ui output$ordination_plot <- renderPlot({whichOrdinationPlot()}) output$ordination_caption <- renderText({whichOrdinationCaption()}) #returns adonis results whichPermanova <- reactive({ taxa <- input$taxon_beta_test my_list <- get(taxa) if(input$rawrare_beta=="Raw Data"){ perm_method <- my_list$OTU_t } else{ perm_method <- my_list$OTU_rarefy_t } perm_out <- adonis(perm_method~treatment,data=meta_global,method=input$distmeasure) return(as.data.frame(perm_out$aov.tab)) }) output$permanova_caption <- renderText({"PERMANOVA results"}) output$beta_permanova <- renderTable({whichPermanova()},rownames=TRUE) } } else if(input$tabs=="tab9"){ if(exists("Phylum")) { HeatMapHeight <- reactive({ taxa <- input$taxon_heatmap my_list <- get(taxa) my_data <- my_list$LongData num_taxa <- nrow(unique(my_data[,1])) height <-max(c(225,num_taxa*15)) return(height) }) HeatMap <-reactive({ taxa <- input$taxon_heatmap my_list <- get(taxa) if(input$rawrare_heatmap=="Raw Data"){ my_data <- my_list$LongData } else{ #Rarified Data my_data <- my_list$LongDataRare } ggplot(my_data,aes(x=Sample,y=.data[[taxa]],fill=Abundance))+geom_tile(color="gray")+ theme(legend.justification = "top", axis.text.x=element_text(angle = -90, hjust = 0.5))+scale_x_discrete(position = "top") }) observe({output$heatmap<-renderPlot({HeatMap()},height=HeatMapHeight())}) } } #close tab observe event }) } shinyApp(ui, server)