FET demo

README.md

Shiny app for exploring Fisher's exact tests.

See a running demo at https://daler.shinyapps.io/lcdb_fet/

Run locally with the following command in this directory:

Rscript -e "shiny::runApp()"

Or follow directions on shinyapps.io to upload to your account there.

include.md

Input controls

Input controls can be used to build a table. Tables, pval, odds ratio, and figures are updated automatically.

Tables

Tables are shown both as total counts as well as in versions that show row percentages and column percentages to help with interpretation.

Mosaic plot

The mosaic plot shows the relative size of cells in table. When a table is not different from the null hypothesis, the gaps in the mosaic plot will form a cross-shape. The more divergent from a cross-shape, the higher the odds ratio.

Line plot

The black line shows the expected distribution of the top-left cell of the table ("upregulated AND peak"), given the marginal values of the table (total genes assayed, total genes with a peak, total upregulated genes). Try modifying those values and watch the position of the null distribution shift.

The blue vertical line shows the current value of "upregulated AND peak". When the blue line is close to the middle of the black distribution, it is close to the expected value and therefore we get a non-significant p-value and odds ratios close to 1.

server.R

library(shiny)
library(vcd)
library(Cairo)
options(shiny.usecairo=TRUE)

get.table <- function(input){
    r1c1 <- input$interesting
    r1c2 <- input$genes_with_peak - r1c1
    r2c1 <- input$upregulated - r1c1
    r2c2 <- input$genes - r1c1 - r1c2 - r2c1
    m <- matrix(
        c(
          r1c1, r1c2,
          r2c1, r2c2),
        nrow=2,
        dimnames=list(
            peak=c('peak in promoter', 'no peak'),
            genes=c('upregulated', 'not upregulated')
            )
        )
}


fet <- function(X){
    m <- sum(X[, 1L])  # first column sum
    n <- sum(X[, 2L])  # second column sum
    k <- sum(X[1L, ])  # first row sum
    x <- X[1L, 1L]     # value of upper left cell
    lo <- max(0L, k - n)
    hi <- min(k, m)
    support <- lo:hi
    logdc <- dhyper(support, m, n, k, log = TRUE)
    dnhyper <- function(ncp) {
        d <- logdc + log(ncp) * support
        d <- exp(d - max(d))
        d/sum(d)
    }
    oddsratio <- (X[1,1]/X[1,2])/(X[2,1]/X[2,2])
    d <- dnhyper(1)
    PVAL <- sum(d[d <= d[x - lo + 1]])
    return(list(pval=PVAL, d=d, oddsratio=oddsratio, x=x))
}

formatDecimal <- function(x, k) format(round(x, k), trim=T, nsmall=k)


# Define server logic required to draw a histogram
shinyServer(
    function(input, output) {
        output$distPlot <- renderPlot({
            X <- get.table(input)
            res <- fet(X)
            LIM <- input$xmax
            d <- res[['d']]
            x <- res[['x']]
            plot(d[1:LIM], type='l', lwd=2, xlab='Number of upregulated genes with peak', ylab='density')
            lines(d[d<=d[x+1]], col='red', lwd=2)
            abline(v=x, col='blue', lwd=2)
        })
        output$mosaicplot <- renderPlot({
            X <- get.table(input)
            mosaic(X, main="")
        })
        output$counts <- renderTable({
            X <- get.table(input)
            y <- as.data.frame(addmargins(X))
            y[, c(1, 2,3)] <- apply(y[, c(1,2,3)], 1, as.integer)
            y
        }, caption='2x2 table of counts, with marginal sums',
           caption.placement='bottom',
           caption.width=NULL
        )

        output$rowperc <- renderTable({
            X <- addmargins(get.table(input))
            X <- X / X[,3] * 100
            X[1:2,]
        }, caption='Fraction of genes that are upregulated or not',
           caption.placement='bottom',
           caption.width=NULL
        )
        output$colperc <- renderTable({
            X <- addmargins(t(get.table(input)))
            X <- X / X[,3] * 100
            X[1:2,]
            t(X)[,1:2]
        }, caption='Fraction of genes that have a peak or not',
           caption.placement='bottom',
           caption.width=NULL
        )
        output$pval <- renderText({
            X <- get.table(input)
            res <- fet(X)
            paste('pval:', signif(res[['pval']], 4))
        })
        output$or <- renderText({
            X <- get.table(input)
            res <- fet(X)
            paste('odds ratio:', signif(res[['oddsratio']], 4))
        })
        output$results <- renderText({
            X <- get.table(input)
            res <- fet(X)
            or <- res[['oddsratio']]
            pval <- res[['pval']]
            if (or < 1) {direction <- 'depleted'} else {direction <- 'enriched'}
            if (pval > 0.05){ direction <- 'n.s'}
            direction
        })
})

ui.R

library(shiny)

# Define UI for application that draws a histogram
shinyUI(
  fluidPage(
      # Application title
      fluidRow(
        column(
          2,
          wellPanel(
            h3('LCDB data club FET demo'),
            numericInput("genes",
                        "total genes assayed",
                        min = 1,
                        max = 25000,
                        value = 15000,
                        step=100
                        ),
            numericInput("genes_with_peak",
                        "total genes with a peak",
                        min = 1,
                        max = 10000,
                        value = 2130,
                        step=10),
            numericInput("upregulated",
                        "total upregulated genes",
                        min = 1,
                        max = 10000,
                        value = 930,
                        step=10
                        ),
            numericInput("interesting",
                        "upregulated AND peak",
                        min = 1,
                        max = 5000,
                        value = 130,
                        step=5),
               sliderInput("xmax",
                           "xmax",
                           min=1,
                           max=5000,
                           value=300)
            )
          ),
          column(
            5,
            wellPanel(
                      h3(textOutput("pval")),
                      h3(textOutput("or")),
                      h3(textOutput("results"))
            ),
            wellPanel(
                h4(tableOutput("counts"))
            ),
                      tableOutput("rowperc"),
                      tableOutput("colperc")
           ),

          column(
            5,
                      plotOutput("mosaicplot"),
                      plotOutput("distPlot")
          )
      )
  )
)