inst/shinyEg/ShinyTestPMD/server.R

In YunhuiQi/TestPMD: Post Selection Inference on Results from Penalized Matrix Decomposition (SCCA only)

#
# This is the server logic of a Shiny web application. You can run the
# application by clicking 'Run App' above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#

library(shiny)
library(TestPMD)
library(PMA)
library(ggplot2)
library(cowplot)
library(plyr)

# Define server logic required to draw a histogram
shinyServer(function(input, output, session) {

   X <- reactive({
    if(input$datamethod == "covid multi-omics"){
      X <- covid$metabolite
    }
    else{
      Xcsv <- input$Xcsv
      ext <- tools::file_ext(Xcsv$datapath)
      req(Xcsv)
      validate(need(ext == "csv", "Please upload a csv file"))
      X <- read.csv(Xcsv$datapath, header = T)
    }
     X
  })
   Y <- reactive({
     if(input$datamethod == "covid multi-omics"){
       Y <- covid$protein
     }
     else{
       Ycsv <- input$Ycsv
       ext <- tools::file_ext(Ycsv$datapath)
       req(Ycsv)
       validate(need(ext == "csv", "Please upload a csv file"))
       Y <- read.csv(Ycsv$datapath, header = T)
     }
     Y
   })
   SG <- reactive({
     if(input$datamethod == "covid multi-omics"){
       SG <- ifelse(covid$meta$icu == 1, "ICU", "Non-ICU")
     }
     else{
       if(input$asksamplegroup == "Yes"){
         samplecsv <- input$samplecsv
         ext <- tools::file_ext(samplecsv$datapath)
         req(samplecsv)
         validate(need(ext == "csv", "Please upload a csv file"))
         SG <- read.csv(samplecsv$datapath, header = T)$Group}
       else{SG <- NULL}
     }
     SG
   })
   XFG <- reactive({
     if(input$datamethod == "covid multi-omics"){
       XFG <- c(rep(paste("Group", 1:3, sep = "-"), rep(37, 3)))
     }
     else{
       if(input$askfeaturegroup == "Yes"){
         Xfeaturecsv <- input$Xfeaturecsv
         ext <- tools::file_ext(Xfeaturecsv$datapath)
         req(Xfeaturecsv)
         validate(need(ext == "csv", "Please upload a csv file"))
         XFG <- read.csv(Xfeaturecsv$datapath, header = T)$Group
       }
       else{XFG <- NULL}
      }
     XFG
     })
   YFG <- reactive({
     if(input$datamethod == "covid multi-omics"){
       YFG <- c(rep(paste("Group", 1:8, sep = "-"), rep(57, 8)), rep("Group-9", 61))
     }
     else{
       if(input$askfeaturegroup == "Yes"){
         Yfeaturecsv <- input$Yfeaturecsv
         ext <- tools::file_ext(Yfeaturecsv$datapath)
         req(Yfeaturecsv)
         validate(need(ext == "csv", "Please upload a csv file"))
         YFG <- read.csv(Yfeaturecsv$datapath, header = T)$Group
         }
       else{YFG <- NULL}
     }
     YFG
   })
   out <- reactive({
     x <- standsdmu(X())
     z <- standsdmu(Y())
     withProgress(message = 'Estimating CCA', value = 0,{
     outCCA <- CCA(x = x, z = z, penaltyx = input$rhox, penaltyz = input$rhoy, typex = "standard", typez = "standard", standardize = F, trace = F, K = 10)
     })
     outCCA
     })
   corpvalue <- reactive({
     x <- X()
     y <- Y()
     x <- standsdmu(x)
     y <- standsdmu(y)
     withProgress(message = 'Testing canonical correlations', value = 0,{
       pvalue <- PPTvalue(X = x, Y = y, l_range = 1:input$valuerange, penalty = "Fixed", rho_x = input$rhox, rho_y = input$rhoy, permutation_no = input$valueperm)
     })
       pvalue
   })
   loadingqvalue <- reactive({
     x <- X()
     y <- Y()
     x <- standsdmu(x)
     y <- standsdmu(y)
     withProgress(message = 'Testing loadings', value = 0,{
       pvalue <- CPTloading(X = x, Y = y, side = input$side, K=input$loadingrange, r = 15, penalty = "Fixed", rho_x = input$rhox, rho_y = input$rhoy, permutation_no = input$loadingperm)
     }
     )

     qvalue <- p.adjust(pvalue, method="fdr")
     qvalue
   })

  output$dataXYeg <- renderTable({
    eg <- log10(covid$protein[1:5,1:3])
    colnames(eg) <- c("feature-1", "feature-2", "feature-3")
    eg
  })
  output$dataSGeg <- renderTable({
    head(data.frame(Sample= 1:nrow(covid$metabolite), Group = ifelse(covid$meta$icu == 1, "ICU", "Non-ICU")), 5)
  })
  output$dataFGeg <- renderTable({
    head(data.frame(Feature= paste("feature",1:ncol(covid$metabolite), sep="-"), Group = c(rep(paste("Group", 1:3, sep = "-"), rep(37, 3)))), 5)
  })

  output$dataSUM <- renderText({
    if(input$datamethod == "covid multi-omics"){
      paste("Using multi-omics data for patients with COVID-19, there are 102 samples, 111 metabolites and 517 proteins,
            two sample groups - ICU or Non-ICU, to facilitate the visualization, groups for metabolites are manually created to be 3,
            groups for metabolites are manually created to be 9. For more details, see
            Overmyer K A, Shishkova E, Miller I J, et al. Large-scale multi-omic analysis of COVID-19 severity[J]. Cell systems, 2021, 12(1): 23-40. e7.")
    }
    else{
      paste("There are", nrow(X()), "samples,", ncol(X()),"features in X; ", ncol(Y()),"features in Y.", sep = " ")
    }
  })
  output$dataSGplot <-renderPlot({
    if(input$datamethod == "covid multi-omics" | input$asksamplegroup == "Yes"){
      ggplot(data = data.frame(table(SG())))+
        geom_bar(aes(x = Var1, y = Freq, fill = Var1), stat = "identity")+
        labs(x = "Groups", y = "Freq", title = "Sample groups", fill = "Groups")+
        theme_linedraw()+
        theme(plot.title = element_text(hjust = 0.5))+
        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
              panel.background = element_blank(), axis.line = element_line(colour = "black"))
      }
    })
  output$dataFGplot <-renderPlot({
    if(input$datamethod == "covid multi-omics" | input$askfeaturegroup == "Yes"){
      p_XFG <- ggplot2::ggplot(data = data.frame(table(XFG())))+
        geom_bar(aes(x = Var1, y = Freq, fill = Var1), stat = "identity")+
        labs(x = "Groups", y = "Freq", title = "Metabolite feature groups", fill = "Groups")+
        theme_linedraw()+
        theme(plot.title = element_text(hjust = 0.5))+
        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
              panel.background = element_blank(), axis.line = element_line(colour = "black"))
      p_YFG <- ggplot2::ggplot(data = data.frame(table(YFG())))+
        geom_bar(aes(x = Var1, y = Freq, fill = Var1), stat = "identity")+
        labs(x = "Groups", y = "Freq", title = "Protein feature groups", fill = "Groups")+
        theme_linedraw()+
        theme(plot.title = element_text(hjust = 0.5))+
        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
              panel.background = element_blank(), axis.line = element_line(colour = "black"))
      plot_grid(p_XFG, p_YFG, nrow = 1)
    }
  })

  ### PMD estimation
  output$PMDcor <- renderPlot({
    x <- X()
    y <- Y()
    x <- standsdmu(x)
    y <- standsdmu(y)
    u <- out()$u
    v <- out()$v
    Xvariates <- as.matrix(x)%*%as.matrix(u)
    Yvariates <- as.matrix(y)%*%as.matrix(v)
    ggplot(data.frame(Component = 1:input$K, Corr = diag(round(cor(Xvariates[,1:input$K], Yvariates[,1:input$K]), 4))))+
      geom_point(aes(x = Component, y = Corr), size = 4)+
      geom_text(aes(x = Component, y = Corr - 0.05, label = paste("corr=",as.character(Corr), sep="")))+
      ylim(0, 1)+
      xlim(0, input$K + 1)+
      labs(title="Canonical correlations estimation")+
      theme_linedraw()+
      theme(plot.title = element_text(hjust = 0.5))+
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            panel.background = element_blank(), axis.line = element_line(colour = "black"))
  })
  output$PMDloading <- renderPlot({
    x <- X()
    y <- Y()
    x <- standsdmu(x)
    y <- standsdmu(y)
    u <- out()$u
    v <- out()$v
    p_Xloading <- ggplot(data.frame(Features = 1:ncol(x), Loadings = u[,as.numeric(input$visKloading)]))+
      geom_point(aes(x = Features, y = Loadings))+
      geom_hline(yintercept = 0)+
      labs(title = paste("Canonical loadings estimation wrt X for K = ", input$visKloading, sep=" "))+
      theme_linedraw()+
      theme(plot.title = element_text(hjust = 0.5))+
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            panel.background = element_blank(), axis.line = element_line(colour = "black"))
    p_Yloading <- ggplot(data.frame(Features = 1:ncol(y), Loadings = v[,as.numeric(input$visKloading)]))+
      geom_point(aes(x = Features, y = Loadings))+
      labs(title = paste("Canonical loadings estimation wrt Y for K = ", input$visKloading, sep=" "))+
      theme_linedraw()+
      geom_hline(yintercept = 0)+
      theme(plot.title = element_text(hjust = 0.5))+
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            panel.background = element_blank(), axis.line = element_line(colour = "black"))
    plot_grid(p_Xloading, p_Yloading, nrow = 2)
  })
  output$PMDscatter <- renderPlot({
    x <- X()
    y <- Y()
    p_list <- plotScatter(X = x, Y = y, CCA_out = out(), groups = SG(), K = as.numeric(input$visKcor) )
    ggpubr::ggarrange(plotlist = p_list, nrow = 1, common.legend = TRUE, legend = "right")
    })
  ## Inference
  output$Infercor <- renderPlot({
    x <- X()
    y <- Y()
    x <- standsdmu(x)
    y <- standsdmu(y)
    u <- out()$u
    v <- out()$v
    Xvariates <- as.matrix(x)%*%as.matrix(u)
    Yvariates <- as.matrix(y)%*%as.matrix(v)
    pvalue <- corpvalue()
    ggplot(data.frame(Component = 1:input$valuerange, Corr = diag(round(cor(Xvariates[,1:input$valuerange], Yvariates[,1:input$valuerange]), 4))))+
      geom_point(aes(x = Component, y = Corr), size = 2)+
      geom_text(aes(x = Component, y = Corr - 0.05, label = paste("corr=",as.character(Corr), sep="")))+
      geom_text(aes(x = Component, y = Corr - 0.1, label = paste("pvalue=",as.character(pvalue), sep = "")))+
      ylim(0, 1)+
      xlim(0, input$valuerange + 1)+
      labs(title="Canonical correlations tests")+
      theme_linedraw()+
      theme(plot.title = element_text(hjust = 0.5))+
      theme_linedraw()+
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                     panel.background = element_blank(), axis.line = element_line(colour = "black"))
  })
  output$Inferloading <- renderPlot({
    x <- X()
    y <- Y()
    x <- standsdmu(x)
    y <- standsdmu(y)
    u <- out()$u
    v <- out()$v
    Xvariates <- as.matrix(x)%*%as.matrix(u)
    Yvariates <- as.matrix(y)%*%as.matrix(v)
    qvalue <- loadingqvalue()
    if(input$side == "X"){
      data = x
      loadings = u[, input$loadingrange]
      feature_groups = XFG()
    }
    else{
      data = y
      loadings = v[, input$loadingrange]
      feature_groups = YFG()
    }
    plotBar(data = data, loadings = loadings, qvalues = qvalue, feature_groups = feature_groups, alpha = input$loadingsig)
  })
  output$Inferfisher <- renderPlot({
    x <- X()
    y <- Y()
    x <- standsdmu(x)
    y <- standsdmu(y)
    u <- out()$u
    v <- out()$v
    qvalue <- loadingqvalue()

    if(input$side == "X"){
      data = x
      loadings = u[, input$loadingrange]
      feature_groups = XFG()
    }
    else{
      data = y
      loadings = v[, input$loadingrange]
      feature_groups = YFG()
    }
    groupFisher(data = data, loadings = loadings, qvalues = qvalue,
                feature_groups = feature_groups, alpha = input$loadingsig, draw = T)
  })
  ### Searching
  output$Searchsource <- renderTable({
    qvalue <- loadingqvalue()
    if(input$searchside == "X"){features <- colnames(X())}
    else{features <- colnames(Y())}
    data.frame(Features = features[order(qvalue)], Qvalues = sort(qvalue), Sig = ifelse(sort(qvalue)<input$loadingsig,"Sig","Insig"))[1:input$dispnum,]
  })
  output$Searchres <- renderTable({
    getInfo(main = input$main, abstract = T)
  })

})