An R shiny app that displays both ggplot plots and plotly plots

Question

I'm trying to set up an R shiny app that will enable viewing three types of plots relating to gene expression data.

The data are comprised of:

A data.frame which has the output of the differential expression analysis (each row is a gene and the columns are the effect sizes and their p-values):

set.seed(1)
model.df <- data.frame(id = paste0("g",1:30),symbol = sample(LETTERS[1:5],30,replace=T),
                       group.effect.size = rnorm(30), group.p.value = runif(30,0,1),
                       sex.effect.size = rnorm(30), sex.p.value = runif(30,0,1),
                       stringsAsFactors = F)

A data.frame which has the design of the study (each row is a sample and the columns are the factors that the sample is associated with):

set.seed(1)
design.df <- data.frame(group = c(rep("A",6),rep("B",6)), sex = rep(c(rep("F",3),rep("M",3)),2), replicate = rep(1:6,2)) %>% 
  dplyr::mutate(sample = paste0(group,".",sex,"_",replicate))
design.df$group <- factor(design.df$group, levels = c("A","B"))
design.df$sex <- factor(design.df$sex, levels = c("F","M"))

A matrix which has the abundance (each row is a gene and each column is a sample):

set.seed(1)
abundance.mat <- matrix(rnorm(30*12), nrow=30, ncol=12, dimnames=list(model.df$id,design.df$sample))

A data.frame which has the results of a gene set enrichment analysis (each row is a set name and the columns are the enrichment test p-values for each factor in design.df):

set.seed(1)
gsea.df <- data.frame(set.name = paste0("S",1:4), group.p.value = format(round(runif(4,0,1),2),scientific = T), sex.p.value = format(round(runif(4,0,1),2),scientific = T), stringsAsFactors = F)

And finally, a data.frame which associates the genes with each set.name in gsea.df:

set.seed(1)
gene.sets.df <- do.call(rbind,lapply(1:4,function(s) data.frame(set.name = paste0("S",s), id = sample(model.df$id,10,replace = F),stringsAsFactors = F)))

I would like the shiny app to enable viewing these types of plots:

1. Feature Plot - plotting expression level of a single user-selected gene on the y-axis and sample on the x-axis, and that would be combined with an inset of a caterpillar plot showing the estimated effects:
2. Feature User-Defined Sets Plot - same as Feature Plot, however rather than showing a single -selected gene this will show a set of user-selected-genes and hence rather than points it will show violins of the distributions:
3. Feature Sets GSEA Plot - a combined list of volcano plots, where in each one the x-axis is the effect size of the factor, the y-axis is the -log10(p-value) of the effect, and the genes are colored red if they belong to the selected gene set:

Here are the three functions for generating these figures given the user selection:

featurePlot <- function(selected.id)
{
  replicate.df <- reshape2::melt(abundance.mat[which(rownames(abundance.mat) == selected.id),,drop=F], varnames=c("id","sample")) %>%
    dplyr::left_join(design.df)
  effects.df <- data.frame(factor.name = c("group","sex"), 
                           effect.size = c(dplyr::filter(model.df,id == selected.id)$group.effect.size,dplyr::filter(model.df,id == selected.id)$sex.effect.size),
                           p.value = c(dplyr::filter(model.df,id == selected.id)$group.p.value,dplyr::filter(model.df,id == selected.id)$sex.p.value),
                           stringsAsFactors = F)
  effects.df$factor.name <- factor(effects.df$factor.name, levels = c("group","sex"))
  main.plot <- ggplot(replicate.df,aes(x=replicate,y=value,color=group,shape=sex))+
    geom_point(size=3)+facet_grid(~group,scales="free_x")+
    labs(x="Replicate",y="TPM")+theme_minimal()
  xlims <- c(-1*max(abs(effects.df$effect.size))-0.1*max(abs(effects.df$effect.size)),max(abs(effects.df$effect.size))+0.1*max(abs(effects.df$effect.size)))
  effects.plot <- ggplot(effects.df,aes(x=effect.size,y=factor.name,color=factor.name))+
    geom_point()+
    geom_vline(xintercept=0,linetype="longdash",colour="black",size=0.25)+theme_minimal()+xlim(xlims)+
    theme(legend.position="none")+ylab("")+xlab("Effect Size")

  null.plot <- ggplot(data.frame())+geom_point()+geom_blank()+theme_minimal()
  combined.plot <- gridExtra::arrangeGrob(main.plot,gridExtra::arrangeGrob(null.plot,effects.plot,ncol=1),nrow=1,ncol=2,widths=c(5,2.5))
  return(combined.plot)
}


