inst/old/ggshine10.R

In vjcitn/gQTLbrowser: interactive exploration of gQTL results

gQTLbrowse3 = function( store, baseSE, 
   stateGR, phenGR, FDRsupp, orgDbObj=Homo.sapiens, selector=selectizeInput ) {
#
# interface to shiny/ggvis eqtl exploration
# assumes central identifier is the gene symbol
#
   #
   # get all available gene symbols
   #
   
   allsyms = keys(orgDbObj, keytype="SYMBOL")
   
#
#  filter the symbols relevant to the input SummarizedExperiment baseSE,
#  which will generally not use gene symbols as rownames ... in fact
#  FIXME -- we are assuming existence of gene_name and gene_type in
#  rowData, as for geuvPack geuFPKM
#
#  this sort of symbol mapping exercise is common and should be
#  substantially abstracted -- perhaps a mapping between GSEABase entities
#
   
   availProbes = store@probemap$probeid
   availSyms = rowData(baseSE[ availProbes, ])$gene_name
   sorts = sort(availSyms)
   symok = which(availSyms %in% allsyms)
   availTypes = rowData(baseSE[ availProbes, ])$gene_type
   p2g = availSyms = availSyms[symok]
   p2t = availTypes[symok]
   availProbes = availProbes[symok]
   names(p2g) = availProbes
   names(p2t) = availProbes

input = list()
input$sym = availSyms[2]
   
   #
   # we want to get the symbols from shiny, but use
   # ggvis for tooltips
   #
   
#   server = function(input, output) {
   
   ## first, extract the ENSEMBL ID for selected symbol
   
                         paste0("GEUVADIS ENSEMBL ID: ",
                            availProbes[ which(availSyms == input$sym)[1] ] )
   
   ## second, acquire the appropriate eQTL testing results
   ##   and bind information on chromatin state, transcript location
   
   #
   # get the GRanges with eQTL results
   #
        n1 = extractByProbes( store, 
                  availProbes[ which(availSyms == input$sym)[1] ] )
   #
   # obtain chromatin state labels
   #
        n1$st878 = rep("none", length(n1))
        fo = findOverlaps(n1, stateGR)
        n1$st878[ queryHits(fo) ] = as.character(stateGR$name)[ subjectHits(fo) ]
   #     uniqst = unique(stateGR$name) # useless effort at persistent colormap
   #     nuniqst = length(uniqst)
   #     cmap = colorRampPalette(c("red", "blue"))(nuniqst) # vector of codes
   #     n1$col878 = cmap[as.numeric(factor(n1$st878))]
   #
   # execute the FDR filter
   #
        n1 = FDRsupp@filterUsed(n1)
   #
   # compute FDR
   #
        n1$ml10FDR = pmin(6, -log10(getFDRfunc(FDRsupp)(n1$chisq)))
   #     n1 <<- n1
   #
   # build data frame for visualization
   #
        mydf <- data.frame(chr=as.character(seqnames(n1)), pos=start(n1),
     	  MAF = n1$MAF,
               probeid=n1$probeid, snp=n1$snp, ml10FDR = n1$ml10FDR,
               stringsAsFactors=FALSE, Mb=start(n1)/1e6, st878=n1$st878)
   #
   # use global maps to recover symbol and GEUVADIS 'gene type'
   #
        mydf$gene = as.character( p2g[ mydf$probeid ] )
        mydf$type = as.character( p2t[ mydf$probeid ] )
   #
   # obtain a 'gene model' for the selected symbol, so that
   # locations of transcripts can be given
   #
        mod = gmod2( input$sym )
   #
   # add the location information, "faking" fields for eQTL results
   #
        extra = tail(mydf, length(mod))
        extra$st878 = paste0("TXLOC(", input$sym,")")
        extra$Mb = start(mod)/1e6
        extra$ml10FDR = -.25
        extra$snp = extra$MAF = NA
        mydf = rbind(mydf, extra)
   #
   # get disease loci
   #
        disdat = phenGR[ which(phenGR$geuvvid %in% mydf$snp) ]
        if (length(disdat) > 0) {
           extra2 = tail(mydf, length(disdat))
           extra2$ml10FDR = -.4
           extra2$MAF = NA
           extra2$Mb = start(disdat)/1e6
           extra2$st878 = paste0("  trait: ", disdat$Disease.Trait)
           mydf = rbind(mydf, extra2)
           }
   #     mydf = mydf[ order(mydf$st878), ]
   #
   # construct the key for tooltip
   #
        mydf$rowid = 1:nrow(mydf)
        mydf <<- mydf
        vals = mydf %>% filter( gene == input$sym )  # do earlier !?!
   
  #    P1 = reactive( {
  #      validate( 
  #        need( input$sym != "", FALSE )
  #      )
  #       all_values <- function(x) {
  #           if(is.null(x)) return(NULL)
  #           row <- mydf[mydf$rowid == x$rowid, ]
  #           paste0(names(row), ": ", format(row), collapse = "<br />")
  #         }
  #       filteredData %>% ggvis(~Mb, ~ml10FDR, key := ~rowid,
  #                fill = ~st878) %>% 
  #             layer_points() %>%
  #             add_tooltip(all_values, "hover") %>% layer_points() %>%
  ##             add_legend("fill", title=paste0(deparse(substitute(stateGR)), " state"), values=unique(mydf$st878)) %>% 
  # 
  #             add_axis("y", title=paste0("-log10 FDR assoc w/ ", input$sym, " expr" ))
  #       } )
  #    P1 %>% bind_shiny("p")
  # }
  # 
  # shinyApp(ui = ui, server=server)
}