EDA Shiny Rmarkdown script

EDA.RMD

---
title: "Shiny App"
author: "jpaulson"
output: html_document
runtime: shiny
---

This R Markdown document is made interactive using Shiny and allows you to explore an eSet object.

# Examples

Here are a bunch of examples for various shiny apps and the utility of the app.

- [This app!](https://teaching.shinyapps.io/shinyExample/)
- [Gallery](http://shiny.rstudio.com/gallery/)
- [MSD16S](http://epiviz.cbcb.umd.edu/shiny/MSD1000/) ... shameless plug
- [Malnutrition](https://jpaulson.shinyapps.io/helpmeviz-malnutrition/) or [Hunger](https://jpaulson.shinyapps.io/hunger/)

You can 'publish' apps easily online through shinyapps.io or create your own private server, use a gist or github account, etc.

# Onto "our"" data

```{r,include=TRUE,eval=FALSE,message=FALSE,warning=FALSE}
# First we install necessary packages
# source("https://bioconductor.org/biocLite.R")
library(BiocInstaller)
pkgs = c("Biobase","shiny","RColorBrewer","metagenomeSeq","vegan","rafalib")
biocLite(pkgs,suppressUpdates = TRUE,verbose=FALSE,quiet=TRUE)
```

Replace `data = eSet` with your own eSet object's name.

```{r, echo=TRUE,warning=FALSE}
suppressPackageStartupMessages(library(metagenomeSeq))
suppressPackageStartupMessages(library(RColorBrewer))
suppressPackageStartupMessages(library(parathyroidSE))
data("parathyroidGenesSE")
counts = assay(parathyroidGenesSE)
phenodata = AnnotatedDataFrame(as.data.frame(colData(parathyroidGenesSE)))
paraESet = ExpressionSet(counts,phenoData=phenodata)
data = paraESet
```

# Metadata

## Feature information
```{r,echo=FALSE}
renderDataTable({
          fd = fData(data)
          Index = 1:nrow(fd)
          Rownames = rownames(fd)
          fd = cbind(Index,Rownames,fd)
          return(fd)
        },options = list(iDisplayLength = 10))
```

## Phenotype information
```{r,echo=FALSE}
renderDataTable({
          return(pData(data))
        },options = list(iDisplayLength = 10))
```

# Plots

## Input for all plots
```{r,echo=FALSE}
inputPanel(
  checkboxInput("norm","Log2-transform",value=TRUE),
  numericInput("pd","Phenotype column:",3,min=1)
)
```

### Feature Abundance
```{r,echo=FALSE}
inputPanel(
  numericInput("featureIndex","Display feature (index):",1,min=1)
  )
renderPlot({
  mat = exprs(data)
  if(input$norm){
    mat = log2(exprs(data)+1)
  }
  pd = factor(pData(data)[,input$pd])
  ylabel = ifelse(input$norm,yes=expression(bold("Log"[2]*" Abundance")),no="No. raw reads")
  metagenomeSeq::plotFeature(mat,otuIndex = input$featureIndex,ylab=ylabel,main=rownames(mat)[input$featureIndex],
  classIndex = list(All_samples=1:ncol(mat)),col=pd,font.lab=2,font.axis=1,sort=FALSE,xaxt="n")
  legend("topleft",legend=unique(pd),fill=unique(pd),box.col="NA")
})
```

### Heatmap
```{r,echo=FALSE}
inputPanel(
  numericInput("noFeatures","Number of features:",15,min=1,max =200),
  radioButtons("heatMethod","Choose features by:",c("Median Absolute Deviation"="mad","Variability"="sd"))
)
renderPlot({
  mat = exprs(data)
  if(input$norm){
    mat = log2(exprs(data)+1)
  }
  trials = pData(data)[,input$pd]
  otusToKeep = which(rowSums(mat) > 0)
  otuStats = apply(mat[otusToKeep, ], 1, input$heatMethod)
  otuIndices = otusToKeep[order(otuStats, decreasing = TRUE)[1:input$noFeatures]]
  my_mat <- mat[otuIndices,]
  heatmapColColors=brewer.pal(12,"Set3")[as.integer(factor(trials))];
  heatmapCols = colorRampPalette(brewer.pal(9, "RdBu"))(50)
  gplots::heatmap.2(my_mat,trace="none",cexRow=.8,cexCol=.8,col = heatmapCols,ColSideColors = heatmapColColors)
  legend("left",fill=unique(heatmapColColors),legend=unique(trials))
},height=800)
```


### PCA/PCoA
```{r,echo=FALSE}
inputPanel(
      radioButtons("pcaOrMds","PCA or MDS:",c("PCA"="TRUE","MDS"="FALSE"),selected="TRUE"),
      radioButtons("useDist","Make use of count distances:",c("False"="FALSE","True"="TRUE"),selected="FALSE"),
      conditionalPanel(condition = "input.useDist == 'TRUE'",
        selectInput("distance", "Distance:", 
          choices=c("euclidean","manhattan","canberra","bray",
            "kulczynski","jaccard","gower","altGower","morisita",
            "horn","raup","binomial","chao","cao"),selected="euclidean")
      ),
      numericInput('dimensionx', 'X-axis dimension:', 1,
         min = 1, max = 4),
      numericInput('dimensiony', 'Y-axis dimension:', 2,
         min = 1, max = 4)
)
renderPlot({
  mat = exprs(data)
  if(input$norm){
    mat = log2(exprs(data)+1)
  }
  useDist = input$useDist
  pd  = factor(pData(data)[,input$pd])
  metagenomeSeq::plotOrd(mat,n=100,pch=21,bg=pd,usePCA=input$pcaOrMds,
  comp=c(input$dimensionx,input$dimensiony),
  useDist=useDist,distfun=vegan::vegdist,dist.method=input$distance)
  legend("topleft",levels(pd),fill=factor(levels(pd)),box.col="NA")
})
```

### Diversity
```{r,echo=FALSE}
inputPanel(selectInput("diversity","Diversity index:",choices=c("shannon", "simpson", "invsimpson")))
renderPlot({
  mat = exprs(data)
  if(input$norm){
    mat = log2(exprs(data)+1)
  }
  mat = t(mat)
  pd  = factor(pData(data)[,input$pd])
  H = vegan::diversity(mat,index=input$diversity)
  boxplot(H~pd,ylab=paste(input$diversity,"diversity index"))
})
```