featureSetPlot <- function(selected.ids)
{
  replicate.df <- reshape2::melt(abundance.mat[which(rownames(abundance.mat) %in% selected.ids),,drop=F], varnames=c("id","sample")) %>%
    dplyr::left_join(design.df)
  replicate.df$replicate <- as.factor(replicate.df$replicate)
  effects.df <- data.frame(factor.name = c("group","sex"), 
                           effect.size = c(dplyr::filter(model.df,id %in% selected.ids)$group.effect.size,dplyr::filter(model.df,id %in% selected.ids)$sex.effect.size),
                           p.value = c(dplyr::filter(model.df,id %in% selected.ids)$group.p.value,dplyr::filter(model.df,id %in% selected.ids)$sex.p.value),
                           stringsAsFactors = F)
  effects.df$factor.name <- factor(effects.df$factor.name, levels = c("group","sex"))
  main.plot <- ggplot(replicate.df,aes(x=replicate,y=value,color=group,fill=sex))+
    geom_violin(trim=F,draw_quantiles=c(0.25,0.5,0.75),alpha=0.25)+facet_grid(~group,scales="free_x")+
    labs(x="Replicate",y="TPM")+theme_minimal()
  effects.plot <- ggplot(effects.df,aes(y=effect.size,x=factor.name,color=factor.name,fill=factor.name))+
    geom_violin(trim=F,draw_quantiles=c(0.25,0.5,0.75),alpha=0.25)+coord_flip()+
    geom_hline(yintercept=0,linetype="longdash",colour="black",size=0.25)+theme_minimal()+
    theme(legend.position="none")+xlab("")+ylab("Effect Size Distribution")

  null.plot <- ggplot(data.frame())+geom_point()+geom_blank()+theme_minimal()
  combined.plot <- gridExtra::arrangeGrob(main.plot,gridExtra::arrangeGrob(null.plot,effects.plot,ncol=1),nrow=1,ncol=2,widths=c(5,2.5))
  return(combined.plot)
}

gseaPlot <- function(selected.set)
{
  plot.df <- model.df %>%
    dplyr::left_join(gene.sets.df %>% dplyr::filter(set.name == selected.set))
  plot.df$set.name[which(is.na(plot.df$set.name))] <- "non.selected"
  plot.df$set.name <- factor(plot.df$set.name, levels = c("non.selected",selected.set))
  factor.names <- c("group","sex")
  gsea.volcano.plot <- lapply(factor.names,function(f)
    plotly::plot_ly(type='scatter',mode="markers",marker=list(size=5),color=plot.df$set.name,colors=c("lightgray","darkred"),x=plot.df[,paste0(f,".effect.size")],y=-log10(plot.df[,paste0(f,".p.value")]),showlegend=F) %>%
      plotly::layout(annotations=list(showarrow=F,x=0.5,y=0.95,align="center",xref="paper",xanchor="center",yref="paper",yanchor="bottom",font=list(size=12,color="darkred"),text=paste0(f," (",dplyr::filter(gsea.df,set.name == selected.set)[,paste0(f,".p.value")],")")),
                     xaxis=list(title=paste0(f," Effect"),zeroline=F),yaxis=list(title="-log10(p-value)",zeroline=F))
  ) %>% plotly::subplot(nrows=1,shareX=F,shareY=T,titleX=T,titleY=T) %>%
    plotly::layout(title=selected.set)
  return(gsea.volcano.plot)
}

Thus:

plot.type.choices <- c('Feature User-Defined Set Plot','Feature Sets GSEA Plot','Feature Plot')

So the first two use ggplot2 for generating each of the two figures they combine, which is then achieved using gridExtra::arrangeGrob. The last one uses plotly.

Here's the shiny code part I've been trying out, but with no luck:

server <- function(input, output)
{
  out.plot <- reactive({
    if(input$plotType == "Feature Plot"){
      out.plot <- featurePlot(selected.id=dplyr::filter(model.df,symbol == input$symbol)$id[1])
    } else if(input$plotType == "Feature User-Defined Set Plot"){
      out.plot <- featureSetPlot(selected.ids=unique(dplyr::filter(model.df,symbol == input$set.symbols)$id))
    } else if(input$plotType == "Feature Sets GSEA Plot"){
      out.plot <- gseaVolcanoPlot(selected.set=input$set.name)
    }
  })

  output$out.plot <- renderPlot({
    if(input$plotType != "Feature Sets GSEA Plot"){
      grid::grid.draw(out.plot())
    } else{
      out.plot()
    }
  })

  output$save <- downloadHandler(
    filename = function() {
      paste0("./plot.pdf")
    },
    content = function(file) {
      ggsave(out.plot(),filename=file,width=10,height=5)
    }
  )
}

