working data import & MDS

WEBS/Webs_NMDSandbetadiversity.R

+library(permute)
+library(vegan)
+##remove previous data objects##
+remove(list=ls())
+#set work directory##
+#setwd("C:\\Users\\Joris\\WEBS_dev")
+setwd("~/WEBS_dev/")
+
+##request citation of this package##
+citation("vegan")
+##read datafile (note that you need to safe your data file in excel first as a csv.file)##
+Alldata <-read.csv("WEBS/dummydataPyr.csv", dec=',', sep=';', header=TRUE, stringsAsFactors = F)
+
+##check dataset##
+head(Alldata)
+names(Alldata)
+
+###select columns with biotic data###
+biotic_allgroups<- (Alldata[5:56])
+
+###select columns with abiotic data and environmental data (=treatment description, abiotics)###
+environmentaldata <- Alldata[2:4]
+
+##calculate dissimilarity matrix, e.g. how different are communities between plots##
+a <- capture.output(q <-metaMDS(biotic_allgroups,"bray", permu=1000, zerodist="add", noshare=0))
+
+##note the stress, this value reflects goodness of fit: stress > 0.3 = poor fit, stress 0.2 = fairly good fit, stress = 0.1 or lower = good fit 
+q$stress
+
+
+#plot communities composition per plot in relation to other plots (the more distant, the more dissimilar)
+a <- plot(q, display="sites")
+
+    #check overlap and difference between treatments
+#!!!Do yourself: replace U and I with your own names of the plot types
+ordihull(q, group=Alldata$Treatment, show="U", col="red", lty=1, cex=.8, label=TRUE)
+ordihull(q, group=Alldata$Treatment, show="I", col ="green", lty=1, cex=.8, label=TRUE)  
+
+##test how environmental factors explain communities
+ef<-envfit(q,environmentaldata,permu=1000)
+plot(ef)
+##check P value (= Pr(>r)) of pH and moisture and how these factors correlate with the NMDS1 axis and NMDS2 axis
+ef
+
+##add taxon names, position taxon reflects their correlation with a given plot type or environmental factor
+stems = colSums(biotic_allgroups)
+taxonnames<-orditorp(q, dis="species", priority = stems, col="grey",pch="", cex=.8)
+
+#calculate beta-diversity: the average distance between plots within a 'treatment'
+dissimilarities <- vegdist(biotic_allgroups, index="bray")
+2dissimilarities
+
+##mean dissimilarity treatment "I" and "U" = beta-diversity I and U
+#I
+mean(dissimilarities[c(1,2,6)])
+sd(dissimilarities[c(1,2,6)])
+
+#U
+mean(dissimilarities[c(13,14,15)])
+sd(dissimilarities[c(13,14,15)])
+
+##species richness per plot is a measure for alpha-diversity
+##!!do this yourself in excel##
+
+
+
+

WEBS/dummydataPyr.csv

