R/sconeReport.R
d9a78035
 #' SCONE Report Browser: Browse Evaluation of Normalization Performance
d415a9d8
 #'
 #' This function opens a shiny application session for visualizing performance
d9a78035
 #' of a variety of normalization schemes.
d415a9d8
 #'
e9aed5e1
 #' @param x a \code{SconeExperiment} object
 #' @param methods character specifying the normalizations to report.
 #' @param qc matrix. QC metrics to be used for QC evaluation report. Required.
c9cd07e1
 #' @param bio factor. A biological condition (variation to be preserved).
 #'   Default NULL.
 #' @param batch factor. A known batch variable (variation to be removed).
 #'   Default NULL.
d415a9d8
 #' @param negcon character. Genes to be used as negative controls for
 #'   evaluation. These genes should be expected not to change according to the
d9a78035
 #'   biological phenomenon of interest. Default empty character.
d415a9d8
 #' @param poscon character. Genes to be used as positive controls for
 #'   evaluation. These genes should be expected to change according to the
d9a78035
 #'   biological phenomenon of interest. Default empty character.
d415a9d8
 #' @param eval_proj function. Projection function for evaluation  (see
 #'   \code{\link{score_matrix}} for details). If NULL, PCA is used for
a91387f8
 #'   projection.
 #' @param eval_proj_args list. List of args passed to projection function as
 #'   eval_proj_args.
d415a9d8
 #'
d9a78035
 #' @importFrom RColorBrewer brewer.pal
a91387f8
 #' @importFrom rARPACK svds
af3da307
 #' @importFrom graphics boxplot
 #' @importFrom utils write.csv
 #' @importFrom stats setNames var
d9a78035
 #' @export
d415a9d8
 #'
e9aed5e1
 #' @return An object that represents the SCONE report app.
d415a9d8
 #'
e9aed5e1
 #' @examples
c9cd07e1
 #' set.seed(101)
e9aed5e1
 #' mat <- matrix(rpois(1000, lambda = 5), ncol=10)
 #' colnames(mat) <- paste("X", 1:ncol(mat), sep="")
53885e45
 #' obj <- SconeExperiment(mat)
e9aed5e1
 #' res <- scone(obj, scaling=list(none=identity, uq=UQ_FN, deseq=DESEQ_FN),
d415a9d8
 #'            evaluate=TRUE, k_ruv=0, k_qc=0, eval_kclust=2,
c9cd07e1
 #'            bpparam = BiocParallel::SerialParam())
e9aed5e1
 #' qc = as.matrix(cbind(colSums(mat),colSums(mat > 0)))
 #' rownames(qc) = colnames(mat)
 #' colnames(qc) = c("NCOUNTS","NGENES")
ec2f8f98
 #' \dontrun{
 #' sconeReport(res,rownames(get_params(res)), qc = qc)
 #' }
d415a9d8
 #'