ui <- fluidPage(

  tags$style(type="text/css",".shiny-output-error { visibility: hidden; }",".shiny-output-error:before { visibility: hidden; }"),

  titlePanel("Results Explorer"),

  sidebarLayout(

    sidebarPanel(

      # select plot type
      selectInput("plotType","Plot Type",choices=plot.type.choices),

      #in case Feature User-Defined Set Plot was chosen select the genes
      conditionalPanel(condition="input.plotType=='Feature User-Defined Set Plot'",
                       selectizeInput(inputId="set.symbols",label="Features Set Symbols",choices=unique(model.df$symbol),selected=model.df$symbol[1],multiple=T)),

      #in case Feature Sets GSEA Plot was chosen select the databses
      conditionalPanel(condition="input.plotType=='Feature Sets GSEA Plot'",
                       selectizeInput(inputId="set.name",label="Set Name",choices=unique(gene.sets.df$set.name),selected=gene.sets.df$set.name[1],multiple=F)),

      #in case Feature Plot was chosen select the gene
      conditionalPanel(condition="input.plotType=='Feature Plot'",
                       selectizeInput(inputId="symbol",label="Feature Symbol",choices=unique(model.df$symbol),selected=unique(model.df$symbol)[1],multiple=F)),

      downloadButton('save', 'Save to File')
    ),

    mainPanel(
      plotOutput("output.plot")
    )
  )
)

shinyApp(ui = ui, server = server)

I'm suspecting that the renderPlot here may be the issue since I probably have to use plotly::renderPlotly for the Feature Sets GSEA Plot option but I'm not really sure how to tie it all up in the shiny server part.

Another complication that exists and it would be nice to have a solution for is the fact that the gene symbols are not unique WRT gene IDs (as shown in model.df). So it would be nice to have a list that's added if the user selected the Feature Plot option, and that list will show the subset of gene IDs which the selected symbol maps to (dplyr::filter(model.df == input$symbol)$id)

Thanks!

Answer 1

Here's my solution from start to end:

Packages to load:

suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(plotly))
suppressPackageStartupMessages(library(shiny))

Generate example data:

set.seed(1)
model.df <- data.frame(id = paste0("g",1:30),symbol = sample(LETTERS[1:5],30,replace=T),
                       group.effect.size = rnorm(30), group.p.value = runif(30,0,1),
                       sex.effect.size = rnorm(30), sex.p.value = runif(30,0,1),
                       stringsAsFactors = F)
set.seed(1)
design.df <- data.frame(group = c(rep("A",6),rep("B",6)), sex = rep(c(rep("F",3),rep("M",3)),2), replicate = rep(1:6,2)) %>% 
  dplyr::mutate(sample = paste0(group,".",sex,"_",replicate))
design.df$group <- factor(design.df$group, levels = c("A","B"))
design.df$sex <- factor(design.df$sex, levels = c("F","M"))
set.seed(1)
abundance.mat <- matrix(rnorm(30*12), nrow=30, ncol=12, dimnames=list(model.df$id,design.df$sample))
set.seed(1)
gsea.df <- data.frame(set.name = paste0("S",1:4), group.p.value = format(round(runif(4,0,1),2),scientific = T), sex.p.value = format(round(runif(4,0,1),2),scientific = T), stringsAsFactors = F)
set.seed(1)
gene.sets.df <- do.call(rbind,lapply(1:4,function(s) data.frame(set.name = paste0("S",s), id = sample(model.df$id,10,replace = F),stringsAsFactors = F)))
plot.type.choices <- c("Feature Plot","User-Defined Feature Set Plot","Feature Sets GSEA Plot")

Plotting functions:

