diff --git a/benchmark.RData b/benchmark.RData
new file mode 100644
index 0000000000000000000000000000000000000000..7a3c300da4149459b57e9ebfdf87b67678cb18ef
Binary files /dev/null and b/benchmark.RData differ
diff --git a/benchmark.png b/benchmark.png
new file mode 100644
index 0000000000000000000000000000000000000000..6eaeebabde61c1edc6f0527969949b8080312224
Binary files /dev/null and b/benchmark.png differ
diff --git a/internals.R b/internals.R
index f64c46fab8a264c070563ef342c402feaf350f24..a39924b27e08464260657a51ccf31468883926c8 100644
--- a/internals.R
+++ b/internals.R
@@ -19,10 +19,9 @@
 
   ## Happy path: excel version
   if (input$CSVsep == "excel") {
-    ## For efficiency and consistency, maybe replace with readr::read_xl
-    sdf <- RcmdrMisc::readXL(input$fileMeta$datapath,
-                             rownames = TRUE,
-                             header = TRUE)
+    sdf <- as.data.frame(readxl::read_excel(input$fileMeta$datapath))
+    row.names(sdf) <- sdf[, 1]
+    sdf <- sdf[, -1]
     sdf$SampleID <- rownames(sdf)
     return(sdf)
   }
diff --git a/perf.R b/perf.R
new file mode 100644
index 0000000000000000000000000000000000000000..b379388f8a656240fc0f8ae2b9986b016ea38238
--- /dev/null
+++ b/perf.R
@@ -0,0 +1,81 @@
+library(microbenchmark)
+
+mb <- microbenchmark::microbenchmark(
+  "clear" = {
+    rm(list = ls())
+    lapply(paste('package:', names(sessionInfo()$otherPkgs), sep=""), detach, character.only=TRUE, unload=TRUE)
+    },
+  "load packages" = {
+    library(shinydashboard)
+    library(glue)
+    source("https://raw.githubusercontent.com/mahendra-mariadassou/phyloseq-extended/master/R/load-extra-functions.R")
+    source("internals.R")
+  },
+  "load demo Chailloux" = {
+    load("demo/demo.RData")
+    data <- get("food")
+  },
+  "barplot" = {
+    p <- plot_bar(physeq = data, fill = "Phylum", x = "Description", title = "OTU abundance barplot")
+    p <- p + facet_grid(". ~ EnvType", scales = "free_x")
+    plot(p)
+  },
+  "filtered plot" = {
+    p <- plot_composition(physeq = data, taxaRank1 = "Kingdom", taxaSet1 = "Bacteria", taxaRank2 = "Phylum", numberOfTaxa = 10, fill = "Phylum", x = "Description")
+    p <- p + facet_grid(". ~ EnvType", scales = "free_x")
+    plot(p)
+  },
+  "heatmap" = {
+    p <- plot_heatmap(prune_taxa(names(sort(taxa_sums(data), decreasing = TRUE)[1:250]), data), distance = "bray", method = "NMDS", low = "yellow", high = "red", na.value = "white", sample.order = "Description", title = "Taxa heatmap by samples")
+    p <- p + facet_grid(". ~ EnvType", scales = "free_x")
+    plot(p)
+  },
+  "alpha" = {
+    p <- plot_richness(physeq = data, measures = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson", "Fisher"), x = "EnvType", color = "EnvType", shape = "FoodType", title = "Alpha diversity graphics")
+    p <- p + geom_boxplot()
+    p <- p + geom_point()
+    plot(p)
+  },
+  "beta" = {
+    beta <- melt(as(distance(data, method = "bray"), "matrix"))
+    colnames(beta) <- c("x", "y", "distance")
+    new_factor = as.factor(get_variable(data, "EnvType"))
+    variable_sort <- as.factor(get_variable(data, "EnvType")[order(new_factor)])
+    L = levels(reorder(sample_names(data), as.numeric(new_factor)))
+    beta$x <- factor(beta$x, levels = L)
+    beta$y <- factor(beta$y, levels = L)
+    palette <- hue_pal()(length(levels(new_factor)))
+    tipColor <- col_factor(palette, levels = levels(new_factor))(variable_sort)
+    p1 <- ggplot(beta, aes(x = x, y = y, fill = distance))
+    p1 <- p1 + geom_tile()
+    p1 <- p1 + ggtitle("Beta diversity heatmap")
+    p1 <- p1 + theme(axis.text.x = element_text(angle = 90, hjust = 1, color = tipColor), axis.text.y = element_text(color = tipColor), axis.title.x = element_blank(), axis.title.y = element_blank())
+    plot(p1 + scale_fill_gradient2())
+  },
+  "rarefaction" = {
+    p <- ggrare(physeq = data, step = 100, se = FALSE, color = "EnvType", label = "Description")
+    p <- p + facet_grid(". ~ FoodType")
+    p <- p + geom_vline(xintercept = min(sample_sums(data)), color = "gray60")
+    p <- p + ggtitle("Rarefaction curves")
+    plot(p)
+  },
+  "acp" = {
+    p <- plot_samples(physeq = data, ordination = ordinate(data, method = "MDS", distance = "unifrac"), axes = c(1, 2), color = "EnvType", shape = "FoodType", replicate = "EnvType", label = "Description", title = "Samples ordination graphic")
+    p <- p + stat_ellipse(aes_string(group = "EnvType"))
+    plot(p + theme_bw())
+  },
+  "tree" = {
+    p <- plot_tree(physeq = prune_taxa(names(sort(taxa_sums(data), decreasing = TRUE)[1:20]), data), method = "sampledodge", color = "EnvType", size = "abundance", label.tips = "taxa_names", sizebase = 5, ladderize = "left", plot.margin = 0, title = "Phylogenetic tree")
+    plot(p)
+  },
+  "clustering" = {
+    p <- plot_clust(physeq = data, dist = "unifrac", method = "ward.D2", color = "EnvType")
+    plot(p)
+  },
+  times = 100, unit = "s", control = list(order="inorder"))
+
+mb
+save(mb, file = "benchmark.RData")
+
+mb_plot <- microbenchmark::autoplot.microbenchmark(mb)
+ggsave("benchmark.png", plot = mb_plot)
diff --git a/server.R b/server.R
index dd72eab327fd3d34db6308dc3dea109e01c467fd..8c666bb98212440539d95f5a9bd7db8cceb57c90 100644
--- a/server.R
+++ b/server.R
@@ -1,4 +1,5 @@
 library(shinydashboard)
+library(glue)
 
 shinyServer
 (function(input, output, session)
@@ -27,7 +28,7 @@ shinyServer
       options = list(
         dom = "lBtip",
         pageLength = 10,
-        lengthMenu = list(c(10, 25, 50, 100,-1), list('10', '25', '50', '100', 'All')),
+        lengthMenu = list(c(10, 25, 50, 100, -1), list('10', '25', '50', '100', 'All')),
         buttons = list(
           'colvis',
           list(
@@ -45,6 +46,17 @@ shinyServer
     )
   }
   
+  collapsedBox <- function(data, title = title)  {
+    box(
+      title = title,
+      width = NULL,
+      status = "primary",
+      collapsible = TRUE,
+      collapsed = TRUE,
+      data
+    )
+  }
+  
   source({
     "https://raw.githubusercontent.com/mahendra-mariadassou/phyloseq-extended/master/R/load-extra-functions.R"
   })
@@ -112,6 +124,20 @@ shinyServer
     })
   }
   
