From 9d3608fd983610c2ecd0d7a9ca3597dadffb0947 Mon Sep 17 00:00:00 2001
From: Midoux Cedric <cedric.midoux@irstea.fr>
Date: Tue, 11 Dec 2018 17:44:55 +0100
Subject: [PATCH] PCA

---
 server.R | 250 ++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 ui.R     |   4 +
 2 files changed, 251 insertions(+), 3 deletions(-)

diff --git a/server.R b/server.R
index 555e6eb..28916b8 100644
--- a/server.R
+++ b/server.R
@@ -1,8 +1,9 @@
-options(shiny.maxRequestSize = 30*1024^2)
+options(shiny.maxRequestSize = 30 * 1024 ^ 2)
 
 library(shinydashboard)
 library(dplyr)
 library(glue)
+library(factoextra)
 
 shinyServer
 (function(input, output, session)
@@ -31,7 +32,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(
@@ -222,7 +223,7 @@ shinyServer
   })
   
   output$tableGlom <- DT::renderDataTable(server = FALSE, {
-    Glom <- tax_glom(data16S(), input$glomRank, NArm=FALSE)
+    Glom <- tax_glom(data16S(), input$glomRank, NArm = FALSE)
     taxTableGlom <- Glom %>%
       tax_table() %>%
       as.data.frame(stringsAsFactors = FALSE) %>%
@@ -1321,4 +1322,247 @@ shinyServer
     }
     return(p + theme_bw())
   })
