# biobroom Vignette

#### Andrew J. Bass, Emily Nelson, David Robinson and John D. Storey

#### 2016-10-17

# About biobroom

The `biobroom` package contains methods for converting standard objects in Bioconductor into a “tidy format”. It serves as a complement to the popular broom package, and follows the same division (`tidy`/`augment`/`glance`) of tidying methods.

“Tidy data” is a data analysis paradigm that focuses on keeping data formatted as a single observation per row of a data table. For further information, please see Hadley Wickham’s seminal paper on the subject. “Tidy” is not a normative statement about the quality of an object’s structure. Rather, it is a technical specification about the choice of rows and columns. A tidied data frame is not “better” than an S4 object; it simply allows analysis with a different set of tools.

Tidying data makes it easy to recombine, reshape and visualize bioinformatics analyses. Objects that can be tidied include

- `ExpressionSet` object,
- `GRanges` and `GRangesList` objects,
- `RangedSummarizedExperiment` object,
- `MSnSet` object,
- per-gene differential expression tests from `limma`, `edgeR`, and `DESeq2`,
- `qvalue` object for multiple hypothesis testing.

We are currently working on adding more methods to existing Bioconductor objects. If any bugs are found please contact the authors or visit our github page. Otherwise, any questions can be answered on the Bioconductor support site.

# Installation

The `biobroom` package can be installed by typing in a R terminal:

    source("http://bioconductor.org/biocLite.R")
    biocLite("biobroom")

To find out more about the provided objects:

    library(biobroom)
    ?edgeR_tidiers
    ?DESeq2_tidiers
    ?limma_tidiers
    ?ExpressionSet_tidiers
    ?MSnSet_tidiers
    ?qvalue_tidiers

# Examples

## qvalue object

The qvalue package is a popular package to estimate q-values and local false discovery rates. To get started, we can load the `hedenfalk` dataset included in the `qvalue` package:

    library(qvalue)
    data(hedenfalk)

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

    ## [1] "call"       "pi0"        "qvalues"    "pvalues"    "lfdr"
    ## [6] "pi0.lambda" "lambda"     "pi0.smooth"

`qobj` is a `qvalue` object, generated from the p-values contained in the `hedenfalk` dataset. If we wanted to use a package such as `dplyr` or `ggplot`, we would need to convert the results into a data frame. The `biobroom` package makes this conversion easy by using the `tidy`, `augment` and `glance` functions:

- `tidy` returns one row for each choice of the tuning parameter lambda.
- `augment` returns one row for each provided p-value, including the computed q-value and local false discovery rate.
- `glance` returns a single row containing the estimated `pi0`.

Applying these functions to `qobj`:

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

    ## # A tibble: 6 × 3
    ##   lambda smoothed       pi0
    ##    <dbl>    <lgl>     <dbl>
    ## 1   0.05    FALSE 0.8517350
    ## 2   0.10    FALSE 0.8068700
    ## 3   0.15    FALSE 0.7823344
    ## 4   0.20    FALSE 0.7563091
    ## 5   0.25    FALSE 0.7310200
    ## 6   0.30    FALSE 0.7138351

    head(augment(qobj))

    ## # A tibble: 6 × 3
    ##      p.value    q.value      lfdr
    ##        <dbl>      <dbl>     <dbl>
    ## 1 0.01212618 0.08819163 0.1680943
    ## 2 0.07502524 0.20936729 0.4128630
    ## 3 0.99492114 0.66799864 1.0000000
    ## 4 0.04178549 0.16163643 0.3088248
    ## 5 0.84581388 0.63247386 1.0000000
    ## 6 0.25192429 0.37172329 0.6924674

    head(glance(qobj))

    ## # A tibble: 1 × 2
    ##        pi0 lambda
    ##      <dbl>  <dbl>
    ## 1 0.669926   0.95

The original data, or in this example the gene names, can be inputted into `augment` using the `data` argument:

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

    ## # A tibble: 6 × 4
    ##    gene    p.value    q.value      lfdr
    ##   <int>      <dbl>      <dbl>     <dbl>
    ## 1     1 0.01212618 0.08819163 0.1680943
    ## 2     2 0.07502524 0.20936729 0.4128630
    ## 3     3 0.99492114 0.66799864 1.0000000
    ## 4     4 0.04178549 0.16163643 0.3088248
    ## 5     5 0.84581388 0.63247386 1.0000000
    ## 6     6 0.25192429 0.37172329 0.6924674

The tidied data can be used to easily create plots:

    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()