+  scriptHead <- c(
+    "# Loading packages",
+    "source(\"https://raw.githubusercontent.com/mahendra-mariadassou/phyloseq-extended/master/R/load-extra-functions.R\")",
+    "",
+    "# Loading data",
+    glue(
+      "load(\"Easy16S-data.{Sys.Date()}.RData\") # if necessary, adapt the file path"
+    ),
+    "",
+    "# View data",
+    "data",
+    ""
+  )
+  
   output$downloadData <- {
     downloadHandler(
       filename = function() {
@@ -144,12 +170,20 @@ shinyServer
     validate(
       need(
         data16S(),
-        "Firstly, you should select a demo dataset or upload an abundance BIOM file.\nFor example, with Galaxy, a BIOM file can be obtained at the end of FROGS workflow with the 'FROGS BIOM to std BIOM' tool.\n Make sure that the phyloseq object in the RData file is called 'data'."
+        "Firstly, you should select a demo dataset or upload an abundance BIOM file.\nFor example, with Galaxy, a BIOM file can be obtained at the end of FROGS workflow with the 'FROGS BIOM to std BIOM' tool. \nMake sure that the phyloseq object in the RData file is called 'data'."
       )
     )
     data16S()
   })
   
+  output$sampledataTable <- renderUI({
+    validate(need(sample_data(data16S(), errorIfNULL = FALSE), ""))
+    collapsedBox(renderTable({
+      (sapply(sample_data(data16S()), class))
+    }, rownames = TRUE, colnames = FALSE),
+    title = "Class of sample_data")
+  })
+  
   output$summaryTable <- renderUI({
     validate(need(data16S(), ""))
     box(
@@ -220,14 +254,39 @@ shinyServer
         "barGrid",
         label = "Subplot : ",
         choices = c("..." = 0, sample_variables(data16S()))
-      )
-      ,
+      ),
       selectInput(
         "barX",
         label = "X : ",
         choices = c("..." = 0, sample_variables(data16S()))
-      )
+      ),
+      collapsedBox(verbatimTextOutput("histScript"), title = "RCode")
+    )
+  })
+  
+  output$histScript <- renderText({
+    scriptArgs <- c("physeq = data",
+                    glue("fill = \"{input$barFill}\""))
+    if (!is.null(checkNull(input$barX))) {
+      scriptArgs <- c(scriptArgs, glue("x = \"{input$barX}\""))
+    }
+    if (!is.null(checkNull(input$barTitle))) {
+      scriptArgs <- c(scriptArgs, glue("title = \"{input$barTitle}\""))
+    }
+    script <- c(
+      scriptHead,
+      "# Plot barplot",
+      glue("p <- plot_bar({glue_collapse(scriptArgs, sep=', ')})")
     )
+    if (!is.null(checkNull(input$barGrid))) {
+      script <- c(script,
+                  glue(
+                    "p <- p + facet_grid(\". ~ {input$barGrid}\", scales = \"free_x\")"
+                  ))
+    }
+    script <- c(script, "", "plot(p)")
+    
+    return(glue_collapse(script, sep = "\n"))
   })
   
   output$histo <- renderPlot({
@@ -257,7 +316,8 @@ shinyServer
   })
   
   output$histFocusUIfocusTaxa <- renderUI({
-    validate(need(data16S(), ""))
+    validate(need(data16S(), ""),
+             need(input$focusRank, ""))
     selectInput(
       "focusTaxa",
       label = "Selected taxa : ",
@@ -279,24 +339,79 @@ shinyServer
       value = 10
     )
   })