featurePlot <- function(selected.id)
{
  replicate.df <- reshape2::melt(abundance.mat[which(rownames(abundance.mat) == selected.id),,drop=F], varnames=c("id","sample")) %>%
    dplyr::left_join(design.df)
  effects.df <- data.frame(factor.name = c("group","sex"), 
                           effect.size = c(dplyr::filter(model.df,id == selected.id)$group.effect.size,dplyr::filter(model.df,id == selected.id)$sex.effect.size),
                           p.value = c(dplyr::filter(model.df,id == selected.id)$group.p.value,dplyr::filter(model.df,id == selected.id)$sex.p.value),
                           stringsAsFactors = F)
  effects.df$factor.name <- factor(effects.df$factor.name, levels = c("group","sex"))
  main.plot <- ggplot(replicate.df,aes(x=replicate,y=value,color=group,shape=sex))+
    geom_point(size=3)+facet_grid(~group,scales="free_x")+
    labs(x="Replicate",y="TPM")+theme_minimal()
  xlims <- c(-1*max(abs(effects.df$effect.size))-0.1*max(abs(effects.df$effect.size)),max(abs(effects.df$effect.size))+0.1*max(abs(effects.df$effect.size)))
  effects.plot <- ggplot(effects.df,aes(x=effect.size,y=factor.name,color=factor.name))+
    geom_point()+
    geom_vline(xintercept=0,linetype="longdash",colour="black",size=0.25)+theme_minimal()+xlim(xlims)+
    theme(legend.position="none")+ylab("")+xlab("Effect Size")

  null.plot <- ggplot(data.frame())+geom_point()+geom_blank()+theme_minimal()
  combined.plot <- gridExtra::arrangeGrob(main.plot,gridExtra::arrangeGrob(null.plot,effects.plot,ncol=1),nrow=1,ncol=2,widths=c(5,2.5))
  return(combined.plot)
}


featureSetPlot <- function(selected.ids)
{
  replicate.df <- reshape2::melt(abundance.mat[which(rownames(abundance.mat) %in% selected.ids),,drop=F], varnames=c("id","sample")) %>%
    dplyr::left_join(design.df)
  replicate.df$replicate <- as.factor(replicate.df$replicate)
  effects.df <- data.frame(factor.name = c("group","sex"), 
                           effect.size = c(dplyr::filter(model.df,id %in% selected.ids)$group.effect.size,dplyr::filter(model.df,id %in% selected.ids)$sex.effect.size),
                           p.value = c(dplyr::filter(model.df,id %in% selected.ids)$group.p.value,dplyr::filter(model.df,id %in% selected.ids)$sex.p.value),
                           stringsAsFactors = F)
  effects.df$factor.name <- factor(effects.df$factor.name, levels = c("group","sex"))
  main.plot <- ggplot(replicate.df,aes(x=replicate,y=value,color=group,fill=sex))+
    geom_violin(trim=F,draw_quantiles=c(0.25,0.5,0.75),alpha=0.25)+facet_grid(~group,scales="free_x")+
    labs(x="Replicate",y="TPM")+theme_minimal()
  effects.plot <- ggplot(effects.df,aes(y=effect.size,x=factor.name,color=factor.name,fill=factor.name))+
    geom_violin(trim=F,draw_quantiles=c(0.25,0.5,0.75),alpha=0.25)+coord_flip()+
    geom_hline(yintercept=0,linetype="longdash",colour="black",size=0.25)+theme_minimal()+
    theme(legend.position="none")+xlab("")+ylab("Effect Size Distribution")

  null.plot <- ggplot(data.frame())+geom_point()+geom_blank()+theme_minimal()
  combined.plot <- gridExtra::arrangeGrob(main.plot,gridExtra::arrangeGrob(null.plot,effects.plot,ncol=1),nrow=1,ncol=2,widths=c(5,2.5))
  return(combined.plot)
}

gseaPlot <- function(selected.set)
{
  plot.df <- model.df %>%
    dplyr::left_join(gene.sets.df %>% dplyr::filter(set.name == selected.set))
  plot.df$set.name[which(is.na(plot.df$set.name))] <- "non.selected"
  plot.df$set.name <- factor(plot.df$set.name, levels = c("non.selected",selected.set))
  factor.names <- c("group","sex")
  gsea.plot <- lapply(factor.names,function(f)
    plotly::plot_ly(type='scatter',mode="markers",marker=list(size=5),color=plot.df$set.name,colors=c("lightgray","darkred"),x=plot.df[,paste0(f,".effect.size")],y=-log10(plot.df[,paste0(f,".p.value")]),showlegend=F) %>%
      plotly::layout(annotations=list(showarrow=F,x=0.5,y=0.95,align="center",xref="paper",xanchor="center",yref="paper",yanchor="bottom",font=list(size=12,color="darkred"),text=paste0(f," (",dplyr::filter(gsea.df,set.name == selected.set)[,paste0(f,".p.value")],")")),
                     xaxis=list(title=paste0(f," Effect"),zeroline=F),yaxis=list(title="-log10(p-value)",zeroline=F))
  ) %>% plotly::subplot(nrows=1,shareX=F,shareY=T,titleX=T,titleY=T) %>%
    plotly::layout(title=selected.set)
  return(gsea.plot)
}