Additionally, we can extract out important information such as significant genes under a false discovery rate threshold:

    library(dplyr)

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

    ## # A tibble: 6 × 3
    ##        p.value    q.value       lfdr
    ##          <dbl>      <dbl>      <dbl>
    ## 1 0.0007129338 0.02315186 0.03652658
    ## 2 0.0017507886 0.03645186 0.05965671
    ## 3 0.0026498423 0.04292489 0.07511774
    ## 4 0.0011640379 0.03014667 0.04757577
    ## 5 0.0029842271 0.04479572 0.08021691
    ## 6 0.0023533123 0.03966387 0.07033209

## DESeq2 objects

To demonstrate tidying on `DESeq2` objects we have used the published `airway` RNA-Seq experiment, available as a package from Bioconductor:

    source("https://bioconductor.org/biocLite.R")
    biocLite("airway")

Import the `airway` dataset:

    library(DESeq2)
    library(airway)

    data(airway)
    airway_se <- airway

`airway_se` is a `SummarizedExperiment` object, which is a type of object used by the `DESeq2` package. Next, we create a `DESeqDataSet` object and show the output of tidying this object:

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

    head(tidy(airway_dds))

    ## # A tibble: 6 × 3
    ##              gene     sample count
    ##             <chr>      <chr> <int>
    ## 1 ENSG00000000003 SRR1039508   679
    ## 2 ENSG00000000005 SRR1039508     0
    ## 3 ENSG00000000419 SRR1039508   467
    ## 4 ENSG00000000457 SRR1039508   260
    ## 5 ENSG00000000460 SRR1039508    60
    ## 6 ENSG00000000938 SRR1039508     0

Only the gene counts are outputted since there has been no analysis performed. We perform an analysis on the data and then `tidy` the resulting object:

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

    ## # A tibble: 6 × 7
    ##              gene    baseMean    estimate   stderror  statistic
    ##             <chr>       <dbl>       <dbl>      <dbl>      <dbl>
    ## 1 ENSG00000000003 708.6021697  0.37415246 0.09884435  3.7852692
    ## 2 ENSG00000000005   0.0000000          NA         NA         NA
    ## 3 ENSG00000000419 520.2979006 -0.20206175 0.10974241 -1.8412367
    ## 4 ENSG00000000457 237.1630368 -0.03616686 0.13834540 -0.2614244
    ## 5 ENSG00000000460  57.9326331  0.08445399 0.24990709  0.3379415
    ## 6 ENSG00000000938   0.3180984  0.08413901 0.15133424  0.5559813
    ## # ... with 2 more variables: p.value <dbl>, p.adjusted <dbl>

As an example to show how easy it is to manipulate the resulting object, `tidy_results`, we can use `ggplot2` to create a volcano plot of the p-values:

    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()

## edgeR objects

Here we use the `hammer` dataset included in `biobroom` package. `edgeR` can be used to perform differential expression analysis as follows:

    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)

The results of the analysis are stored in `et`, which is an `DGEExact` object. We can `tidy` this object using `biobroom`:

    head(tidy(et))

    ## # A tibble: 6 × 4
    ##                 gene    estimate      logCPM      p.value
    ##                <chr>       <dbl>       <dbl>        <dbl>
    ## 1 ENSRNOG00000000001  2.64635814  1.49216267 1.309933e-06
    ## 2 ENSRNOG00000000007 -0.40869816 -0.22616605 1.000000e+00
    ## 3 ENSRNOG00000000008  2.22296029 -0.40665547 1.288756e-01
    ## 4 ENSRNOG00000000009  0.00000000 -1.31347471 1.000000e+00
    ## 5 ENSRNOG00000000010  0.03307909  1.79448965 1.000000e+00
    ## 6 ENSRNOG00000000012 -3.39210151  0.07939132 3.745676e-03

`glance` shows a summary of the experiment: the number of genes found significant (at a specified `alpha`), and which contrasts were compared to get the results.

    glance(et, alpha = 0.05)

    ##   significant     comparison
    ## 1        6341 control/L5 SNL

Additionally, we can can easily manipulate the resulting object and create a volcano plot of the p-values using `ggplot2`:

    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()

## limma objects

To demonstrate how `biobroom` works with `limma` objects, we generate some simulated data to test the tidier for `limma` objects.

    # 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))