-
+  
   output$histFocusUIfocusGrid <- renderUI({
     validate(need(data16S(), ""))
     selectInput("focusGrid",
                 label = "Subplot : ",
                 choices = c("..." = 0, sample_variables(data16S())))
   })
-
+  
   output$histFocusUIfocusX <- renderUI({
     validate(need(data16S(), ""))
     selectInput("focusX",
                 label = "X : ",
                 choices = c("..." = 0, sample_variables(data16S())))
   })
-
-  output$histoFocus <- renderPlot({
-    validate(need(data16S(),
-                  "Requires an abundance dataset"))
+  
+  output$histFocusUI <- renderUI({
+    validate(need(data16S(), ""))
+    box(
+      title = "Setting : ",
+      width = NULL,
+      status = "primary",
+      uiOutput("histFocusUIfocusRank"),
+      uiOutput("histFocusUIfocusTaxa"),
+      uiOutput("histFocusUIfocusNbTaxa"),
+      uiOutput("histFocusUIfocusGrid"),
+      uiOutput("histFocusUIfocusX"),
+      collapsedBox(verbatimTextOutput("histFocusScript"), title = "RCode")
+    )
+  })
+  
+  output$histFocusScript <- renderText({
+    scriptArgs <- c(
+      "physeq = data",
+      glue("taxaRank1 = \"{input$focusRank}\""),
+      glue("taxaSet1 = \"{input$focusTaxa}\""),
+      glue(
+        "taxaRank2 = \"{rank_names(data16S())[which(rank_names(data16S()) == input$focusRank) + 1]}\""
+      ),
+      glue("numberOfTaxa = {input$focusNbTaxa}"),
+      glue(
+        "fill = \"{rank_names(data16S())[which(rank_names(data16S()) == input$focusRank) + 1]}\""
+      )
+    )
+    if (!is.null(checkNull(input$focusX))) {
+      scriptArgs <- c(scriptArgs, glue("x = \"{input$focusX}\""))
+    }
+    script <- c(
+      scriptHead,
+      "# Plot filtered barplot",
+      glue(
+        "p <- plot_composition({glue_collapse(scriptArgs, sep=', ')})"
+      )
+    )
+    if (!is.null(checkNull(input$focusGrid))) {
+      script <- c(
+        script,
+        glue(
+          "p <- p + facet_grid(\". ~ {input$focusGrid}\", scales = \"free_x\")"
+        )
+      )
+    }
+    script <- c(script, "", "plot(p)")
+    
+    return(glue_collapse(script, sep = "\n"))
+  })
+  
+  output$histFocus <- renderPlot({
+    validate(
+      need(data16S(),
+           "Requires an abundance dataset"),
+      need(input$focusRank, ""),
+      need(input$focusTaxa, "")
+    )
     p <- plot_composition(
       physeq = data16S(),
       taxaRank1 = input$focusRank,
@@ -350,8 +465,28 @@ shinyServer
         "clustCol",
         label = "Color : ",
         choices = c("..." = 0, sample_variables(data16S()))
-      )
+      ),
+      collapsedBox(verbatimTextOutput("clustScript"), title = "RCode")
+    )
+  })
+  
+  output$clustScript <- renderText({
+    scriptArgs <- c(
+      "physeq = data",
+      glue("dist = \"{input$clustDist}\""),
+      glue("method = \"{input$clustMethod}\"")
     )
+    if (!is.null(checkNull(input$clustCol))) {
+      scriptArgs <- c(scriptArgs, glue("color = \"{input$clustCol}\""))
+    }
+    script <- c(
+      scriptHead,
+      "# Plot samples clustering",
+      glue("p <- plot_clust({glue_collapse(scriptArgs, sep=', ')})")
+    )
+    script <- c(script, "", "plot(p)")
+    
+    return(glue_collapse(script, sep = "\n"))
   })
   
   output$clust <- renderPlot({
@@ -422,8 +557,61 @@ shinyServer
         "richnessShape",
         label = "Shape : ",
         choices = c("..." = 0, sample_variables(data16S()))
+      ),
+      collapsedBox(verbatimTextOutput("richnessAScript"), title = "RCode")
+    )
+  })
+  
+  output$richnessAScript <- renderText({
+    if (!is.null(checkNull(input$richnessMeasures))) {
+      measures <-
+        glue("measures = c(\"{glue_collapse(input$richnessMeasures, sep='\", \"')}\")")
+    } else {
+      measures <- NULL
+    }
+    scriptArgs <- c("physeq = data", measures)
+    if (!is.null(checkNull(input$richnessX))) {
+      scriptArgs <- c(scriptArgs, glue("x = \"{input$richnessX}\""))
+    }
+    if (!is.null(checkNull(input$richnessColor))) {
+      scriptArgs <-
+        c(scriptArgs, glue("color = \"{input$richnessColor}\""))
+    }
+    if (!is.null(checkNull(input$richnessShape))) {
+      scriptArgs <-
+        c(scriptArgs, glue("shape = \"{input$richnessShape}\""))
+    }
+    if (!is.null(checkNull(input$richnessTitle))) {
+      scriptArgs <-
+        c(scriptArgs, glue("title = \"{input$richnessTitle}\""))
+    }
+    
+    script <- c(
+      scriptHead,
+      "# Plot boxplot of alpha diversity",
+      glue(
+        "p <- plot_richness({glue_collapse(scriptArgs, sep=', ')})"
       )
     )
+    if (input$richnessBoxplot >= 2) {
+      script <- c(script,
+                  "p <- p + geom_boxplot()")
+    }
+    if (input$richnessBoxplot <= 2) {
+      script <- c(script,
+                  "p <- p + geom_point()")
+    }
+    script <- c(script, "", "plot(p)")
+    script <- c(script, "", "# Tables")
+    script <- c(
+      script,
+      glue(
+        "t <- estimate_richness({glue_collapse(c(\"data\", measures), sep=', ')})"
+      ),
+      "write.table(t, file = \"richness.tsv\", sep = \"\\t\", col.names = NA)"
+    )
+    
+    return(glue_collapse(script, sep = "\n"))
   })
   
   output$richnessA <- renderPlot({
@@ -457,6 +645,7 @@ shinyServer
   })
   
   output$richnessBUI <- renderUI({
+    validate(need(data16S(), ""))
     box(
       title = "Setting : " ,
       width = NULL,
@@ -468,10 +657,66 @@ shinyServer
       ),
       textInput("richnessBTitle",
                 label = "Title : ",
-                value = "Beta diversity heatmap")
+                value = "Beta diversity heatmap"),
+      collapsedBox(verbatimTextOutput("richnessBScript"), title = "RCode")
     )
   })
   
