comparisons for GSEA

comparison_app.R

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                           GLOBAL                          # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  

if (!require(shiny)) install.packages("shiny")
if (!require(shinydashboard)) install.packages("shinydashboard")
require(shiny)
require(shinydashboard)
require(stringr)
require(data.table)
require(dplyr)

# primary samplesheet
primary_ss_csv <- "/Volumes/MOLECPATH/shared/TPU/Runs/NovaSeq/815/run_815_SAMPLESHEET.csv"
primary_ss <- fread(primary_ss_csv)
sample_info_dt <- unique(primary_ss[,.(Sample_ID, Sample_Name)])

# specify files and dir for analysis
analysis_template <- "GSEA_report_RNASeq_template.Rmd"
projdir <- "../"

# specify all the factors to be used in comparisons
factors <- list(
  level1 = c("786-O","A498","U2OS"),
  level2 = c("WT","PBRM1"),
  level3 = c("IFN")
)

# parse sample names (i.e. remove replicate numbers etc.)
# helper functions
replace_NAs_DT <- function(DT, replacement='', cols=colnames(DT)) {
  for (j in cols)
    set(DT,which(is.na(DT[[j]])),j,replacement)
}
paste_cols_DT <- function(DT, sep='', cols=colnames(DT), colname='newcolumn') {
  DT[,I(colname):=Reduce(function(...) paste(..., sep=sep), .SD[,..cols])]
}
# extract factors from sample names
fv <- unlist(unname(factors))
factors_dt <- sample_info_dt[,tstrsplit(Sample_Name,split=" ")][
  ,lapply(.SD,str_extract,pattern=paste(fv,collapse="|"))]
# remove NAs for empty factors
replace_NAs_DT(factors_dt)
# make unique groups
paste_cols_DT(factors_dt, sep=' ', colname='group')
# replace spaces in group labels
sample_info_dt[,sampleConditions:=factors_dt$group][
  ,sampleConditions:=factors_dt$group][,sampleConditions:=str_trim(sampleConditions)][
    ,sampleConditions:=str_replace_all(sampleConditions,"\\s+","_")]
# separate comparison levels
sample_info_dt[,names(factors):=lapply(names(factors), function(x) 
  str_extract(sampleConditions, paste(factors[[x]],collapse="|")))]
# get all the groups that will be compared
comparisons <- function(DT, factors, group_level, compare_level, inter=NULL) 
{
  groups <- unlist(factors[group_level])
  compare <- unlist(factors[compare_level])
  
  if (is.numeric(group_level)) group_level <- names(factors)[group_level]
  if (is.numeric(compare_level)) compare_level <- names(factors)[compare_level]
  
  if (!any((compare_level %in% group_level))) {
    
    DT_list <- split(DT, by=group_level)
    not_used <- !(names(factors) %in% c(group_level, compare_level))
    
    out_list <- lapply(DT_list, function(x) {
      DT <- x[,names(factors)[not_used]:=NULL]
      group <- unique(unlist(DT[,..group_level]))
      if (length(group)>1)
        group <- paste(group,collapse="_")
      treatment <- unique(unlist(DT[,..compare_level]))
      if (any(is.na(treatment)))
        treatment <- treatment[!is.na(treatment)]
      if (length(treatment)==1) {
        control_gp <- group
        treated_gp <-  sprintf("%s_%s",group,treatment)
      } else {
        control_gp <- sprintf("%s_%s",group,treatment[1])
        treated_gp <-  sprintf("%s_%s",group,treatment[2])
      }
      Exp_Name <- sprintf("%s_vs_%s",control_gp,treated_gp)
      Exp_Name
      #list(Exp_Name=Exp_Name, control_gp=control_gp, treated_gp=treated_gp, sample_info=DT)
    })
    
  } else if (any(compare_level %in% group_level)) {
  
    out_list <- lapply(compare_level, function(cl) {
      
      constant_level <- group_level[group_level != cl]
      DT_list <- split(DT, by=constant_level)
      not_used <- !(names(factors) %in% c(group_level, compare_level))
      
      lapply(DT_list, function(x) {
        ns <- x[,..group_level]
        replace_NAs_DT(DT=ns, replacement='')
        names <- unique(paste_cols_DT(ns, sep='_', colname='name')[,name])
        names <- str_replace(str_trim(names), pattern='_$', replacement='')
        
        dm <- matrix(nrow=length(names),ncol=length(names))
        for (i in seq_along(names)) {
          for (j in seq_along(names)) {
            g1 <- names[i]
            g2 <- names[j]
            dm[i,j] <- sprintf("%s_vs_%s", g1, g2)
          }
        }
        dm[upper.tri(dm) & !is.na(dm)]
      })
      
    })
    
  }
  
  names(out_list) <- compare_level
  out_list
  
}

# function to generate reports
render_report <- function(Exp_Name, sample_info, analysis_template) {
  rmarkdown::render(
    analysis_template, 
    params = list(
      Exp_Name = Exp_Name,
      sample_info = sample_info
    ),
    output_file = paste0("GSEA_report_",Exp_Name,".docx")
  )
}


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                           SERVER                          # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  

server <- function(input, output, session) {
  
  # table of level options
  output$levels <- renderTable({
    cbind.data.frame(
      levels = rep(names(factors),sapply(factors, length),each=TRUE), 
      values = c(
        c("786-O","A498","U2OS"),
        c("WT","PBRM1"),
        c("IFN")
      )
    )
  })
  
  # table of comparisons i.e. experiments
  clist <- reactive({
    comparisons(sample_info_dt,factors=factors,group_level=input$group_level,compare_level=input$compare_level)
  })
  
  output$comparisons <- renderTable({
    arrange(cbind.data.frame(Exp_Name = unname(unlist(clist()))), Exp_Name)
  })
  
  # analyses
  eventReactive(
    eventExpr = input$analysis,
    valueExpr = cat("Mock analysis for comparisons starting with ", output$comparisons[1,]),
    ignoreNULL = FALSE
  )
  
  # # report download
  # output$downloa_report <- downloadHandler(
  #   
  #   filename = function() {
  #     paste('report', sep = '.', 'docx') #switch(input$format, PDF='pdf', HTML='html', Word='docx'))
  #   },
  #   
  #   content = function(file) {
  #     src <- normalizePath('report.Rmd')
  #     # temporarily switch to the temp dir, in case you do not have write
  #     # permission to the current working directory
  #     owd <- setwd(tempdir())
  #     on.exit(setwd(owd))
  #     file.copy(src, 'report.Rmd', overwrite=TRUE)
  #     out <- render_report(
  #       Exp_Name = input$Exp_Name, sample_info = input$sample_info, analysis_template = switch(input$analysis, GSEA='')
  #     )
  #     file.rename(out, file)
  #   }
  # )
  
}

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                           UI                              # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
ui <- pageWithSidebar(
  headerPanel(' '),
  sidebarPanel(
    tableOutput('levels'),
    selectInput('group_level', 'grouping levels', names(factors), selected=names(factors)[1], multiple=TRUE, selectize=TRUE),
    selectInput('compare_level', 'comparison levels', names(factors), selected=names(factors)[1], multiple=TRUE, selectize=TRUE),
    submitButton("analysis")
    # downloadButton('download_report')
  ),
  mainPanel(
    tableOutput('comparisons')
  )
)

shinyApp(ui = ui, server = server)