+Plot ID;Treatment;pH;moisture;Bact_tot;Actino;Acido;Bact;Fung_tot;Asc;Basi;Fil;Aphe;Acho;Ceph;Monh;Prst;Nem_tot;Plant 1;Plant 2;Plant 3;Plant 4;Plant 5;Plant 6;Plant 7;Plant 8;Plant 9;Plant 10;Plant 11;Plant 12;Plant 13;Plant 14;Plant 15;Plant 16;Plant 17;Plant 18;Plant 19;Plant 20;Insect 1;Insect 2;Insect 3;Insect 4;Insect 5;Insect 6;Insect 7;Insect 8;Insect 9;Insect 10;Insect 11;Insect 12;Insect 13;Insect 14;Insect 15;Insect 16;Insect 17;Insect 18;Insect 19;Insect 20
1;I;6,5;20;200;24;65;1;1000;234;65;43;32;2;300;67;0;900;5;3;0;4;3;0;1;0;4;3;2;0;5;4;0;0;1;1;4;0;1;1;4;4;0;3;3;2;1;0;0;1;1;3;2;5;5;2;3;3
2;I;6,7;21;210;34;54;2;1200;432;66;64;32;2;433;89;0;678;4;3;2;5;4;3;0;5;4;2;3;0;4;1;1;2;5;4;5;1;3;4;1;0;2;2;0;4;0;0;1;5;5;1;4;2;2;0;1;2
3;I;6,4;23;250;43;76;4;1244;343;77;18;12;1;244;87;0;890;2;5;2;4;2;5;3;0;0;2;2;1;2;3;1;0;2;2;1;3;2;4;5;1;2;5;4;3;1;2;4;4;2;3;2;0;5;0;0;3
4;U;7,1;19;500;12;123;10;344;23;43;8;9;6;400;78;123;1000;0;5;0;4;3;4;1;0;2;3;4;4;0;1;3;2;5;3;5;5;5;5;2;0;5;0;3;3;3;3;2;0;2;2;2;0;4;5;5;4
5;U;7,3;23;340;32;121;21;545;100;32;0;7;21;344;78;76;1200;3;0;2;5;4;0;4;3;0;5;2;0;2;0;0;3;1;2;0;0;0;4;2;1;4;4;5;3;5;0;5;3;0;0;3;1;4;3;1;2
6;U;7;21;460;21;150;16;654;212;33;0;8;12;433;67;0;900;5;1;0;0;5;5;1;1;1;3;4;1;4;1;1;3;2;4;0;5;1;3;4;5;5;4;4;0;3;3;5;5;2;4;1;4;5;0;2;3
\ No newline at end of file

Webs2020.Rproj

+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 4
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX

app.R

+#
+# This is a Shiny web application. You can run the application by clicking
+# the 'Run App' button above.
+#
+# Find out more about building applications with Shiny here:
+#
+#    http://shiny.rstudio.com/
+#
+
+library(shiny)
+
+# Define UI for application that draws a histogram
+ui <- fluidPage(
+   
+   # Application title
+   titlePanel("Old Faithful Geyser Data"),
+   
+   # Sidebar with a slider input for number of bins 
+   sidebarLayout(
+      sidebarPanel(
+         sliderInput("bins",
+                     "Number of bins:",
+                     min = 1,
+                     max = 50,
+                     value = 30)
+      ),
+      
+      # Show a plot of the generated distribution
+      mainPanel(
+         plotOutput("distPlot")
+      )
+   )
+)
+
+# Define server logic required to draw a histogram
+server <- function(input, output) {
+   
+   output$distPlot <- renderPlot({
+      # generate bins based on input$bins from ui.R
+      x    <- faithful[, 2] 
+      bins <- seq(min(x), max(x), length.out = input$bins + 1)
+      
+      # draw the histogram with the specified number of bins
+      hist(x, breaks = bins, col = 'darkgray', border = 'white')
+   })
+}
+
+# Run the application 
+shinyApp(ui = ui, server = server)
+

markdown/mds_composition.md

+
+#### Community composition
+
+We have calculated how distant the plots are from each other, now let's see what variable (e.g. treatment) separates them the best.

markdown/mds_stress.md

+
+#### Stress
+
+Stress is an indication for the goodness of fit. As it is a measure between the observed similarity versus the expected similarity, a lower stress value means a better fit
+
+|Stress value|Conclusion|
+|--------------|-----------------------------|
+| 20% or above | Very Poor (not worth doing) |
+| 10%-19.9% | Fair |
+| 5%-9.9% | Good |
+| 2.5%-4.9% | Excellent |
+| 0% - 2.4% | Near Perfect Fit |
\ No newline at end of file

server.R