e9aed5e1
 sconeReport = function(x, methods,
d415a9d8
                        qc,
d9a78035
                        bio = NULL, batch = NULL,
a91387f8
                        poscon = character(), negcon = character(),
                        eval_proj = NULL,
                        eval_proj_args = NULL){
d415a9d8
 
   if (!requireNamespace("shiny", quietly = TRUE)) {
     stop("shiny package needed for sconeReport()")
   }
 
   if (!requireNamespace("DT", quietly = TRUE)) {
     stop("DT package needed for sconeReport()")
   }
 
   if (!requireNamespace("NMF", quietly = TRUE)) {
     stop("NMF package needed for sconeReport()")
   }
 
   if (!requireNamespace("reshape2", quietly = TRUE)) {
     stop("reshape2 package needed for sconeReport()")
   }
 
af3da307
   if (!requireNamespace("plotly", quietly = TRUE)) {
d415a9d8
     stop("plotly package needed for sconeReport()")
   }
 
af3da307
   if (!requireNamespace("visNetwork", quietly = TRUE)) {
d415a9d8
     stop("visNetwork package needed for sconeReport()")
   }
 
af3da307
   if (!requireNamespace("ggplot2", quietly = TRUE)) {
d415a9d8
     stop("ggplot2 package needed for sconeReport()")
   }
 
e9aed5e1
   scone_res = list()
d415a9d8
 
e9aed5e1
   # List of normalized matrices
d415a9d8
   scone_res$normalized_data = lapply(as.list(methods),
c9cd07e1
                                      FUN = function(z){
                                        get_normalized(x,method = z,log=TRUE)
                                      })
e9aed5e1
   names(scone_res$normalized_data) = methods
d415a9d8
 
e9aed5e1
   # Parameter matrix
   scone_res$params = get_params(x)[methods,]
d415a9d8
 
e9aed5e1
   # Merged score matrix
   scone_res$scores = cbind(get_scores(x),get_score_ranks(x))[methods,]
   colnames(scone_res$scores)[ncol(scone_res$scores)] = "mean_score_rank"
d415a9d8
 
 
d9a78035
   ## ----- If NULL classifications, Replace with NA ------
d415a9d8
 
d9a78035
   if(is.null(bio)){
     bio = factor(rep("NA",ncol(scone_res$normalized_data[[1]])))
   }
d415a9d8
 
d9a78035
   if(is.null(batch)){
     batch = factor(rep("NA",ncol(scone_res$normalized_data[[1]])))
   }
d415a9d8
 
d9a78035
   ## ----- Network Data for Visualization -----
d415a9d8
 
d9a78035
   # Matrix nodes in scone_res
d415a9d8
   leaves = rownames(scone_res$params)
 
d9a78035
   # Parents
   lay1 = unique(gsub(",[^,]*$","",leaves))
   lay2 = unique(gsub(",[^,]*$","",lay1))
   lay3 = unique(gsub(",[^,]*$","",lay2))
   lay4 = unique(gsub(",[^,]*$","",lay3))
   all_nodes = c(leaves,lay1,lay2,lay3,lay4)
   rm(list = c("lay1","lay2","lay3","lay4"))
d415a9d8
 
d9a78035
   # Collapse no-op names for all node names
c9cd07e1
   all_nodes = unique(gsub("^,*|,*$","",
                           gsub("\\,+",",",
                                gsub("none|no_uv|no_bio|no_batch","",
                                     all_nodes))))
d415a9d8
 
d9a78035
   # Order by decreasing number of operations (root at the end)
c9cd07e1
   all_nodes = all_nodes[order(paste0(gsub(".+",",",all_nodes),
                                      gsub("[^,]*","",all_nodes)),
                               decreasing = TRUE)]
d415a9d8
 
d9a78035
   # Collapse no-op names for all leaves
c9cd07e1
   leaves = gsub("^,*|,*$","",
                 gsub("\\,+",",",
                      gsub("none|no_uv|no_bio|no_batch","",
                           leaves)))
d415a9d8
 
d9a78035
   # Re-title nodes in scone_res and root (instead of collapsed names)
   titles = all_nodes
   names(titles) =  paste0("_",all_nodes)
   titles[paste0("_",leaves)] = rownames(scone_res$params)
   titles[titles == ""] = "none"
d415a9d8
 
d9a78035
   # Create Nodes
c9cd07e1
   nodes <- data.frame(id = 1:length(all_nodes),
                       title = titles,
                       group = c("NoData",
d415a9d8
                                 "Loaded")[1 +
c9cd07e1
                                             (all_nodes  %in% leaves)])
d415a9d8
 
d9a78035
   # Create Edges
   edges = data.frame()
   for(i in 1:length(all_nodes)){
c9cd07e1
     parents = as.vector(sapply(all_nodes,
                                function(x){
                                  grepl(paste0("^",x),all_nodes[i])
                                }))
d9a78035
     parents = parents & (all_nodes != all_nodes[i])
     if(any(parents)){
       edges = rbind(edges,cbind(min(which(parents )),i))
     }
   }
   colnames(edges)= c("from","to")
d415a9d8
 
d9a78035
   # Re-map indices according to performance
   new_title_list = rownames(scone_res$params) # Ordered by perfomance
c9cd07e1
   new_title_list = c(new_title_list,
                      as.character(
                        nodes$title[!nodes$title %in%
                                      new_title_list])
   )
   old_ids = nodes[paste0("_",
                          gsub("^,*|,*$","",
                               gsub("\\,+",",",
                                    gsub("none|no_uv|no_uv|no_bio|no_batch",
                                         "",new_title_list)))),]$id
d9a78035
   old_nodes = nodes
   old_edges = edges
   for(i in 1:length(old_ids)){
     nodes$id[old_nodes$id == old_ids[i]] = i
     edges[old_edges == old_ids[i]] = i
   }
   nodes = nodes[order(nodes$id),]
d415a9d8
 
d9a78035
   ## ----- Color Variable -----
   cc <- c(brewer.pal(9, "Set1"), brewer.pal(8, "Set2"),brewer.pal(12,"Set3"),
           brewer.pal(8, "Accent"))
d415a9d8
 
d9a78035
   ## ----- UI Definition -----
d415a9d8
   ui <- shiny::fluidPage(
     shiny::titlePanel("SCONE Report Browser"),
 
     shiny::sidebarLayout(
       shiny::sidebarPanel(
         shiny::selectInput("norm_code", label = "Normalization",
c9cd07e1
                     choices = structure(as.list(rownames(scone_res$params)),
                                         names = rownames(scone_res$params)),
d9a78035
                     selected = rownames(scone_res$params)[1]),
d415a9d8
         shiny::selectInput("color_code", label = "Stratify Plots by",
c9cd07e1
                     choices = list("Batch"="batch",
d415a9d8
                                    "Biological Condition"="bio",
c9cd07e1
                                    "PAM Cluster"="clust"),
d9a78035
                     selected = "batch"),
d415a9d8
         shiny::sliderInput("dim", label = "Reduced Dimension (PCA)",
d9a78035
                     min=2, max=min(length(batch)-1,10),
                     value = 3),
d415a9d8
         shiny::sliderInput("k", label = "Number of Clusters (PAM)",
d9a78035
                     min=2, max=10,
                     value = 2),
d415a9d8
         shiny::helpText(paste0("Please note that with many cells (>500),",
c9cd07e1
                         " some plots may take several seconds ",
                         "(sometimes minutes) to appear.")),
d415a9d8
         shiny::downloadLink('downloadData', 'Download')
d9a78035
       ),
d415a9d8
 
       shiny::mainPanel(shiny::tabsetPanel(type = "tabs",
                             shiny::tabPanel("Overview",
                                      shiny::br(),
c9cd07e1
                                      visNetworkOutput("norm_net",
                                                       width = "650px"),
d415a9d8
                                      shiny::br(),
d9a78035
                                      DT::dataTableOutput('tbl_score')
                             ),
d415a9d8
                             shiny::tabPanel("PCA",
                                      shiny::br(),
af3da307
                                 shiny::p(paste0("This panel shows principal ",
c9cd07e1
                                               "component analysis ",
                                               "(PCA) results",
                                               " on different subsets ",
                                               "of genes.")),
d415a9d8
                                      shiny::br(),
af3da307
                                   shiny::h6("Variance Explained (All Genes)"),
d415a9d8
                                      shiny::p(paste0("Use this plot",
c9cd07e1
                                               " to decide on the ",
                                               "dimensionality of the reduced",
                                               " spaced used for evaluation.")),
af3da307
                                 shiny::plotOutput("plot_scree",width = "650px",
c9cd07e1
                                                 height = "400px"),
d415a9d8
                                      shiny::br(),
                                      shiny::h6("2D (All Genes)"),
                                      shiny::plotOutput("plot_base",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px"),
d415a9d8
                                      shiny::br(),
                                      shiny::h6("3D Interactive (All Genes)"),
c9cd07e1
                                      plotlyOutput("plot3d_base",
                                                   width = "650px",
                                                   height = "650px"),
d415a9d8
                                      shiny::br(),
                                      shiny::selectInput("gene_set",
c9cd07e1
                                                  label = "Gene selection",
d415a9d8
                                                  choices =
c9cd07e1
                                                    list("Negative Controls"=
                                                           "neg",
                                                         "Positive Controls"=
                                                           "pos"),
d9a78035
                                                  selected = "neg"),
d415a9d8
                                      shiny::h6("2D (Select Genes)"),
                                      shiny::p(paste0("Visualize cells by PCs",
c9cd07e1
                                               " of control gene sets.")),
d415a9d8
                                      shiny::plotOutput("plot_select",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px")
d9a78035
                             ),
d415a9d8
                             shiny::tabPanel("QC",
                                      shiny::br(),
af3da307
                               shiny::p(paste0("This panel shows the absolute",
c9cd07e1
                                               " correlations between",
                                               " Principal",
                                               " Components (PCs) of",
                                               " expression",
                                               " matrix vs. PCs of the QC",
                                               " (Quality Control) measure ",
                                               "matrix. It also shows the PCA ",
                                               "of the QC measures. Note that",
                                               " the PCA score signs are not",
                                               " meaningful (e.g.,",
                                               " 'good' cells",
                                               " could have a low value of ",
                                               "PC1 and 'bad' cells a",
                                               " high value).")),
d415a9d8
                                      shiny::br(),
c9cd07e1
                                      plotlyOutput('qccorPlot',
                                                   width = "650px",
                                                   height = "450px"),
d415a9d8
                                      shiny::h5("PCA of QC metrics"),
                                      shiny::h6("2D"),
                                      shiny::plotOutput("plot_qc",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px"),
d415a9d8
                                      shiny::h6("3D Interactive"),
c9cd07e1
                                      plotlyOutput("plot3d_qc",
                                                   width = "650px",
                                                   height = "650px")
d9a78035
                             ),
d415a9d8
                             shiny::tabPanel("Silhouette",
                                      shiny::br(),
                                      shiny::p(paste0("Silhouette Width per",
c9cd07e1
                                               " Sample and Contingency",
                                               " Tables")),
d415a9d8
                                      shiny::br(),
                                      shiny::plotOutput("plotsil",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px"),
d415a9d8
                                      shiny::selectInput("cat1",
c9cd07e1
                                                  label =
                                                    "Row Class",
                                                  choices =
                                                    list("Batch"=
                                                           "batch",
                                                         "Biological Condition"=
                                                           "bio", "PAM Cluster"=
                                                           "clust"),
d9a78035
                                                  selected = "batch"),
d415a9d8
                                      shiny::selectInput("cat2",
c9cd07e1
                                                  label = "Column Class",
                                                  choices =
                                                    list("Batch"=
                                                           "batch",
                                                         "Biological Condition"=
                                                           "bio",
                                                         "PAM Cluster"=
                                                           "clust"),
d9a78035
                                                  selected = "bio"),
d415a9d8
                                      shiny::tableOutput("cat_tab")
d9a78035
                             ),
d415a9d8
                             shiny::tabPanel("Control Genes",
                                      shiny::br(),
af3da307
                                   shiny::p(paste0("Heatmap of control genes, ",
c9cd07e1
                                               "colored by all",
                                               " three categories.")),
d415a9d8
                                      shiny::br(),
                                      shiny::p("Positive controls:"),
                                      shiny::plotOutput("hmposcon",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px"),
d415a9d8
                                      shiny::p("Negative controls:"),
                                      shiny::plotOutput("hmnegcon",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px")
d9a78035
                             ),
d415a9d8
 
                             shiny::tabPanel("Stratified PCA",
                                      shiny::br(),
                                      shiny::p(paste0("Violin plots showing",
c9cd07e1
                                               " conditional ",
                                               "distributions of",
                                               " selected Principal",
                                               " Component.")),
d415a9d8
                                      shiny::br(),
                                      shiny::sliderInput("pcsel",
c9cd07e1
                                                  label =
                                                    paste0("Selected",
                                                           " Principal ",
                                                           "Component"),
                                                  min=1, max=
                                                    min(length(batch)-1,10),
d9a78035
                                                  value = 1),
d415a9d8
                                      shiny::h6("All Genes:"),
                                      shiny::plotOutput("violin_base",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px"),
d415a9d8
                                      shiny::selectInput("gene_set2",
c9cd07e1
                                                  label = "Gene selection",
                                                  choices =
                                                    list("Negative Controls"=
                                                           "neg",
                                                         "Positive Controls"=
                                                           "pos"),
d9a78035
                                                  selected = "neg"),
d415a9d8
                                      shiny::h6("Select Genes:"),
                                      shiny::plotOutput("violin_select",
c9cd07e1
                                                 width = "650px",
                                                 height = "450px")
d9a78035
                             ),
d415a9d8
 
                             shiny::tabPanel("Relative Log-Expression",
                                      shiny::br(),
                                      shiny::p(paste0("Relative ",
c9cd07e1
                                               "Log-Expression Plot ",
                                               "for top 100 Most ",
                                               "Variable Genes.")),
d415a9d8
                                      shiny::plotOutput("rle",
c9cd07e1
                                                 width = "850px",
                                                 height = "450px")
d9a78035
                             )
       ))))
d415a9d8
 
d9a78035
   server <- function(input, output, session) {
d415a9d8
 
d9a78035
     ## ------ Overview Tab ------
d415a9d8
 
d9a78035
     # Render Network
     output$norm_net <- renderVisNetwork({
d415a9d8
 
d9a78035
       # Awk: Check if any non-loaded methods are included
       if(any(!all_nodes  %in% leaves)){
d415a9d8
         visNetwork(nodes, edges,width = "100%",main = "Tree of Methods") %>%
           visHierarchicalLayout(sortMethod = "directed" )  %>%
           visGroups(groupname = "NoData", shape = "dot",
c9cd07e1
                     size = 10, color = list(background = "lightgrey",
                                             border = "darkblue",
d415a9d8
                                             highlight =
c9cd07e1
                                               list(background =
                                                      "black",
d415a9d8
                                                    border = "red")),
                     shadow = list(enabled = TRUE)) %>%
           visGroups(groupname = "Loaded", shape = "dot",
c9cd07e1
                     size = 20, color = list(background = "lightblue",
d415a9d8
                                             border = "darkblue",
                                             highlight =
c9cd07e1
                                               list(background = "cyan",
                                                    border = "red")),
d9a78035
                     shadow = list(enabled = TRUE)) %>%
           visEdges(shadow = TRUE, arrows = list(to = list(enabled = TRUE)),
                    color = list(color = "darkblue", highlight = "red")) %>%
c9cd07e1
           visOptions(nodesIdSelection =
d415a9d8
                        list(enabled = TRUE,
c9cd07e1
                             values = nodes$id[nodes$group == "Loaded"],
                             selected =
                               nodes$id[nodes$title ==
                                          rownames(scone_res$params)[1]])) %>%
d9a78035
           visLegend(width = 0.1, position = "right", main = "Status")
       }else{
d415a9d8
         visNetwork(nodes, edges,width = "100%",main = "Tree of Methods") %>%
           visHierarchicalLayout(sortMethod = "directed" )  %>%
d9a78035
           visEdges(shadow = TRUE, arrows = list(to = list(enabled = TRUE)),
                    color = list(color = "darkblue", highlight = "red")) %>%
c9cd07e1
           visOptions(nodesIdSelection =
d415a9d8
                        list(enabled = TRUE,
c9cd07e1
                             values = nodes$id[nodes$group == "Loaded"],
d415a9d8
                             selected =
c9cd07e1
                               nodes$id[nodes$title ==
                                          rownames(scone_res$params)[1]]))
d9a78035
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
d9a78035
     # Render Table
     output$tbl_score <- DT::renderDataTable({
       DT::datatable(scone_res$scores,
                     extensions = 'Buttons',
                     selection = list( mode = 'single',
                                       target = 'row',
                                       selected = c(1)),
c9cd07e1
                     options=
                       list(columnDefs =
                              list(list(visible=FALSE,
                                        targets=1:(ncol(scone_res$scores)-1))),
d415a9d8
                            dom = 'Bfrtip',
c9cd07e1
                            buttons = list(
                              list(extend = 'colvis',
                                   columns = c(1:(ncol(scone_res$scores)-1))))
                       )) %>%
d9a78035
         DT::formatSignif(columns=c(1:(ncol(scone_res$scores))), digits=3)
     })
d415a9d8
 
d9a78035
     # Update Menu Upon Network Selection of (Loaded) Normalization
d415a9d8
     shiny::observeEvent(input$norm_net_selected,{
       if(is.character(input$norm_net_selected)){ # probably unnecessary
d9a78035
         if(!input$norm_net_selected == ""){ # probably unnecessary
c9cd07e1
           if(as.integer(input$norm_net_selected) %in%
              nodes$id[nodes$group == "Loaded"]){
d415a9d8
             shiny::updateSelectInput(
c9cd07e1
               session,
               "norm_code",
               selected = as.character(
                 nodes$title[nodes$id ==
                               as.integer(input$norm_net_selected)]
               ))
d9a78035
           }
         }
       }
     })
d415a9d8
 
d9a78035
     # Update Menu Upon Table Selection
d415a9d8
     shiny::observeEvent(input$tbl_score_rows_selected,{
       if(length(input$tbl_score_rows_selected) > 0 ){ # probably unnecessary
         shiny::updateSelectInput(session,
                           "norm_code",
c9cd07e1
                           selected = as.character(
                             nodes$title[nodes$id ==
                                           as.integer(
                                             input$tbl_score_rows_selected
                                           )]
                           )
         )
d9a78035
       }
     })
d415a9d8
 
d9a78035
     # Upon Menu Selection of Normalization
d415a9d8
     shiny::observeEvent(input$norm_code,{
 
       # Update Network
d9a78035
       if(is.character(input$norm_net_selected)){
d415a9d8
 
d9a78035
         # If selection is empty, then no valid node is selected.
         if(input$norm_net_selected == ""){
           visNetworkProxy("norm_net") %>%
c9cd07e1
             visSelectNodes(id =
                              array(nodes$id[nodes$title == input$norm_code]))
d9a78035
         }else{
d415a9d8
 
d9a78035
           # If selection is different from menu, then network must be updated.
d415a9d8
           if(as.character(nodes$title[nodes$id ==
c9cd07e1
                                       as.integer(
                                         input$norm_net_selected
                                         )]) != input$norm_code){
d9a78035
             visNetworkProxy("norm_net") %>%
c9cd07e1
               visSelectNodes(id =
                                array(nodes$id[nodes$title == input$norm_code]))
d9a78035
           }
         }
       }
d415a9d8
 
d9a78035
       # Update Table
d415a9d8
 
d9a78035
       # If selection is empty, then no row is selected.
       if(length(input$tbl_score_rows_selected) == 0){
d415a9d8
         DT::dataTableProxy("tbl_score") %>%
           DT::selectRows(nodes$id[nodes$title == input$norm_code])
d9a78035
       }else{
d415a9d8
 
d9a78035
         # If selection is different from menu, then row must be updated.
c9cd07e1
         if(input$tbl_score_rows_selected !=
            nodes$id[nodes$title == input$norm_code]){
d415a9d8
           DT::dataTableProxy("tbl_score") %>%
             DT::selectRows(nodes$id[nodes$title == input$norm_code])
d9a78035
         }
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
d9a78035
     ## ------ PCA Tab ------
d415a9d8
 
     normle <- shiny::reactive({
d9a78035
       as.matrix(scone_res$normalized_data[[input$norm_code]])
     })
d415a9d8
 
     strat_col <- shiny::reactive({
d9a78035
       switch(input$color_code,
              bio = bio,
              batch = batch,
              clust = pam_obj()$clust)
     })
d415a9d8
 
     pc_col <- shiny::reactive({
d9a78035
       cc[strat_col()]
     })
d415a9d8
 
     pc_gene <- shiny::reactive({
d9a78035
       switch(input$gene_set,
              pos = intersect(poscon,rownames(normle())),
              neg = intersect(negcon,rownames(normle())))
     })
d415a9d8
 
     pc_obj_base <- shiny::reactive({
 
a91387f8
       # If user has not provided evaluation projection
       if(is.null(eval_proj)){
d415a9d8
 
a91387f8
         prcomp(t(normle()),
                center = TRUE, scale = TRUE)
d415a9d8
 
a91387f8
       } else {
d415a9d8
 
a91387f8
         proj = eval_proj(normle(),
                          eval_proj_args = eval_proj_args)
         colnames(proj) = paste0("PC",1:ncol(proj))
d415a9d8
 
a91387f8
         pc_out = list(x = proj, sdev = apply(proj,2,sd))
d415a9d8
 
a91387f8
         pc_out
d415a9d8
 
a91387f8
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
     pc_obj_select <- shiny::reactive({
d9a78035
       if(length(pc_gene()) > 0){
d415a9d8
 
a91387f8
         # If user has not provided evaluation projection
         if(is.null(eval_proj)){
d415a9d8
 
a91387f8
           prcomp(t(normle()[pc_gene(),]),
                  center = TRUE, scale = TRUE)
d415a9d8
 
a91387f8
         } else {
d415a9d8
 
a91387f8
           proj = eval_proj(normle()[pc_gene(),],
                            eval_proj_args = eval_proj_args)
           colnames(proj) = paste0("PC",1:ncol(proj))
d415a9d8
 
a91387f8
           pc_out = list(x = proj, sdev = apply(proj,2,sd))
d415a9d8
 
a91387f8
           pc_out
d415a9d8
 
a91387f8
         }
d415a9d8
 
d9a78035
       }else{
         list()
       }
     })
d415a9d8
 
     pam_obj <- shiny::reactive({
d9a78035
       cluster::pam(x = pc_obj_base()$x[,1:input$dim],k = input$k)
     })
d415a9d8
 
     output$plot_scree <- shiny::renderPlot({
c9cd07e1
       plot(pc_obj_base()$sdev[1:min(input$dim*2,
                                     length(batch)-1)]^2/
d415a9d8
              sum(pc_obj_base()$sdev^2),
c9cd07e1
            typ = "l",
d415a9d8
            ylab = "Proportion of Variance",
c9cd07e1
            xlab = "PC Index",
            log = "y")
d9a78035
       abline(v = input$dim, col = "red",lty = 2)
     })
d415a9d8
 
     output$plot_base <- shiny::renderPlot({
d9a78035
       par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
       plot(pc_obj_base()$x[,1:2],col =  pc_col(), pch = 16)
d415a9d8
       legend("topright", inset=c(-0.3,0),
              legend=levels(factor(strat_col())),
c9cd07e1
              fill = cc[sort(unique(factor(strat_col())))])
d9a78035
     })
d415a9d8
 
37faf222
     output$plot3d_base <- plotly::renderPlotly({
d9a78035
       PC1 <- PC2 <- PC3 <- NULL
c9cd07e1
       df <- setNames(data.frame(pc_obj_base()$x[,1:3]),
                      c("PC1", "PC2", "PC3"))
       plot_ly(df, x = ~PC1, y = ~PC2, z = ~PC3,
               type = "scatter3d", mode = "markers",
               marker = list(color=pc_col() ))
d9a78035
     })
d415a9d8
 
     output$plot_select <- shiny::renderPlot({
 
d9a78035
       par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
d415a9d8
 
d9a78035
       if(length(pc_gene()) > 0){
d415a9d8
 
c9cd07e1
         plot(pc_obj_select()$x[,1:2],
              col =  pc_col(), pch = 16)
d415a9d8
         legend("topright", inset=c(-0.3,0),
c9cd07e1
                legend=levels(factor(strat_col())),
                fill = cc[sort(unique(factor(strat_col())))])
d415a9d8
 
d9a78035
       }else{
d415a9d8
 
d9a78035
         plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
         text(0,labels = "Control gene set is empty.")
d415a9d8
 
d9a78035
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
     pam_obj <- shiny::reactive({
d9a78035
       cluster::pam(x = pc_obj_base()$x[,1:input$dim],k = input$k)
     })
d415a9d8
 
d9a78035
     ## ------ QC Tab ------
d415a9d8
 
     cor_qc <- shiny::reactive({
d9a78035
       abs(cor(qc,pc_obj_base()$x))
     })
d415a9d8
 
37faf222
     output$qccorPlot <- plotly::renderPlotly({
d9a78035
       metric <- value <- PC <- NULL
d415a9d8
 
       df = reshape2::melt(cor_qc()[,1:input$dim])
d9a78035
       colnames(df) = c("metric","PC","value")
       df$metric = factor(df$metric,levels = colnames(qc))
c9cd07e1
       df$PC = factor(df$PC,
                      levels = paste0("PC",1:input$dim))
d415a9d8
       p <- ggplot(data = df, aes(x = PC,
                                  fill = metric,
c9cd07e1
                                  weight=value)) +
d415a9d8
         geom_bar(position = "dodge") +
         ylim(0, 1) +
c9cd07e1
         labs(x="PC of Expression",
d415a9d8
              y="Absolute Correlation") +
c9cd07e1
         theme(legend.title=element_blank())
d9a78035
       ggplotly(p)
     })
d415a9d8
 
     pc_obj_qc <- shiny::reactive({
d9a78035
       prcomp(as.matrix(qc),center = TRUE, scale = TRUE)
     })
d415a9d8
 
     output$plot_qc <- shiny::renderPlot({
d9a78035
       par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
       plot(pc_obj_qc()$x[,1:2],col =  pc_col(), pch = 16)
c9cd07e1
       legend("topright", inset=c(-0.3,0),
              legend=levels(factor(strat_col())),
              fill = cc[sort(unique(factor(strat_col())))])
d9a78035
     })
d415a9d8
 
37faf222
     output$plot3d_qc <- plotly::renderPlotly({
e9aed5e1
       if(ncol(pc_obj_qc()$x) >= 3){
         PC1 <- PC2 <- PC3 <- NULL
c9cd07e1
         df <- setNames(data.frame(pc_obj_qc()$x[,1:3]),
                        c("PC1", "PC2", "PC3"))
d415a9d8
         plot_ly(df, x = ~PC1, y = ~PC2, z = ~PC3,
                 type = "scatter3d",
c9cd07e1
                 mode = "markers",
                 marker = list(color=pc_col() ))
e9aed5e1
       }else{
         plot_ly(data.frame(), type = "scatter3d", mode = "markers")
       }
d9a78035
     })
d415a9d8
 
d9a78035
     ## ------ Silhouette Tab ------
d415a9d8
 
     sil_obj <- shiny::reactive({
c9cd07e1
       cluster::silhouette(x = as.numeric(strat_col()),
                           dist = dist(pc_obj_base()$x[,1:input$dim]))
d9a78035
     })
d415a9d8
 
     output$plotsil <- shiny::renderPlot({
 
d9a78035
       par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
d415a9d8
 
d9a78035
       if(length(unique(strat_col())) > 1){
d415a9d8
 
d9a78035
         sil = sil_obj()
         o1 = order(sil[,3])
         o = o1[order(-sil[,1][o1])]
c9cd07e1
         barplot(sil[,3][o], main = paste0("ASW = ",
                                           signif(mean(sil[,3]),3)),
                 xlab = "silhouette width",horiz=TRUE,
d9a78035
                 xlim = 1.5*range(c(sil[,3],-sil[,3])),
                 col = pc_col()[o],
                 border = pc_col()[o])
d415a9d8
         legend("topright", inset=c(-0.3,0),
                legend=levels(factor(strat_col())),
c9cd07e1
                fill = cc[sort(unique(factor(strat_col())))])
d415a9d8
 
d9a78035
       }else{
d415a9d8
 
d9a78035
         plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
         text(0,labels = "Stratify plots by a multi-level classification.")
d415a9d8
 
d9a78035
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
     cat1 <- shiny::reactive({
d9a78035
       switch(input$cat1,
              bio = bio,
              batch = batch,
              clust = pam_obj()$clust)
     })
d415a9d8
 
     cat2 <- shiny::reactive({
d9a78035
       switch(input$cat2,
              bio = bio,
              batch = batch,
              clust = pam_obj()$clust)
     })
d415a9d8
 
     output$cat_tab <- shiny::renderTable({
       table(cat1(), cat2())
d9a78035
     })
d415a9d8
 
 
d9a78035
     ## ------ Control Genes Tab ------
d415a9d8
 
     silo <- shiny::reactive({
d9a78035
       sil = sil_obj()
       o1 = order(sil[,3])
d415a9d8
       o1[order(-sil[,1][o1])]
d9a78035
     })
d415a9d8
 
     output$hmnegcon <- shiny::renderPlot({
 
d9a78035
       if(length(negcon) > 0){
d415a9d8
 
d9a78035
         if(length(unique(strat_col())) > 1){
d415a9d8
 
           NMF::aheatmap(normle()[negcon,], Colv = silo(),
c9cd07e1
                    annCol = list(batch = batch, bio = bio,
                                  pam =  as.factor(pam_obj()$clust)),
d9a78035
                    annColors = list(batch = cc, bio = cc, pam =  cc))
d415a9d8
 
d9a78035
         }else{
d415a9d8
 
d9a78035
           plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
           text(0,labels = "Stratify plots by a multi-level classification.")
d415a9d8
 
d9a78035
         }
d415a9d8
 
d9a78035
       }else{
d415a9d8
 
d9a78035
         plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
         text(0,labels = "Control gene set is empty.")
d415a9d8
 
d9a78035
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
     output$hmposcon <- shiny::renderPlot({
 
d9a78035
       if(length(poscon) > 0){
d415a9d8
 
d9a78035
         if(length(unique(strat_col())) > 1){
d415a9d8
 
           NMF::aheatmap(normle()[poscon,], Colv = silo(),
                    annCol = list(batch = batch,
                                  bio = bio,
c9cd07e1
                                  pam = as.factor(pam_obj()$clust)),
d9a78035
                    annColors = list(batch = cc, bio = cc, pam =  cc))
d415a9d8
 
d9a78035
         }else{
d415a9d8
 
d9a78035
           plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
           text(0,labels = "Stratify plots by a multi-level classification.")
d415a9d8
 
d9a78035
         }
d415a9d8
 
d9a78035
       }else{
d415a9d8
 
d9a78035
         plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
         text(0,labels = "Control gene set is empty.")
d415a9d8
 
d9a78035
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
d9a78035
     ## ------ Stratified PCA Tab ------
d415a9d8
 
     output$violin_base <- shiny::renderPlot({
 
d9a78035
       Class = factor(strat_col())
d415a9d8
       Val = pc_obj_base()$x[,input$pcsel]
 
       ggplot(data.frame(Class,Val ),aes(x = Class,y = Val))   +
         geom_violin(scale = "width", trim = TRUE, aes(fill = Class))+
d9a78035
         labs(x = "Plot Stratification", y = "PC") +
         coord_cartesian(ylim = max(abs(range(Val)))*c(-1.5,1.5))  +
         scale_fill_manual(values=cc[sort(unique(factor(strat_col())))]) +
         geom_point(colour = "black") + guides(fill=FALSE)
d415a9d8
 
d9a78035
     })
d415a9d8
 
     pc_gene2 <- shiny::reactive({
d9a78035
       switch(input$gene_set2,
              pos = intersect(poscon,rownames(normle())),
              neg = intersect(negcon,rownames(normle())))
     })
d415a9d8
 
     pc_obj_select2 <- shiny::reactive({
 
a91387f8
       # If selected genes exist
d9a78035
       if(length(pc_gene2()) > 0){
d415a9d8
 
a91387f8
         # If user has not provided evaluation projection
         if(is.null(eval_proj)){
d415a9d8
 
a91387f8
           prcomp(t(normle()[pc_gene2(),]),
                  center = TRUE, scale = TRUE)
d415a9d8
 
a91387f8
         } else {
d415a9d8
 
a91387f8
           proj = eval_proj(normle()[pc_gene2(),],
                            eval_proj_args = eval_proj_args)
           colnames(proj) = paste0("PC",1:ncol(proj))
d415a9d8
 
a91387f8
           pc_out = list(x = proj, sdev = apply(proj,2,sd))
d415a9d8
 
a91387f8
           pc_out
d415a9d8
 
a91387f8
         }
d415a9d8
 
d9a78035
       }else{
         list()
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
     output$violin_select <- shiny::renderPlot({
 
d9a78035
       if(length(pc_gene2()) > 0){
d415a9d8
 
d9a78035
         Class = factor(strat_col())
d415a9d8
         Val = pc_obj_select2()$x[,input$pcsel]
 
         ggplot(data.frame(Class,Val ),aes(x = Class,y = Val))   +
           geom_violin(scale = "width", trim = TRUE, aes(fill = Class))+
d9a78035
           labs(x = "Plot Stratification", y = "PC") +
           coord_cartesian(ylim = max(abs(range(Val)))*c(-1.5,1.5))  +
           scale_fill_manual(values=cc[sort(unique(factor(strat_col())))]) +
           geom_point(colour = "black") + guides(fill=FALSE)
d415a9d8
 
d9a78035
       }else{
d415a9d8
 
d9a78035
         plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
         text(0,labels = "Control gene set is empty.")
d415a9d8
 
d9a78035
       }
     })
d415a9d8
 
d9a78035
     ## ------ Relative Log-Expression ------
d415a9d8
 
     output$rle <- shiny::renderPlot({
 
d9a78035
       par(mar=c(5.1, 4.1, 4.1, 10.1), xpd=TRUE)
d415a9d8
 
d9a78035
       if(length(unique(strat_col())) > 1){
d415a9d8
 
d9a78035
         vars = apply(exp(normle()),1,var)
         is_var = rank(-vars) <= 100
         median <- apply(normle()[is_var,], 1, median)
         rle <- apply(normle()[is_var,], 2, function(x) x - median)
d415a9d8
         boxplot(rle[,rev(silo())],col = pc_col()[rev(silo())],
d9a78035
                 outline = FALSE, names = rep("",ncol(rle)))
         abline(h=0, lty=2)
c9cd07e1
         legend("topright", inset=c(-0.2,0),
d415a9d8
                legend=levels(factor(strat_col())),
c9cd07e1
                fill = cc[sort(unique(factor(strat_col())))])
d415a9d8
 
d9a78035
       }else{
d415a9d8
 
d9a78035
         plot(0, type = 'n',xlab = "",ylab = "",xaxt = 'n',yaxt = 'n')
         text(0,labels = "Stratify plots by a multi-level classification.")
d415a9d8
 
d9a78035
       }
d415a9d8
 
d9a78035
     })
d415a9d8
 
d9a78035
     ## ----- Download Button -----
d415a9d8
 
     output$downloadData <- shiny::downloadHandler(
d9a78035
       filename = function() {
         paste('scone_out-', Sys.Date(), '.csv', sep='')
       },
       content = function(con) {
c9cd07e1
         datt = cbind(as.character(bio),
                      as.character(batch),
                      as.character(pam_obj()$clust))
d9a78035
         colnames(datt) = c("Class-Bio","Class-Batch","Class-Pam")
         rownames(datt) = colnames(normle())
         datt = rbind(t(datt),normle())
         write.csv(datt, con)
       }
     )
d415a9d8
 
d9a78035
   }
d415a9d8
 
d9a78035
   # Shiny App
d415a9d8
   shiny::shinyApp(ui = ui, server = server)
 
0db504ea
 }