Server:

server <- function(input, output)
{
  out.plot <- reactive({
    if(input$plotType == "Feature Plot"){
      out.plot <- featurePlot(selected.id=dplyr::filter(model.df,symbol == input$symbol)$id[1])
    } else if(input$plotType == "User-Defined Feature Set Plot"){
      out.plot <- featureSetPlot(selected.ids=unique(dplyr::filter(model.df,symbol == input$set.symbols)$id))
    } else if(input$plotType == "Feature Sets GSEA Plot"){
      out.plot <- gseaPlot(selected.set=input$set.name)
    }
  })

  output$feature.plot <- renderPlot({
    req(input$plotType == "Feature Plot")
    grid::grid.draw(out.plot())
  })

  output$user.defined.feature.set.plot <- renderPlot({
    req(input$plotType == "User-Defined Feature Set Plot")
    grid::grid.draw(out.plot())
  })

  output$feature.set.gsea.plot <- renderPlotly({
    req(input$plotType == "Feature Sets GSEA Plot")
    out.plot()
  })

  output$save <- downloadHandler(
    filename = function() {
      paste0("./plot.pdf")
    },
    content = function(file) {
      if(input$plotType != "Feature Sets GSEA Plot"){
        ggsave(out.plot(),filename=file,width=10,height=5)
      } else{
        plotly::export(out.plot(),file=file)
      }
    }
  )
}

UI:

ui <- fluidPage(

  tags$style(type="text/css",".shiny-output-error { visibility: hidden; }",".shiny-output-error:before { visibility: hidden; }"),

  titlePanel("Results Explorer"),

  sidebarLayout(

    sidebarPanel(

      # select plot type
      selectInput("plotType","Plot Type",choices=plot.type.choices),

      #in case User-Defined Feature Set Plot was chosen select the genes
      conditionalPanel(condition="input.plotType == 'User-Defined Feature Set Plot'",
                       selectizeInput(inputId="set.symbols",label="Features Set Symbols",choices=unique(model.df$symbol),selected=model.df$symbol[1],multiple=T)),

      #in case Feature Sets GSEA Plot was chosen select the databses
      conditionalPanel(condition="input.plotType == 'Feature Sets GSEA Plot'",
                       selectizeInput(inputId="set.name",label="Set Name",choices=unique(gene.sets.df$set.name),selected=gene.sets.df$set.name[1],multiple=F)),

      #in case Feature Plot was chosen select the gene
      conditionalPanel(condition="input.plotType == 'Feature Plot'",
                       selectizeInput(inputId="symbol",label="Feature Symbol",choices=unique(model.df$symbol),selected=unique(model.df$symbol)[1],multiple=F)),

      downloadButton('save', 'Save to File')
    ),

    mainPanel(
      conditionalPanel(
        condition = "input.plotType == 'User-Defined Feature Set Plot'",
        plotOutput("user.defined.feature.set.plot")
      ),
      conditionalPanel(
        condition =  "input.plotType == 'Feature Sets GSEA Plot'",
        plotly::plotlyOutput("feature.set.gsea.plot")
      ),
      conditionalPanel(
        condition =  "input.plotType == 'Feature Plot'",
        plotOutput("feature.plot")
      )
    )
  )
)

Call:

shinyApp(ui = ui, server = server)

Answer 2

I also guess the problem is "renderPlot". One, not so very elegant way that should solve this problem would be to instead of one output, split it in two, but only ever display one of both using "req()".

This piece of code would become:

 output$out.plot <- renderPlot({
     ....
})

This:

output$out.plot1 <- renderPlot({
     req(input$plotType != "Feature Sets GSEA Plot")
     grid::grid.draw(out.plot())
})
output$out.plot2 <- renderPlotly({
     req(input$plotType == "Feature Sets GSEA Plot")
     out.plot()
})

You can now just add the the plots below each other in you UI. "req()" makes sure absolultely nothing is plotted when the statement inside it is not "truthy" (see ?req), in this case "TRUE". The user would not see a difference between this and replacing one output like you tried.