+  output$richnessBScript <- renderText({
+    script <- c(
+      scriptHead,
+      "# Plot heatmap of beta diversity",
+      glue(
+        "beta <- melt(as(distance(data, method = \"{input$richnessBDist}\"), \"matrix\"))"
+      ),
+      "colnames(beta) <- c(\"x\", \"y\", \"distance\")"
+    )
+    if (!is.null(checkNull(input$richnessBOrder))) {
+      script <- c(
+        script,
+        glue(
+          "new_factor = as.factor(get_variable(data, \"{input$richnessBOrder}\"))"
+        ),
+        glue(
+          "variable_sort <- as.factor(get_variable(data, \"{input$richnessBOrder}\")[order(new_factor)])"
+        ),
+        "L = levels(reorder(sample_names(data), as.numeric(new_factor)))",
+        "beta$x <- factor(beta$x, levels = L)",
+        "beta$y <- factor(beta$y, levels = L)",
+        "palette <- hue_pal()(length(levels(new_factor)))",
+        "tipColor <- col_factor(palette, levels = levels(new_factor))(variable_sort)"
+      )
+    } else {
+      script <- c(script, "tipColor <- NULL")
+    }
+    script <-
+      c(
+        script,
+        "",
+        "p1 <- ggplot(beta, aes(x = x, y = y, fill = distance))",
+        "p1 <- p1 + geom_tile()"
+      )
+    if (!is.null(checkNull(input$richnessBTitle))) {
+      script <- c(script,
+                  glue("p1 <- p1 + ggtitle(\"{input$richnessBTitle}\")"))
+    }
+    script <- c(
+      script,
+      glue(
+        "p1 <- p1 + theme(axis.text.x = element_text(angle = 90, hjust = 1, color = tipColor), axis.text.y = element_text(color = tipColor), axis.title.x = element_blank(), axis.title.y = element_blank())"
+      )
+    )
+    script <- c(script, "", "plot(p1 + scale_fill_gradient2())")
+    script <- c(script, "", "# Tables")
+    script <- c(
+      script,
+      glue("t <- distance(data, method = \"{input$richnessBDist}\")"),
+      "write.table(t, file = \"distance.tsv\", sep = \"\\t\", col.names = NA)"
+    )
+    
+    return(glue_collapse(script, sep = "\n"))
+  })
+  
   output$richnessB <- renderPlot({
     validate(need(data16S(),
                   "Requires an abundance dataset"))
@@ -542,8 +787,38 @@ shinyServer
           "Sample name" = "value",
           sample_variables(data16S())
         )
-      )
+      ),
+      collapsedBox(verbatimTextOutput("networkBScript"), title = "RCode")
+    )
+  })
+  
+  output$networkBScript <- renderText({
+    scriptArgs <- c("g",
+                    "physeq = data",
+                    "hjust = 2")
+    if (!is.null(checkNull(input$netwCol))) {
+      scriptArgs <- c(scriptArgs, glue("color = \"{input$netwCol}\""))
+    }
+    if (!is.null(checkNull(input$netwShape))) {
+      scriptArgs <- c(scriptArgs, glue("shape = \"{input$netwShape}\""))
+    }
+    if (!is.null(checkNull(input$netwLabel))) {
+      scriptArgs <- c(scriptArgs, glue("label = \"{input$netwLabel}\""))
+    }
+    if (!is.null(checkNull(input$netwTitle))) {
+      scriptArgs <- c(scriptArgs, glue("title = \"{input$netwTitle}\""))
+    }
+    script <- c(
+      scriptHead,
+      "# Plot samples network",
+      glue(
+        "g <- make_network(data, distance = \"{input$richnessBDist}\", max.dist = {input$netwMax}, keep.isolates = {input$netwOrphan})"
+      ),
+      glue("p <- plot_network({glue_collapse(scriptArgs, sep=', ')})")
     )
+    script <- c(script, "", "plot(p)")
+    
+    return(glue_collapse(script, sep = "\n"))
   })
   
   output$networkB <- renderPlot({
@@ -598,7 +873,6 @@ shinyServer
                        color = "gray60")
     }
     return(p + ggtitle(input$rarefactionTitle))
