Shiny Example #4 for Bioconductor

README.md

Bioconductor Shiny Example #4

Example Shiny App showing an example of being embedded in a larger, non-Shiny workflow.

This is an example Shiny app featuring some basic analysis of Ovarian Cancer gene expression data collected on the Affymetrix U133A platform. We filter the available genes and samples down to a much smaller matrix to make reproducibility simpler for a broader audience. The R code involved in sampling the data is available in this Gist as an R-Markdown file, and the sampled data are available in this Gist as Rds files.

To run the application, install shiny (install.packages("shiny")) then run the following command:

library(shiny)
runGist(5944876)

The relevant citation for the original data is available below:

Benjamin Frederick Ganzfried, Markus Riester, Benjamin Haibe-Kains, Thomas Risch, Svitlana Tyekucheva, Ina Jazic, Xin Victoria Wang, Mahnaz Ahmadifar, Michael Birrer, Giovanni Parmigiani, Curtis Huttenhower, Levi Waldron. curatedOvarianData: Clinically Annotated Data for the Ovarian Cancer Transcriptome, Database 2013: bat013 doi:10.1093/database/bat013 published online April 2, 2013.

normalExpr.Rds

ReconstructGRN.Rmd

Compute Gene Network
========================================================

```{r, results='hide', message=FALSE}
library(Biobase)
library(ENA)
library(shiny)
```

## Reconstruct the Full Network

We'll take advantage of the ENA (Ensemble Network Aggregation) package to reconstruct the gene expression network we already loaded.

```{r}
# Load in the sampled matrices we've generated ahead of time.
tumor <- readRDS("tumorExpr.Rds")
# Distill the expression data away from the ExpressionSet object.
tumorExp <- exprs(tumor)

# Reconstruct the gene network. Note that this could potentially be
# distributed across a compute cluster.
gn <- buildWgcna(tumorExp)
```

## Filter the Graph Edges

We'll take a quick look at the distribution of the edge weights in the generated network.

```{r}
weights <- abs(gn)[upper.tri(gn)]
hist(weights)
summary(weights)
```

We'll now filter the edges to only include strong connections.

```{r}
cutoff <- runApp()
filtered <- abs(gn) > cutoff
```

## Write the Output

We'll finally write out the adjacency matrix as a CSV file.

```{r}
write.csv(mat2adj(filtered), file="geneNet.csv")
```

server.R

library(shiny)
library(Biobase)

shinyServer(function(input, output, session) {
  
  filter <- reactive({
    if (!exists("gn"))
      stop(paste("'gn' var doesn't exist. This Shiny App is intended to be run",
                 "as a part of a larger workflow in which some objects would ",
                 "already be defined in this environment prior to `runApp` being", 
                 "executed. Try evaluating the code in `ReconstructGRN.Rmd` which",
                 "wraps this Shiny app in a larger workflow."))
    abs(gn) > input$cutoff
  })
  
  output$tbl <- renderTable({
    filtered <- filter()
    suppressWarnings({head(mat2adj(filtered), n=10)})
  })
  
  output$dim <- renderText({
    filtered <- filter()
    filtered <- filtered[upper.tri(filtered)]
    paste("Edges: ", sum(filtered), sep="")
  })
  
  observe({
    if (input$submit == 0)
      return()
    
    stopApp(input$cutoff)
  })
  
})

setData.Rmd

Extract and Store Some Data
========================================================

We'll first want to install the `curatedOvarianData` package from Bioconductor. Note that this package is a few hundred MB.

```{r, cache=TRUE, results='hide', message=FALSE}
source("http://bioconductor.org/biocLite.R")
biocLite("curatedOvarianData")
```

Now we'll load the package and load in one of the datasets.

```{r, results='hide', message=FALSE}
library(curatedOvarianData)
data(TCGA_eset)
```

Let's write out the clinical data so we have it for later. We'll trim out the "uncurated" metadata column to save some space and make the printing a bit easier.

```{r}
clinical <- phenoData(TCGA_eset)
pData(clinical) <- pData(clinical)[,colnames(pData(clinical))!="uncurated_author_metadata"]
```

We won't want to store all of the tumor samples available in the dataset, so let's grab 50 tumor samples at random and keep all of the normal samples.

```{r}
# Get all normal phenotypic data
normalClinical <- clinical
pData(normalClinical) <- pData(clinical)[pData(clinical)$sample_type == "adjacentnormal",]

# Sample 20 samples from the tumor phenotypic data
tumorIndices <- sample(which(pData(clinical)$sample_type == "tumor"), 50)
tumorClinical <- clinical
pData(tumorClinical) <- pData(clinical)[tumorIndices,]
```

Filtering
---------

For the sake of demonstration, we'll want to keep the dataset fairly small. So let's filter the data to the few dozen samples we selected above and the first 1000 genes.

```{r}
# Extract the relevant probeset to gene mappings
featureList <- fData(TCGA_eset)

# Provide the HGNC symbols of 1000 genes
geneList <- featureList$gene[1:1000]

featureList <- featureList[featureList$gene %in% geneList, ]
featureList <- featureList[match(geneList, rownames(featureList)), ] #order
featureList <- AnnotatedDataFrame(featureList)

# Filter the tumor expression data to only relevant genes and samples
tumor <- exprs(TCGA_eset[,tumorIndices])
tumor <- tumor[match(geneList, rownames(tumor)),]
tumor <- ExpressionSet(tumor, tumorClinical, featureList, experimentData(TCGA_eset))
saveRDS(tumor, "tumorExpr.Rds")
```

And grab all of the normal tissue samples:

```{r}
# Filter the normal expression data to only relevant genes and samples
normal <- exprs(TCGA_eset[,TCGA_eset@phenoData@data$sample_type == "adjacentnormal"])
normal <- normal[match(geneList, rownames(normal)),]
normal <- ExpressionSet(normal, normalClinical, featureList, experimentData(TCGA_eset))
saveRDS(normal, "normalExpr.Rds")
```

Now we should be able to use this sampled data from our apps!

tumorExpr.Rds

ui.R

library(shiny)

shinyUI(pageWithSidebar(
  # Application title
  headerPanel("Hello Shiny Bioconductor Survival!"),
  
  # Sidebar with a selector to choose a gene
  sidebarPanel(
    sliderInput("cutoff", "Cutoff", min=.1, max=1, step=.05, value=.8),
    actionButton("submit", "Use Selected Cutoff")
  ),
  
  # Show a plot of the generated distribution
  mainPanel(
    verbatimTextOutput("dim"),
    tableOutput("tbl")
  )
))

thanks for this example...very useful to create "external configuration" from a shiny app