We then use `lmFit` and `eBayes` (functions included in `limma`) to fit a linear model and use `tidy` to convert the resulting object into tidy format:

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

    head(tidy(lfit))

    ## # A tibble: 6 × 3
    ##    gene       term   estimate
    ##   <chr>      <chr>      <dbl>
    ## 1    g1 treatmentb -3.8761430
    ## 2    g2 treatmentb -3.7243752
    ## 3    g3 treatmentb  1.3402339
    ## 4    g4 treatmentb -4.7675767
    ## 5    g5 treatmentb  0.2044367
    ## 6    g6 treatmentb -4.3279820

    head(tidy(eb))

    ## # A tibble: 6 × 6
    ##    gene       term   estimate  statistic    p.value       lod
    ##   <chr>      <chr>      <dbl>      <dbl>      <dbl>     <dbl>
    ## 1    g1 treatmentb -3.8761430 -2.4656893 0.07989961 -4.591616
    ## 2    g2 treatmentb -3.7243752 -2.5342733 0.07478014 -4.591510
    ## 3    g3 treatmentb  1.3402339  0.9860901 0.38853994 -4.595164
    ## 4    g4 treatmentb -4.7675767 -1.5222121 0.21421312 -4.593618
    ## 5    g5 treatmentb  0.2044367  0.1298449 0.90395019 -4.596995
    ## 6    g6 treatmentb -4.3279820 -2.9632357 0.05033173 -4.590939

Analysis can easily be performed from the tidied data. The package `ggplot2` can be used to make a volcano plot of the p-values:

    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()

## ExpressionSet objects

`tidy` can also be run directly on `ExpressionSet` objects, as described in another popular `Bioconductor` package `Biobase.` The `hammer` dataset we used above (which is included in the `biobroom` package) is an `ExpressionSet` object, so we’ll use that to demonstrate.

    library(Biobase)

    head(tidy(hammer))

    ## # A tibble: 6 × 3
    ##                 gene    sample value
    ##                <chr>     <chr> <int>
    ## 1 ENSRNOG00000000001 SRX020102     2
    ## 2 ENSRNOG00000000007 SRX020102     4
    ## 3 ENSRNOG00000000008 SRX020102     0
    ## 4 ENSRNOG00000000009 SRX020102     0
    ## 5 ENSRNOG00000000010 SRX020102    19
    ## 6 ENSRNOG00000000012 SRX020102     7

We can add the phenotype information by using the argument `addPheno`:

    head(tidy(hammer, addPheno = TRUE))

    ## # A tibble: 6 × 8
    ##                 gene    sample sample.id num.tech.reps protocol
    ##                <chr>     <chr>    <fctr>         <dbl>   <fctr>
    ## 1 ENSRNOG00000000001 SRX020102 SRX020102             1  control
    ## 2 ENSRNOG00000000007 SRX020102 SRX020102             1  control
    ## 3 ENSRNOG00000000008 SRX020102 SRX020102             1  control
    ## 4 ENSRNOG00000000009 SRX020102 SRX020102             1  control
    ## 5 ENSRNOG00000000010 SRX020102 SRX020102             1  control
    ## 6 ENSRNOG00000000012 SRX020102 SRX020102             1  control
    ## # ... with 3 more variables: strain <fctr>, Time <fctr>, value <int>

Now we can easily visualize the distribution of counts for each protocol by using `ggplot2`:

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

## MSnSet Objects

`tidy` can also be run directly on `MSnSet` objects from `MSnbase`, which as specialised containers for quantitative proteomics data.

    library(MSnbase)
    data(msnset)

    head(tidy(msnset))

    ## # A tibble: 6 × 3
    ##   protein  sample.id      value
    ##     <chr>      <chr>      <dbl>
    ## 1      X1 iTRAQ4.114  1347.6158
    ## 2     X10 iTRAQ4.114   739.9861
    ## 3     X11 iTRAQ4.114 27638.3582
    ## 4     X12 iTRAQ4.114 31892.8928
    ## 5     X13 iTRAQ4.114 26143.7542
    ## 6     X14 iTRAQ4.114  6448.0829

    head(tidy(msnset, addPheno = TRUE))

    ## # A tibble: 6 × 5
    ##   protein     sample    mz reporters      value
    ##    <fctr>     <fctr> <dbl>    <fctr>      <dbl>
    ## 1      X1 iTRAQ4.114 114.1    iTRAQ4  1347.6158
    ## 2     X10 iTRAQ4.114 114.1    iTRAQ4   739.9861
    ## 3     X11 iTRAQ4.114 114.1    iTRAQ4 27638.3582
    ## 4     X12 iTRAQ4.114 114.1    iTRAQ4 31892.8928
    ## 5     X13 iTRAQ4.114 114.1    iTRAQ4 26143.7542
    ## 6     X14 iTRAQ4.114 114.1    iTRAQ4  6448.0829

# Note on returned values

All biobroom `tidy` and `augment` methods return a `tbl_df` by default (this prevents them from printing many rows at once, while still acting like a traditional `data.frame`). To change this to a `data.frame` or `data.table`, you can set the `biobroom.return` option:

    options(biobroom.return = "data.frame")
    options(biobroom.return = "data.table")