inst/doc/biobroom_vignette.R

In biobroom: Turn Bioconductor objects into tidy data frames

## ----global_options, include=FALSE, echo=FALSE--------------------------------
knitr::opts_chunk$set(fig.width=12, fig.height=6, out.width='700in', out.height='350in', 
                      echo=TRUE, warning=FALSE, message=FALSE, cache=FALSE, dev='png')

## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("biobroom")

## ----eval=FALSE---------------------------------------------------------------
#  library(biobroom)
#  ?edgeR_tidiers
#  ?DESeq2_tidiers
#  ?limma_tidiers
#  ?ExpressionSet_tidiers
#  ?MSnSet_tidiers
#  ?qvalue_tidiers

## -----------------------------------------------------------------------------
library(qvalue)
data(hedenfalk)

qobj <- qvalue(hedenfalk$p)
names(qobj)

## -----------------------------------------------------------------------------
library(biobroom)
# use tidy/augment/glance to restructure the qvalue object
head(tidy(qobj))
head(augment(qobj))
head(glance(qobj))

## -----------------------------------------------------------------------------
# create sample names
df <- data.frame(gene = 1:length(hedenfalk$p))
head(augment(qobj, data = df))

## -----------------------------------------------------------------------------
library(ggplot2)
# use augmented data to compare p-values to q-values
ggplot(augment(qobj), aes(p.value, q.value)) + geom_point() +
  ggtitle("Simulated P-values versus Computed Q-values") + theme_bw()

## -----------------------------------------------------------------------------
library(dplyr)

# Find significant genes under 0.05 threshold
sig.genes <- augment(qobj) %>% filter(q.value < 0.05)
head(sig.genes)

## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("airway")

## -----------------------------------------------------------------------------
library(DESeq2)
library(airway)

data(airway)
airway_se <- airway

## -----------------------------------------------------------------------------
airway_dds <- DESeqDataSet(airway_se, design = ~cell + dex)

head(tidy(airway_dds))

## -----------------------------------------------------------------------------
# differential expression analysis
deseq <- DESeq(airway_dds)
results <- results(deseq)
# tidy results
tidy_results <- tidy(results)
head(tidy_results)

## -----------------------------------------------------------------------------
ggplot(tidy_results, aes(x=estimate, y=log(p.value),
                         color=log(baseMean))) + geom_point(alpha=0.5) +
  ggtitle("Volcano Plot For Airway Data via DESeq2") + theme_bw()

## -----------------------------------------------------------------------------
library(edgeR)
data(hammer)

hammer.counts <- Biobase::exprs(hammer)[, 1:4]
hammer.treatment <- Biobase::phenoData(hammer)$protocol[1:4]

y <- DGEList(counts=hammer.counts,group=hammer.treatment)
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)
et <- exactTest(y)

## -----------------------------------------------------------------------------
head(tidy(et))

## -----------------------------------------------------------------------------
glance(et, alpha = 0.05)

## -----------------------------------------------------------------------------
ggplot(tidy(et), aes(x=estimate, y=log(p.value), color=logCPM)) +
  geom_point(alpha=0.5) + ggtitle("Volcano Plot for Hammer Data via EdgeR") +
  theme_bw()

## -----------------------------------------------------------------------------
# create random data and design
dat <- matrix(rnorm(1000), ncol=4)
dat[, 1:2] <- dat[, 1:2] + .5  # add an effect
rownames(dat) <- paste0("g", 1:nrow(dat))
des <- data.frame(treatment = c("a", "a", "b", "b"),
                  confounding = rnorm(4))

## -----------------------------------------------------------------------------
lfit <- lmFit(dat, model.matrix(~ treatment + confounding, des))
eb <- eBayes(lfit)

head(tidy(lfit))
head(tidy(eb))

## -----------------------------------------------------------------------------
ggplot(tidy(eb), aes(x=estimate, y=log(p.value), color=statistic)) + 
  geom_point() + ggtitle("Nested Volcano Plots for Simulated Data Processed with limma") +
  theme_bw()
  

## -----------------------------------------------------------------------------
library(Biobase)

head(tidy(hammer))

## -----------------------------------------------------------------------------
head(tidy(hammer, addPheno = TRUE))

## -----------------------------------------------------------------------------
ggplot(tidy(hammer, addPheno=TRUE), aes(x=protocol, y=log(value))) +
  geom_boxplot() + ggtitle("Boxplot Showing Effect of Protocol on Expression")

## -----------------------------------------------------------------------------
library(MSnbase)
data(msnset)

head(tidy(msnset))

head(tidy(msnset, addPheno = TRUE))

## ----eval=FALSE---------------------------------------------------------------
#  options(biobroom.return = "data.frame")
#  options(biobroom.return = "data.table")