-    
   })
   
   output$rarefactionCurveUI <- renderUI({
@@ -632,10 +906,49 @@ shinyServer
         "rarefactionGrid",
         label = "Subplot : ",
         choices = c("..." = 0, sample_variables(data16S()))
-      )
+      ),
+      collapsedBox(verbatimTextOutput("rarefactionCurveScript"), title = "RCode")
     )
   })
   
+  output$rarefactionCurveScript <- renderText({
+    scriptArgs <- c("physeq = data",
+                    "step = 100",
+                    "se = FALSE")
+    if (!is.null(checkNull(input$rarefactionColor))) {
+      scriptArgs <-
+        c(scriptArgs,
+          glue("color = \"{input$rarefactionColor}\""))
+    }
+    if (!is.null(checkNull(input$rarefactionLabel))) {
+      scriptArgs <-
+        c(scriptArgs,
+          glue("label = \"{input$rarefactionLabel}\""))
+    }
+    script <- c(
+      scriptHead,
+      "# Plot rarefaction curves",
+      glue("p <- ggrare({glue_collapse(scriptArgs, sep=', ')})")
+    )
+    if (!is.null(checkNull(input$rarefactionGrid))) {
+      script <- c(script,
+                  glue("p <- p + facet_grid(\". ~ {input$rarefactionGrid}\")"))
+    }
+    if (input$rarefactionMin) {
+      script = c(
+        script,
+        "p <- p + geom_vline(xintercept = min(sample_sums(data)), color = \"gray60\")"
+      )
+    }
+    if (!is.null(checkNull(input$rarefactionTitle))) {
+      script <- c(script,
+                  glue("p <- p + ggtitle({input$rarefactionTitle})"))
+    }
+    script <- c(script, "", "plot(p)")
+    
+    return(glue_collapse(script, sep = "\n"))
+  })
+  
   output$HeatmapUI <- renderUI({
     validate(need(data16S(), ""))
     box(
@@ -691,8 +1004,48 @@ shinyServer
           "median",
           "centroid"
         )
+      ),
+      collapsedBox(verbatimTextOutput("heatmapScript"), title = "RCode")
+    )
+  })
+  
+  output$heatmapScript <- renderText({
+    scriptArgs <-
+      c(
+        glue(
+          "prune_taxa(names(sort(taxa_sums(data), decreasing = TRUE)[1:{input$heatmapTopOtu}]), data)"
+        ),
+        glue("distance = \"{input$heatmapDist}\""),
+        glue("method = \"{input$heatmapMethod}\""),
+        "low = \"yellow\"",
+        "high = \"red\"",
+        "na.value = \"white\""
       )
+    if (!is.null(checkNull(input$heatmapX))) {
+      scriptArgs <-
+        c(scriptArgs, glue("sample.order = \"{input$heatmapX}\""))
+    }
+    if (!is.null(checkNull(input$heatmapTitle))) {
+      scriptArgs <-
+        c(scriptArgs, glue("title = \"{input$heatmapTitle}\""))
+    }
+    
+    script <- c(
+      scriptHead,
+      "# Plot heatmap",
+      glue("p <- plot_heatmap({glue_collapse(scriptArgs, sep=', ')})")
     )
+    if (!is.null(checkNull(input$heatmapGrid))) {
+      script <- c(
+        script,
+        glue(
+          "p <- p + facet_grid(\". ~ {input$heatmapGrid}\", scales = \"free_x\")"
+        )
+      )
+    }
+    script <- c(script, "", "plot(p)")
+    
+    return(glue_collapse(script, sep = "\n"))
   })
   
   output$Heatmap <- renderPlot({
@@ -718,7 +1071,7 @@ shinyServer
   })
   
   output$treeUI <- renderUI({
-    validate(need(data16S(), ""))
+    validate(need(phy_tree(data16S(), errorIfNULL = FALSE), ""))
     box(
       title = "Setting : " ,
       width = NULL,
@@ -752,8 +1105,52 @@ shinyServer
         "treeShape",
         label = "Shape : ",
         choices = c("..." = 0, sample_variables(data16S()))
+      ),
+      collapsedBox(verbatimTextOutput("treeScript"), title = "RCode")
+    )
+  })
+  
+  output$treeScript <- renderText({
+    scriptArgs <- c(
+      glue(
+        "physeq = prune_taxa(names(sort(taxa_sums(data), decreasing = TRUE)[1:{input$treeTopOtu}]), data)"
       )
     )
+    if (input$treeSample) {
+      scriptArgs <- c(scriptArgs, "method = \"sampledodge\"")
+    } else {
+      scriptArgs <- c(scriptArgs, "method = \"treeonly\"")
+    }
+    if (!is.null(checkNull(input$treeCol))) {
+      scriptArgs <- c(scriptArgs, glue("color = \"{input$treeCol}\""))
+    }
+    if (!is.null(checkNull(input$treeShape))) {
+      scriptArgs <- c(scriptArgs, glue("shape = \"{input$treeShape}\""))
+    }
+    scriptArgs <- c(scriptArgs, "size = \"abundance\"")
+    if (!is.null(checkNull(input$treeRank))) {
+      scriptArgs <-
+        c(scriptArgs, glue("label.tips = \"{input$treeRank}\""))
+    }
+    scriptArgs <- c(scriptArgs,
+                    "sizebase = 5",
+                    "ladderize = \"left\"",
+                    "plot.margin = 0")
+    if (!is.null(checkNull(input$treeTitle))) {
+      scriptArgs <- c(scriptArgs, glue("title = \"{input$treeTitle}\""))
+    }
+    script <- c(
+      scriptHead,
+      "# Plot phylogenetic tree",
+      glue("p <- plot_tree({glue_collapse(scriptArgs, sep=', ')})")
+    )
+    if (input$treeRadial) {
+      script <- c(script,
+                  "p <- p + coord_polar(theta = \"y\")")
+    }
+    script <- c(script, "", "plot(p)")
+    
+    return(glue_collapse(script, sep = "\n"))
   })
   
   output$tree <- renderPlot({
@@ -845,8 +1242,52 @@ shinyServer
         "acpRep",
         label = "Barycenters : ",
         choices = c("..." = 0, sample_variables(data16S()))
-      )
+      ),
+      collapsedBox(verbatimTextOutput("acpScript"), title = "RCode")
+    )
+  })
+  
+  output$acpScript <- renderText({
+    scriptArgs <- c(
+      "physeq = data",
+      glue(
+        "ordination = ordinate(data, method = \"{input$acpMethod}\", distance = \"{input$acpDist}\")"
+      ),
+      glue("axes = c({glue_collapse(input$acpAxes, sep = ', ')})")
+    )
+    if (!is.null(checkNull(input$acpCol))) {
+      scriptArgs <- c(scriptArgs, glue("color = \"{input$acpCol}\""))
+    }
+    if (!is.null(checkNull(input$acpShape))) {
+      scriptArgs <- c(scriptArgs, glue("shape = \"{input$acpShape}\""))
+    }
+    if (!is.null(checkNull(input$acpRep))) {
+      scriptArgs <- c(scriptArgs, glue("replicate = \"{input$acpRep}\""))
+    } else {
+      scriptArgs <- c(scriptArgs, glue("replicate = NULL"))
+    }
+    if (!is.null(checkNull(input$acpLabel))) {
+      scriptArgs <- c(scriptArgs, glue("label = \"{input$acpLabel}\""))
+    }
+    if (!is.null(checkNull(input$acpTitle))) {
+      scriptArgs <- c(scriptArgs, glue("title = \"{input$acpTitle}\""))
+    }
+    script <- c(
+      scriptHead,
+      "# MultiDimensional scaling",
+      glue("p <- plot_samples({glue_collapse(scriptArgs, sep=', ')})")
     )
+    if (!is.null(checkNull(input$acpEllipse))) {
+      script <- c(
+        script,
+        glue(
+          "p <- p + stat_ellipse(aes_string(group = \"{input$acpEllipse}\"))"
+        )
+      )
+    }
+    script <- c(script, "", "plot(p + theme_bw())")
+    
+    return(glue_collapse(script, sep = "\n"))
   })
   
   output$acp <- renderPlot({
@@ -864,7 +1305,11 @@ shinyServer
       axes = as.numeric(input$acpAxes),
       title = checkNull(input$acpTitle),
       color = checkNull(input$acpCol),
-      replicate = checkNull(input$acpRep),
+      replicate = if (is.null(checkNull(input$acpRep))) {
+        NULL
+      } else {
+        checkNull(input$acpRep)
+      },
       shape = checkNull(input$acpShape),
       label = checkNull(input$acpLabel)
     )
@@ -873,4 +1318,4 @@ shinyServer
     }
     return(p + theme_bw())
   })
