Data transforms
-Jeff Leek
-Dependencies
-This document depends on the following packages:
- library(devtools)
- library(Biobase)
- library(dplyr)To install these packages you can use the code (or if you are compiling the document, remove the eval=FALSE from the chunk.)
install.packages(c("devtools"))
-source("http://www.bioconductor.org/biocLite.R")
-biocLite(c("Biobase"))General principles
--
-
- Make sure data is on the right scale for -
- Visualization -
- Statistical modeling -
Load some data
-We will use this expression set to look at how we use plots and tables to check for different characteristics
-con =url("http://bowtie-bio.sourceforge.net/recount/ExpressionSets/bodymap_eset.RData")
-load(file=con)
-close(con)
-bm = bodymap.eset
-pdata=pData(bm)
-edata=as.data.frame(exprs(bm))
-fdata = fData(bm)
-ls()## [1] "bm" "bodymap.eset" "con" "edata"
-## [5] "fdata" "pdata" "tropical"Skewed distributions
-We would like continuous data to be nice and symmetric like a normal distribution for two reasons:(1) plots are easier to see this way and (2) most statistical methods are designed to work better for non-skewed data.
-hist(rnorm(1000),col=2)
Realistically though most genomic data is skewed
-hist(edata[,1],col=2,breaks=100)
One way to address this skew is to use a transformation. One common transformation is the log transform.
hist(log(edata[,1]),col=2,breaks=100)
One thing to be careful with is that values of zero become -inf when you apply the log transform because the log of zero isn’t defined.
min(log(edata))## [1] -InfYou can remove this problem by adding a small number to all the counts before taking the log. This doesn’t change the overall distribution much but resolves the zero problem
-min(log(edata[,1] + 1))## [1] 0hist(log(edata[,1] + 1),breaks=100,col=2)
A common choice is to use the log2 transform. The reason is that data on the log2 scale can be expressed as fold changes. If log2(x) - log2(y) = z then there is a z fold change difference between x and y.
hist(log2(edata[,1] + 1),breaks=100,col=2)
Let’s zoom in on the values greater than zero (we’ll come back to the zero values in a minute)
-hist(log2(edata[,1] + 1),breaks=100,col=2,xlim=c(1,15),ylim=c(0,400))
This is starting to look a little better but we still have the zeros problem.
-Removing features without much data
-In many types of sequencing data we observe that there is not much data for many features.
-hist(rowSums(edata==0),col=2)
A common pre-processing technique is to remove features that don’t have much data. One way to identify and remove them is to use rowMeans and rowMedians
low_genes = rowMeans(edata) < 5
-table(low_genes)## low_genes
-## FALSE TRUE
-## 11413 41167filt_edata = filter(edata,!low_genes)
-dim(filt_edata)## [1] 11413 19low_genes2 = rowMedians(as.matrix(edata)) < 5
-table(low_genes2,low_genes)## low_genes
-## low_genes2 FALSE TRUE
-## FALSE 9844 2
-## TRUE 1569 41165filt_edata2 = filter(edata,!low_genes2)
-dim(filt_edata2)## [1] 9846 19After filtering we observe that using the log transform makes the distribution much more symmetric but there are still zeros. We could increase the filter to remove this issue if we needed to.
hist(log2(filt_edata[,1] + 1),col=2)
More transforms
-There are a large number of other transforms depending on the type of data you are using. Some common ones are:
--
-
- Variance stabilizing transforms which seek to remove a mean variance relationship among the data -
- Box-Cox transforms which seek to make the data approximately Normally distributed -
- rlog transform - unique to genomics count data, this is a regularized version of the log transform that seeks to minimize differences at low count levels. -
Session information
-Here is the session information
-devtools::session_info()## setting value
-## version R version 3.2.2 (2015-08-14)
-## system x86_64, darwin13.4.0
-## ui X11
-## language (EN)
-## collate en_GB.UTF-8
-## tz Europe/London
-## date 2015-09-17
-##
-## package * version date source
-## assertthat 0.1 2013-12-06 CRAN (R 3.2.0)
-## Biobase * 2.28.0 2015-04-17 Bioconductor
-## BiocGenerics * 0.14.0 2015-04-17 Bioconductor
-## DBI 0.3.1 2014-09-24 CRAN (R 3.2.0)
-## devtools * 1.9.1 2015-09-11 CRAN (R 3.2.0)
-## digest 0.6.8 2014-12-31 CRAN (R 3.2.0)
-## dplyr * 0.4.3 2015-09-01 CRAN (R 3.2.2)
-## evaluate 0.7.2 2015-08-13 CRAN (R 3.2.0)
-## formatR 1.2 2015-04-21 CRAN (R 3.2.0)
-## htmltools 0.2.6 2014-09-08 CRAN (R 3.2.0)
-## knitr * 1.11 2015-08-14 CRAN (R 3.2.2)
-## lazyeval 0.1.10 2015-01-02 CRAN (R 3.2.0)
-## magrittr 1.5 2014-11-22 CRAN (R 3.2.0)
-## memoise 0.2.1 2014-04-22 CRAN (R 3.2.0)
-## R6 2.1.1 2015-08-19 CRAN (R 3.2.2)
-## Rcpp 0.12.1 2015-09-10 CRAN (R 3.2.0)
-## rmarkdown 0.8 2015-08-30 CRAN (R 3.2.2)
-## stringi 0.5-5 2015-06-29 CRAN (R 3.2.0)
-## stringr 1.0.0 2015-04-30 CRAN (R 3.2.0)
-## yaml 2.1.13 2014-06-12 CRAN (R 3.2.0)It is also useful to compile the time the document was processed. This document was processed on: 2015-09-17.
-