+library(shiny)
+library(DT)
+library(shinythemes)
+library(magrittr)
+library(shinydashboard)
+library(stringr)
+library(vegan)
+library(reshape2)
+library(ggplot2)
+
+server <- function(input, output, session) {
+    
+
+    student_data <- reactiveVal(NULL)
+    MDS_data <- reactiveVal(NULL)
+    env_fit <- reactiveVal(NULL)
+    
+    
+    observeEvent(input$file, {
+        student_data(import_data())
+        
+
+    })
+    
+    observeEvent(input$sep,{
+        student_data(import_data())
+
+    })
+    
+    observeEvent(student_data(),{
+      MDS_data(run_MDS())
+    })
+
+    observeEvent(input$input_colnames,{
+      tryCatch(
+        {
+        MDS_data(run_MDS())
+        },
+        error = function(cond){
+          print(cond)
+        }
+      )
+    })
+    
+    import_data <- function(){
+        req(input$file)
+        return(read.csv(input$file$datapath, sep = input$sep, stringsAsFactors = F))
+    }
+    
+    run_MDS <- function(out_log = F){
+        selection <- input$input_colnames
+        if(is.null(selection)){
+          return(NULL)
+        }else{
+          biotic_allgroups<- student_data()[selection]
+          
+          log <- capture.output(mds <-metaMDS(biotic_allgroups,"bray", permu=1000, zerodist="add", noshare=0))
+          
+          if(out_log){ return(log) }
+          else{ 
+            return(mds)
+          }  
+        }
+    }
+
+    output$input_dataset <- DT::renderDataTable({
+        req(input$file)
+    
+        return(DT::datatable(student_data()))
+    })
+    
+    output$UIinput_colnames <- renderUI({
+
+
+        data <- student_data()
+        
+        if(!is.null(data)){
+
+            selection_part <- data %>% colnames
+            
+            checkboxGroupInput('input_colnames', 'Deselect everything that is not a biotic variable:',
+                               choices = selection_part, selected = selection_part)
+        } else {
+          checkboxGroupInput('input_colnames', 'Variables:',
+                             choices = c(), selected = c())
+        }
+
+    })
+    
+    output$Variable_select <- renderUI({
+      req(student_data())
+      
+      selection <- input$input_colnames
+      
+      var_names <- student_data() %>% colnames
+      abiotic_allgroups <- var_names[-which(var_names %in% selection)]
+      selectInput('selection_input', label = 'Select label', choices = abiotic_allgroups)
+    })
+    
+    output$Env_factors_select <- renderUI({
+      req(student_data())
+      
+      selection <- input$input_colnames
+      
+      var_names <- student_data() %>% colnames
+      abiotic_allgroups <- var_names[-which(var_names %in% selection)]
+      
+      
+      checkboxGroupInput(inputId = 'environmental_factors', label = 'Select environmental factors',choices = abiotic_allgroups)
+    })
+    
+    output$biotic_groups_show <- renderText({
+        req(student_data())
+        
+        rows <- input$input_colnames
+        return(rows)
+        
+        
+   })
+   
+    output$MDS_stress_plt <- renderPlot({
+        req(student_data())
+        req(MDS_data())
+        out <- tryCatch(
+          {
+            MDS_log <- run_MDS(TRUE)
+            
+            stress_vals <- str_extract_all(MDS_log, regex(pattern = '0[.][0-9]+')) %>% 
+              unlist %>%  
+              as.numeric 
+            
+            
+            ##note the stress, this value reflects goodness of fit: stress > 0.3 = poor fit, stress 0.2 = fairly good fit, stress = 0.1 or lower = good fit 
+            plot(stress_vals, type = 'b', xlab = "Run", ylab = "Stress", ylim = c(0,1), main = "MDS stress")
+            
+            
+            abline(h = 0.2)
+            text(x = 5, y = 0.23, "Upper limmit (> 20%)")
+            
+            abline(h = 0.1)
+            text(x = 5.5, y = 0.13, "Good fit (this and lower)")
+            
+          },
+          error = function(cond){
+            return(plot(0,type='n',axes=FALSE,ann=FALSE))
+          }
+        )
+        
+        
+    })
+    
+    output$MDS_composition_plt <- renderPlot({
+        req(student_data())
+        req(MDS_data())
+        
+        data <- student_data()
+        var_names <- data %>% colnames
+        selection <- input$input_colnames
+        
+        
+        biotic_allgroups<- data[which(var_names %in% selection)]
+        abiotic_allgroups <- data[-which(var_names %in% selection)]
+
+        plot(MDS_data(),  display="sites")
+        
+        for(show in unique(abiotic_allgroups[,input$selection_input])){ 
+          ordihull(MDS_data(), group=abiotic_allgroups[,input$selection_input], 
+                   show=show, col="red", 
+                   lty=1, cex=.8, label=TRUE)    
+        }
+        
+        if(! is.null(input$environmental_factors)){
+          
+          env_data <- student_data()[input$environmental_factors]
+          env_data %<>% apply(., 2, function(x){
+            gsub(',','.',x) %>% as.numeric
+          }) 
+          
+          env_fit(envfit(MDS_data(), env_data , permu=1000))
+          plot(env_fit())  
+        }
+
+        ##add taxon names, position taxon reflects their correlation with a given plot type or environmental factor
+        stems = colSums(biotic_allgroups)
+        taxonnames<-orditorp(MDS_data(), dis="species", priority = stems, col="grey",pch="", cex=.8)
+
+    })
+
+    output$MDS_vector <- renderTable({
+      data <- env_fit()
+      
+      tabl <- cbind(data$vectors$arrows, data$vectors$r, data$vectors$pvals)
+      colnames(tabl) <- c("NMDS1","NMDS2","R",'Pvalue')
+      rownames(tabl) <- rownames(data$vectors$arrows)
+      print(tabl)
+      return(tabl)
+    }, rownames = TRUE, bordered = TRUE)
+}
+