+  
+  output$pcaUI <- renderUI({
+    validate(need(data16S(), ""))
+    box(
+      title = "Setting : " ,
+      width = NULL,
+      status = "primary",
+      checkboxGroupInput(
+        "pcaSetting",
+        label = "PCA setting",
+        choices = list(
+          "Ratio normalization" = "norm",
+          "Square root normalization" = "sqrt",
+          "Center" = "center",
+          "Scale" = "scale"
+        ),
+        selected = c("norm", "sqrt", "center", "scale"),
+        inline = TRUE
+      ),
+      radioButtons(
+        "pcaType",
+        label = "Type of graph : ",
+        choices = list(
+          "Biplot of individuals and variables" = "biplot",
+          "Graph of individuals" = "ind",
+          "Graph of variables" = "var"
+        ),
+        selected = "biplot",
+        inline = TRUE
+      ),
+      checkboxGroupInput(
+        "pcaAxes",
+        label = "Axes : ",
+        choices = seq(10),
+        selected = c(1, 2),
+        inline = TRUE
+      ),
+      textInput("pcaTitle",
+                label = "Title : ",
+                value = "Principal Component Analysis"),
+      h4(strong("Individuals ( = Samples)")),
+      checkboxGroupInput(
+        "pcaGeomInd",
+        label = "Geometry for individuals  : ",
+        choices = c("point", "text"),
+        selected = c("point", "text"),
+        inline = TRUE
+      ),
+      selectInput(
+        "pcaHabillage",
+        label = "Group : ",
+        choices = c("..." = 0, sample_variables(data16S()))
+      ),
+      checkboxInput("pcaEllipse",
+                    label = "Add ellipses",
+                    value = FALSE),
+      h4(strong("Variables ( = OTU)")),
+      checkboxGroupInput(
+        "pcaGeomVar",
+        label = "Geometry for variables : ",
+        choices = c("arrow", "text"),
+        selected = c("arrow", "text"),
+        inline = TRUE
+      ),
+      sliderInput(
+        "pcaSelect",
+        label = "Select top contrib variables : ",
+        min = 1,
+        max = ntaxa(data16S()),
+        value = 50,
+        step = 1
+      ),
+      collapsedBox(verbatimTextOutput("pcaScript"), title = "RCode")
+    )
+  })
+  
+  output$pcaScript <- renderText({
+    script <- c(
+      scriptHead,
+      "# PCA",
+      "library(factoextra)",
+      "m <- as.data.frame(t(otu_table(data)))"
+    )
+    if ("perc" %in% input$pcaSetting) {
+      script <- c(script,
+                  "m <- m / rowSums(m)")
+    }
+    if ("sqrt" %in% input$pcaSetting) {
+      script <- c(script,
+                  "m <- sqrt(m)")
+    }
+    script <- c(
+      script,
+      glue(
+        "pca <- prcomp(m[colSums(m) != 0], center = {\"center\" %in% input$pcaSetting}, scale = {\"scale\" %in% input$pcaSetting})"
+      ),
+      ""
+    )
+    scriptArgs <-
+      c("pca",
+        glue("axes = c({glue_collapse(input$pcaAxes, sep = ', ')})"))
+    if (input$pcaType %in% c("biplot", "ind"))
+    {
+      if (length(input$pcaGeomInd) == 0)
+      {
+        scriptArgs <- c(scriptArgs, "geom.ind = \"\"")
+      } else if (length(input$pcaGeomInd) == 1)
+      {
+        scriptArgs <-
+          c(scriptArgs, glue("geom.ind = \"{input$pcaGeomInd}\""))
+      } else if (length(input$pcaGeomInd) == 2)
+      {
+        scriptArgs <- c(scriptArgs, "geom.ind = c(\"point\", \"text\")")
+      }
+      if (!is.null(checkNull(input$pcaHabillage)))
+      {
+        scriptArgs <- c(
+          scriptArgs,
+          glue(
+            "habillage = get_variable(data, \"{input$pcaHabillage}\")"
+          ),
+          "invisible = \"quali\"",
+          glue("addEllipses = {input$pcaEllipses}")
+        )
+      }
+    }
+    if (input$pcaType %in% c("biplot", "var"))
+    {
+      if (length(input$pcaGeomVar) == 0)
+      {
+        scriptArgs <- c(scriptArgs, "geom.var = \"\"")
+      } else if (length(input$pcaGeomVar) == 1)
+      {
+        scriptArgs <-
+          c(scriptArgs, glue("geom.var = \"{input$pcaGeomVar}\""))
+      } else if (length(input$pcaGeomVar) == 2)
+      {
+        scriptArgs <- c(scriptArgs, "geom.var = c(\"arrow\", \"text\")")
+      }
+      scriptArgs <- c(scriptArgs,
+                      glue("select.var = list(contrib = {input$pcaSelect})"))
+    }
+    if (!is.null(checkNull(input$pcaTitle))) {
+      scriptArgs <- c(scriptArgs, glue("title = \"{input$pcaTitle}\""))
+    }
+    if (input$pcaType == "biplot")
+    {
+      script <- c(script,
+                  glue(
+                    "p <- fviz_pca_biplot({glue_collapse(scriptArgs, sep=', ')})"
+                  ))
+    }
+    if (input$pcaType == "ind")
+    {
+      script <- c(script,
+                  glue(
+                    "p <- fviz_pca_ind({glue_collapse(scriptArgs, sep=', ')})"
+                  ))
+    }
+    if (input$pcaType == "var")
+    {
+      script <- c(script,
+                  glue(
+                    "p <- fviz_pca_var({glue_collapse(scriptArgs, sep=', ')})"
+                  ))
+    }
+    script <- c(script, "", "plot(p + theme_bw())")
+    
+    return(glue_collapse(script, sep = "\n"))
+  })
+  
+  output$pca <- renderPlot({
+    validate(
+      need(data16S(), "Requires an abundance dataset"),
+      need(length(input$pcaAxes) == 2, "Requires two projections axes")
+    )
+    m <- as.data.frame(t(otu_table(data16S())))
+    if ("norm" %in% input$pcaSetting) {
+      m <- m / rowSums(m)
+    }
+    if ("sqrt" %in% input$pcaSetting) {
+      m <- sqrt(m)
+    }
+    pca <-
+      prcomp(
+        m[colSums(m) != 0],
+        center = ("center" %in% input$pcaSetting),
+        scale = ("scale" %in% input$pcaSetting)
+      )
+    
+    if (!is.null(checkNull(input$pcaHabillage)))
+    {
+      h <- get_variable(data16S(), input$pcaHabillage)
+    } else {
+      h <- "none"
+    }
+    
+    if (input$pcaType == "biplot")
+    {
+      p <- fviz_pca_biplot(
+        pca,
+        axes = as.numeric(input$pcaAxes),
+        geom.ind = c(input$pcaGeomInd, ""),
+        geom.var = c(input$pcaGeomVar, ""),
+        habillage = h,
+        invisible = "quali",
+        addEllipses = input$pcaEllipse,
+        title = input$pcaTitle,
+        select.var = list(contrib = input$pcaSelect)
+        #select.ind
+        #labelsize = 4,
+        #pointsize = 2
+      )
+    }
+    if (input$pcaType == "ind")
+    {
+      p <- fviz_pca_ind(
+        pca,
+        axes = as.numeric(input$pcaAxes),
+        geom.ind = c(input$pcaGeomInd, ""),
+        habillage = h,
+        invisible = "quali",
+        addEllipses = input$pcaEllipse,
+        title = input$pcaTitle
+        #select.ind
+        #labelsize = 4,
+        #pointsize = 2
+      )
+    }
+    if (input$pcaType == "var")
+    {
+      p <- fviz_pca_var(
+        pca,
+        axes = as.numeric(input$pcaAxes),
+        geom.var = c(input$pcaGeomVar, ""),
+        invisible = "quali",
+        title = input$pcaTitle,
+        select.var = list(contrib = input$pcaSelect)
+        #labelsize = 4,
+      )
+    }
+    return(p + theme_bw())
+  })
 })
diff --git a/ui.R b/ui.R
index d169f7e..9384d42 100644
--- a/ui.R
+++ b/ui.R
@@ -1,5 +1,6 @@
 library(shinydashboard)
 library(shinycustomloader)
+
 shinyUI(dashboardPage(
   dashboardHeader(title = "Easy16S"),
   dashboardSidebar(
@@ -158,6 +159,9 @@ shinyUI(dashboardPage(
         withLoader(plotOutput("mds", height = 700)),
         uiOutput("mdsUI")
       ),
+      tabPanel("PCA",
+               withLoader(plotOutput("pca", height = 700)),
+               uiOutput("pcaUI")),
       tabPanel(
         "Phylogenetic tree",
         withLoader(plotOutput("tree", height = 700)),
-- 
GitLab