Creating annotated output with and ReportingTools
26 July 2018
1 Overview
The package is intended to help people easily create useful output from various analyses. While affycoretools was originally intended for those using Affymetrix microarrays, this is no longer the case. While some functions remain Affy-centric, most are now much more general, and can be used for any microarray or RNA-Seq experiment.
2 Introduction
This package has evolved from my work as a service core biostatistician. I routinely analyze very similar experiments, and wanted to create a way to minimize all the cutting and pasting of code that I found myself doing. In addition, I wanted to come up with a good way to make an analysis reproducible, in the sense that I (or somebody else) could easily re-create the results.
In the past this package relied on the package, and was intended to be used in concert with a ‘Sweave’ document that contained both the code that was used to analyze the data, as well as explanatory text that would end up in a pdf (or HTML page) that could be given to a client. In the intervening period, people have developed other, better packages such as and that make it much easier to create the sort of output I like to present to my clients.
This document was generated using an Rmarkdown document (it’s the file called RefactoredAffycoretools.Rmd that you can find in your R library directory in the affycoretools/doc folder). Part of using Rmarkdown is generating .Rmd files that contain R code, and for each code block there are some arguments that define how the results for that code block are returned. A good resource for the code options is here.
3 Using affycoretools
For this section we will be using the data set that comes with the package. Remember that you can always run this code at home by doing this:
library(knitr)
purl(system.file("doc/RefactoredAffycoretools.Rmd", package="affycoretools"))And then you will have a file called RefactoredAffycoretools.R in your working directory that you can either or open with or , and run by chunk or line by line.
We first load and rename the data:
suppressMessages(library(affycoretools))
data(sample.ExpressionSet)
eset <- sample.ExpressionSet
eset## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 500 features, 26 samples
## element names: exprs, se.exprs
## protocolData: none
## phenoData
## sampleNames: A B ... Z (26 total)
## varLabels: sex type score
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation: hgu95av2This object is a truncated data set, based on an Affymetrix HG-U95av2 array. There are 26 samples and 500 probesets. We will use the to fit a linear model using .
Note that the limma package will retain annotation data that is in the , so we add those data in using
suppressMessages(library(hgu95av2.db))
eset <- annotateEset(eset, hgu95av2.db)## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columnssuppressMessages(library(limma))
pd <- pData(phenoData(eset))
design <- model.matrix(~0+type+sex, pd)
colnames(design) <- gsub("type|sex", "", colnames(design))
contrast <- matrix(c(1,-1,0))
colnames(contrast) <- "Case vs control"
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit)
topTable(fit2, 1)[,1:4]## PROBEID ENTREZID SYMBOL
## 31667_r_at 31667_r_at 10002 NR2E3
## AFFX-HSAC07_X00351_M_st AFFX-HSAC07_X00351_M_st <NA> <NA>
## 31375_at 31375_at <NA> <NA>
## 31466_at 31466_at 3128 HLA-DRB6
## 31597_r_at 31597_r_at 1978 EIF4EBP1
## 31440_at 31440_at 6932 TCF7
## 31396_r_at 31396_r_at 4440 MSI1
## AFFX-hum_alu_at AFFX-hum_alu_at <NA> <NA>
## 31391_at 31391_at 9001 HAP1
## AFFX-HSAC07_X00351_3_at AFFX-HSAC07_X00351_3_at <NA> <NA>
## GENENAME
## 31667_r_at nuclear receptor subfamily 2 group E member 3
## AFFX-HSAC07_X00351_M_st <NA>
## 31375_at <NA>
## 31466_at major histocompatibility complex, class II, DR beta 6 (pseudogene)
## 31597_r_at eukaryotic translation initiation factor 4E binding protein 1
## 31440_at transcription factor 7
## 31396_r_at musashi RNA binding protein 1
## AFFX-hum_alu_at <NA>
## 31391_at huntingtin associated protein 1
## AFFX-HSAC07_X00351_3_at <NA>After adding the annotation data to the object, the output now contains the appropriate annotation data for each probeset. At this point we can output an HTML table that contains these data.
suppressMessages(library(ReportingTools))
htab <- HTMLReport("afile", "My cool results")
publish(topTable(fit2, 1), htab)
finish(htab)## [1] "./afile.html"If you run the code yourself, you will have an HTML file ‘afile.html’ in the working directory, that contains the data for our top 10 genes. This table is not particularly interesting, and the package already has functionality to just do something like
htab <- HTMLReport("afile2", "My cool results, ReportingTools style")
publish(fit2, htab, eset, factor = pd$type, coef = 1, n = 10)
finish(htab)## [1] "./afile2.html"and it will automatically generate another HTML file, ‘afile2.html’, with some extra plots that show the different groups, and we didn’t even have to use directly. However, the default plots in the HTML table are a combination of dotplot and boxplot, which I find weird. Since is easily extensible, we can make changes that are more pleasing.
d.f <- topTable(fit2, 2)
out <- makeImages(df = d.f, eset = eset, grp.factor = pd$type, design = design,
contrast = contrast, colind = 1, repdir = ".")
htab <- HTMLReport("afile3", "My cool results, affycoretools style")
publish(out$df, htab, .mofifyDF = list(entrezLinks, affyLinks))
finish(htab)## [1] "./afile3.html"Note that there are two differences in the way we did things. First, we create a , and then decorate it with the plots using the function. This will by default create dotplots (or you can specify boxplots). For the plots to fit in an HTML table, there are no axis labels. However, each plot is also a link, and if you click on it, a larger plot with axis labels will be presented. If you are running the code yourself, see ‘afile3.html’.
All the little files that get created can get pretty messy, so the default is to put everything into a ‘reports’ subdirectory, so your working directory stays clean. For this example we over-ride the defaults so we do not have to go searching in subdirectories for our tables.
For HTML output the better idea is to generate a table using e.g., the xtable package, with links to the tables that we have generated.
An alternative parameterization that probably makes more sense is to fit coefficients for each sex/treatment combination.
grps <- factor(apply(pd[,1:2], 1, paste, collapse = "_"))
design <- model.matrix(~0+grps)
colnames(design) <- gsub("grps", "", colnames(design))
contrast <- matrix(c(1,-1,0,0,
0,0,1,-1,
1,-1,-1,1),
ncol = 3)
colnames(contrast) <- c("Female_Case vs Female_Control",
"Male_Case vs Male_Control",
"Interaction")
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)With this parameterization we can look at intra-sex differences, as well as the interaction (looking for sex-specific changes). This now means that we have a total of three HTML tables to output, which makes things a bit more complex to present. Luckily, this is pretty simple to accomplish. For this step we will now use the default ‘reports’ subdirectory to keep everything straight. In addition, we will trim down the output a bit.
## get a list containing the output for each comparison
out <- lapply(1:3, function(x) topTable(fit2, x))
## process the output to add images
htab <- lapply(1:3, function(x){
tmp <- HTMLReport(gsub("_", " ", colnames(contrast)[x]), colnames(contrast)[x], "./reports")
tmp2 <- makeImages(out[[x]], eset, grps, design, contrast, x)
publish(tmp2$df, tmp, .modifyDF = list(affyLinks, entrezLinks))
finish(tmp)
return(tmp)
})A reasonable thing to do would be to add a table to this document that lists all the different comparisons we made, with a count of the genes and links to the HTML tables we have generated. To do this we use the xtable package (you could also use as well), and add results = “asis” to the R code chunk in our Rmarkdown document.
d.f <- data.frame(Comparison = sapply(htab, function(x) XML::saveXML(Link(x))),
"Number significant" = sapply(out, nrow), check.names = FALSE)
kable(d.f, caption = "Number of significant genes.",
format = "html", row.names = FALSE)| Comparison | Number significant |
|---|---|
| <a href=“./reports/Female Case vs Female Control.html”>Female_Case vs Female_Control</a> | 10 |
| <a href=“./reports/Male Case vs Male Control.html”>Male_Case vs Male_Control</a> | 10 |
| <a href=“./reports/Interaction.html”>Interaction</a> | 10 |
We are often asked to create a Venn diagram showing overlap between groups. This is pretty simple to do, but it would be nicer to have an HTML version with clickable links, so your PI or end user can see what genes are in each cell of the Venn diagram. As an example, we can generate a Venn diagram comparing overlapping genes between the male and female comparisons.
collist <- list(1:2)
venn <- makeVenn(fit2, contrast, design, collist = collist, adj.meth = "none")
## generate a list of captions for each Venn.
## we only have one, so it's a list of length 1.
## we are using the BiocStyle package to control formatting,
## so the first sentence will be bolded and in blue.
cap <- list(paste("A Venn diagram. Note that the first sentence is bolded",
"and blue, whereas the rest is normal type and black. Usually",
"one would use a more useful caption than this one."))
vennInLine(venn, cap)
Figure 1: A Venn diagram
Note that the first sentence is bolded and blue, whereas the rest is normal type and black. Usually one would use a more useful caption than this one.
There is similar functionality for presenting the results of a GO hypergeometric analysis (), and GSEA analysis, based on the function in ( and ).
4 Session information
The version of R and packages loaded when creating this vignette were:
## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.7-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.7-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ReportingTools_2.20.0 knitr_1.20 limma_3.36.2
## [4] hgu95av2.db_3.2.3 org.Hs.eg.db_3.6.0 AnnotationDbi_1.42.1
## [7] IRanges_2.14.10 S4Vectors_0.18.3 affycoretools_1.52.2
## [10] Biobase_2.40.0 BiocGenerics_0.26.0 BiocStyle_2.8.2
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.2 GOstats_2.46.0
## [3] Hmisc_4.1-1 plyr_1.8.4
## [5] lazyeval_0.2.1 GSEABase_1.42.0
## [7] splines_3.5.1 BiocParallel_1.14.2
## [9] GenomeInfoDb_1.16.0 ggplot2_3.0.0
## [11] digest_0.6.15 foreach_1.4.4
## [13] BiocInstaller_1.30.0 ensembldb_2.4.1
## [15] htmltools_0.3.6 GO.db_3.6.0
## [17] gdata_2.18.0 magrittr_1.5
## [19] checkmate_1.8.5 memoise_1.1.0
## [21] BSgenome_1.48.0 cluster_2.0.7-1
## [23] gcrma_2.52.0 Biostrings_2.48.0
## [25] annotate_1.58.0 matrixStats_0.54.0
## [27] R.utils_2.6.0 ggbio_1.28.4
## [29] prettyunits_1.0.2 colorspace_1.3-2
## [31] blob_1.1.1 xfun_0.3
## [33] dplyr_0.7.6 crayon_1.3.4
## [35] RCurl_1.95-4.11 graph_1.58.0
## [37] genefilter_1.62.0 bindr_0.1.1
## [39] iterators_1.0.10 survival_2.42-6
## [41] VariantAnnotation_1.26.1 glue_1.3.0
## [43] gtable_0.2.0 zlibbioc_1.26.0
## [45] XVector_0.20.0 DelayedArray_0.6.2
## [47] Rgraphviz_2.24.0 scales_0.5.0
## [49] DBI_1.0.0 GGally_1.4.0
## [51] edgeR_3.22.3 Rcpp_0.12.18
## [53] xtable_1.8-2 progress_1.2.0
## [55] htmlTable_1.12 foreign_0.8-71
## [57] bit_1.1-14 OrganismDbi_1.22.0
## [59] preprocessCore_1.42.0 Formula_1.2-3
## [61] AnnotationForge_1.22.1 htmlwidgets_1.2
## [63] httr_1.3.1 gplots_3.0.1
## [65] RColorBrewer_1.1-2 acepack_1.4.1
## [67] ff_2.2-14 pkgconfig_2.0.1
## [69] reshape_0.8.7 XML_3.98-1.12
## [71] R.methodsS3_1.7.1 nnet_7.3-12
## [73] locfit_1.5-9.1 labeling_0.3
## [75] tidyselect_0.2.4 rlang_0.2.1
## [77] reshape2_1.4.3 munsell_0.5.0
## [79] tools_3.5.1 RSQLite_2.1.1
## [81] evaluate_0.11 stringr_1.3.1
## [83] yaml_2.2.0 bit64_0.9-7
## [85] oligoClasses_1.42.0 caTools_1.17.1.1
## [87] purrr_0.2.5 AnnotationFilter_1.4.0
## [89] bindrcpp_0.2.2 RBGL_1.56.0
## [91] R.oo_1.22.0 biomaRt_2.36.1
## [93] compiler_3.5.1 rstudioapi_0.7
## [95] curl_3.2 affyio_1.50.0
## [97] PFAM.db_3.6.0 tibble_1.4.2
## [99] geneplotter_1.58.0 stringi_1.2.4
## [101] highr_0.7 GenomicFeatures_1.32.0
## [103] lattice_0.20-35 ProtGenerics_1.12.0
## [105] Matrix_1.2-14 pillar_1.3.0
## [107] data.table_1.11.4 bitops_1.0-6
## [109] rtracklayer_1.40.3 GenomicRanges_1.32.6
## [111] R6_2.2.2 latticeExtra_0.6-28
## [113] affy_1.58.0 hwriter_1.3.2
## [115] bookdown_0.7 KernSmooth_2.23-15
## [117] gridExtra_2.3 codetools_0.2-15
## [119] dichromat_2.0-0 gtools_3.8.1
## [121] assertthat_0.2.0 SummarizedExperiment_1.10.1
## [123] DESeq2_1.20.0 Category_2.46.0
## [125] rprojroot_1.3-2 GenomicAlignments_1.16.0
## [127] Rsamtools_1.32.2 GenomeInfoDbData_1.1.0
## [129] hms_0.4.2 grid_3.5.1
## [131] rpart_4.1-13 rmarkdown_1.10
## [133] biovizBase_1.28.1 base64enc_0.1-3