-})
\ No newline at end of file
+})
diff --git a/ui.R b/ui.R
index 6d405dec8f10939f4cce4efb43d3ae65a97831b0..94de056697e8f0a02bb575239d8a78b7ed9b9693 100644
--- a/ui.R
+++ b/ui.R
@@ -83,6 +83,7 @@ shinyUI(dashboardPage(
       tabPanel(
         "Summary",
         verbatimTextOutput("phyloseqPrint"),
+        uiOutput("sampledataTable"),
         withLoader(uiOutput("summaryTable")),
         tags$footer(
           "Questions, problems or comments regarding this application should be sent to ",
@@ -103,17 +104,8 @@ shinyUI(dashboardPage(
                uiOutput("histUI")),
       tabPanel(
         "Filtered barplot",
-        withLoader(plotOutput("histoFocus", height = 700)),
-        box(
-          title = "Paramètres",
-          width = NULL,
-          status = "primary",
-          uiOutput("histFocusUIfocusRank"),
-          uiOutput("histFocusUIfocusTaxa"),
-          uiOutput("histFocusUIfocusNbTaxa"),
-          uiOutput("histFocusUIfocusGrid"),
-          uiOutput("histFocusUIfocusX")
-        )
+        withLoader(plotOutput("histFocus", height = 700)),
+        uiOutput("histFocusUI")
       ),
       tabPanel("Heatmap",
                withLoader(plotOutput("Heatmap", height = 700)),
@@ -181,12 +173,12 @@ shinyUI(dashboardPage(
                    Questions, problems or comments regarding this application should be sent to
                    <a href = \"mailto:cedric.midoux@irstea.fr?subject=[Easy16S]\">cedric.midoux@irstea.fr</a>
                    </p>
-
+                   
                    <p>
                    For more information about this tool, you can refer to
                    <a href = \"http://migale.jouy.inra.fr/sites/migale.jouy.inra.fr.drupal7.migale.jouy.inra.fr/files/JOBIM2018_poster.pdf\">this poster</a>.
                    </p>
-
+                   
                    <p>
                    <u>The demo dataset :</u> Chaillou, S., et al. \"
                    <a href = \"https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4409155/\">