ui.R

+library(shiny)
+library(DT)
+library(shinythemes)
+library(shinydashboard)
+
+sidebar <- dashboardSidebar(
+    sidebarMenu(
+        menuItem("Input Data", tabName = 'InputdataPage'),
+        menuItem("NMDS", tabName = 'NMDSPage')
+    )
+)
+
+
+data_input_tab <- tabItem(tabName = 'InputdataPage',
+                          fluidRow(
+                              column(width=4,
+                              h1("Hi there, welcome!")
+                              ),
+                              sidebarPanel(
+                                  fileInput("file", label = h3("Upload your table")),
+                                  tags$hr(),
+                                  radioButtons(inputId = "sep",
+                                               label = "Separator",
+                                               choices = c(Comma = ',', Semicolon = ';', Tab = '\t'),
+                                               selected = ";")
+                              ), 
+                              mainPanel(
+                                  DT::dataTableOutput("input_dataset")
+
+                              ))
+)
+NMDS_tab <- tabItem(tabName = 'NMDSPage',
+                    fluidRow(
+                        sidebarPanel(
+                            fluidRow(uiOutput('UIinput_colnames'))
+                                         
+
+
+                        ), 
+                        mainPanel(
+                            fluidRow(
+                                column(width = 4,
+                                       shiny::includeMarkdown('markdown/mds_stress.md')
+                                ),
+                                column(width = 6,
+                                       plotOutput("MDS_stress_plt")),
+                            ),
+                            fluidRow(
+                                column(width = 4,
+                                       shiny::includeMarkdown('markdown/mds_composition.md'),
+                                       uiOutput("Variable_select"),
+                                       uiOutput("Env_factors_select")
+                                ),
+                                column(width = 6,
+                                       plotOutput("MDS_composition_plt")
+                                       )
+                            ),
+                            fluidRow(
+                                column(width = 4,
+                                       h4("Significance"),
+                                       h5("Now we need to know if the observed differences for each factor are statistically significant.
+                                          Below you will see for each factor how well it correlates (R) and the significance (Pvalue)"),
+                                       tableOutput("MDS_vector")
+                                ),
+                                column(width = 6
+                                       
+                                )
+                            )
+                        ))
+                        
+)
+
+
+body <- dashboardBody(
+        tabItems(
+            data_input_tab, 
+            NMDS_tab
+
+))
+        
+                
+                
+
+ui <- dashboardPage(
+    dashboardHeader(title = 'Webs2020'),
+    sidebar,
+    body
+)