Contents
Package: bacon
Authors: Maarten van Iterson [aut, cre], Erik van Zwet [ctb]
Modified: Fri Mar 25 11:54:05 2016
Compiled: Wed Aug 24 22:50:37 2016
Introduction
bacon can be used to remove inflation and bias often observed in epigenome- and transcriptome-wide association studies(Iterson et al. 2016).
To this end bacon constructs an empirical null distribution using a Gibbs Sampling algorithm by fitting a three-component normal mixture on z-scores. One component is forced, using prior knowledge, to represent the null distribution with mean and standard deviation representing the bias and inflation. The other two components are necessary to capture the amount of true associations present in the data, which we assume unknown but small.
bacon provides functionality to inspect the output of the Gibbs Sampling algorithm, i.e., traces, posterior distributions and the mixture fit, are provided. Furthermore, inflation- and bias-corrected test-statistics or P-values are extracted easily. In addition, functionality for performing fixed-effect meta-analysis are provided.
The function bacon requires a vector or a matrix of z-scores, e.g., those extracted from association analyses using a linear regression approach. For fixed-effect meta-analysis a matrix of effect-sizes and standard errors is required.
The vignette illustrates the use of bacon using simulated z-scores, effect-sizes and standard errors to avoid long runtimes. If multiple sets of test-statisics or effect-sizes and standard errors are provided, the Gibbs Sampler Algorithm can be executed on multiple nodes to reduce computation time using functionality provide by BiocParallel-package.
One set of test-statistics
A vector containing 2000 z-scores is generated from a normal mixture distribution 90% of the z-scores were drawn from \(\mathcal{N}(0,1)\) and the remaining z-scores from \(\mathcal{N}(\mu,1)\), where \(\mu \sim \mathcal{N}(4,1)\).
The function bacon executes the Gibbs Sampler Algorithm and stores all input and output in a Bacon-object. Several accessor-functions are available to access data contained in the Bacon-object.
set.seed(12345)
y <- rnormmix(2000, c(0.9, 0, 1, 0, 4, 1))
bc <- bacon(y)
bc
## Bacon-object containing 1 set(s) of 2000 test-statistics.
## ...estimated bias: -0.0086.
## ...estimated inflation: 0.97.
##
## Emprical null estimates are based on 5000 iterations with a burnin-period of 2000.
estimates(bc)
## p.0 p.1 p.2 mu.0 mu.1 mu.2
## [1,] 0.9085228 0.04003456 0.05144263 -0.008585455 3.022833 -3.034213
## sigma.0 sigma.1 sigma.2
## [1,] 0.9729677 2.438307 2.789957
inflation(bc)
## sigma.0
## 0.9729677
bias(bc)
## mu.0
## -0.008585455
Several methods are provided to inspect the output of the Gibbs Sampler Algorithm, such as traces of all estimates, the posterior distributions, provide as a scatter plot between two parameters, the actual fit of the three component mixture to the histogram of z-scores.
traces(bc, burnin=FALSE)
posteriors(bc)
fit(bc, n=100)
There is also a generic plot function that can generate two types of plots. A histogram of the z-scores with on top the standard normal distribution and the bacon estimated empirical null distribution. Or a quantile-quantile plot of the \(-log_{10}\) transformed P-values.
plot(bc, type="hist")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 40 rows containing non-finite values (stat_bin).
## Warning: Removed 1 rows containing missing values (geom_bar).
## Warning: Removed 40 rows containing missing values (geom_path).
## Warning: Removed 40 rows containing missing values (geom_path).
plot(bc, type="qq")
Multiple sets of test-statistics
Matrices containing \(2000\times6\) effect-sizes and standard errors are generated to simulated data for a fixed-effect meta-analyses.
By default the function bacon detects the number of cores/nodes registered, as described in the BiocParallel, to perform bacon in parallel. To run the vignette in general we set it here for convenience to 1 node.
set.seed(12345)
es <- replicate(6, rnormmix(2000, c(0.9, 0, 1, 0, 4, 1)))
se <- replicate(6, 0.8*sqrt(4/rchisq(2000,df=4)))
colnames(es) <- colnames(se) <- LETTERS[1:ncol(se)]
rownames(es) <- rownames(se) <- 1:2000
head(rownames(es))
## [1] "1" "2" "3" "4" "5" "6"
head(colnames(es))
## [1] "A" "B" "C" "D" "E" "F"
library(BiocParallel)
register(MulticoreParam(1, log=TRUE))
bc <- bacon(NULL, es, se)
## Did you registered a biocparallel back-end?
## Continuing serial!
bc
## Bacon-object containing 6 set(s) of 2000 test-statistics.
## ...estimated bias: -0.0072,-0.02,0.00073,0.027,-0.017,-0.012.
## ...estimated inflation: 1.1,1.1,1.1,1.1,1.1,1.
##
## Emprical null estimates are based on 5000 iterations with a burnin-period of 2000.
estimates(bc)
## p.0 p.1 p.2 mu.0 mu.1 mu.2
## A 0.8798367 0.04596230 0.07420096 -0.0071923032 2.798940 -2.721537
## B 0.8748018 0.03514827 0.09004994 -0.0195838309 2.950528 -2.346152
## C 0.8764848 0.05558323 0.06793199 0.0007300818 2.945500 -2.924680
## D 0.8881519 0.08040985 0.03143825 0.0265941513 1.782361 -2.994063
## E 0.8626044 0.06886578 0.06852984 -0.0165292415 2.726810 -2.704635
## F 0.8556689 0.08026559 0.06406549 -0.0119121793 2.732789 -2.770299
## sigma.0 sigma.1 sigma.2
## A 1.112915 3.314581 4.3679065
## B 1.086432 1.678806 4.6548036
## C 1.071829 3.214286 3.1438970
## D 1.123699 4.807509 0.7085753
## E 1.062347 3.584511 3.4328507
## F 1.037343 3.763897 3.2223009
inflation(bc)
## A B C D E F
## 1.112915 1.086432 1.071829 1.123699 1.062347 1.037343
bias(bc)
## A B C D E
## -0.0071923032 -0.0195838309 0.0007300818 0.0265941513 -0.0165292415
## F
## -0.0119121793
head(tstat(bc))
## A B C D E F
## 1 1.6655396 -1.0920842 -1.27039727 -3.54419344 0.22730396 -0.04754207
## 2 -0.9763604 -1.3866084 0.26764794 0.25283954 0.32951747 0.35145701
## 3 -2.3990357 -3.6003210 0.06117898 -0.88181612 -0.14455838 0.05857065
## 4 -0.6691733 -0.5607974 0.53520168 0.03226887 -0.43766649 1.02497293
## 5 0.2316802 1.2621708 -0.90449418 1.01593285 -0.03631256 -0.08534002
## 6 0.9341088 -2.1999929 0.98707127 1.06925040 0.36631468 1.86003051
head(pval(bc))
## A B C D E F
## 1 0.09580515 0.2747960882 0.2039432 0.0003938162 0.8201874 0.96208121
## 2 0.32888589 0.1655612093 0.7889703 0.8003922129 0.7417646 0.72524552
## 3 0.01643831 0.0003178246 0.9512167 0.3778762551 0.8850596 0.95329409
## 4 0.50338494 0.5749356410 0.5925104 0.9742576369 0.6616281 0.30537596
## 5 0.81678638 0.2068873436 0.3657334 0.3096613667 0.9710331 0.93199108
## 6 0.35024777 0.0278074002 0.3236077 0.2849568533 0.7141303 0.06288121
head(se(bc))
## A B C D E F
## 1 0.9405472 1.3831921 1.2506331 0.8857125 0.9575235 0.5686023
## 2 0.5689726 0.6145793 0.9507878 2.2527801 1.3183756 1.3418700
## 3 2.3020580 0.5050518 0.9296272 1.2309362 0.7592934 0.8475108
## 4 0.6104344 0.7547990 1.2170681 1.4928418 1.2937036 0.5585942
## 5 0.7090617 0.6803576 0.6539196 1.3638215 0.9466390 0.7315781
## 6 0.6987288 0.6218310 0.8613305 1.9490522 1.0951570 0.6209050
head(es(bc))
## A B C D E F
## 1 1.5665186 -1.5105622 -1.58880094 -3.13913641 0.21764888 -0.02703253
## 2 -0.5555224 -0.8521808 0.25447639 0.56959188 0.43442779 0.47160962
## 3 -5.5227195 -1.8183488 0.05687364 -1.08545937 -0.10976221 0.04963926
## 4 -0.4084864 -0.4232893 0.65137691 0.04817231 -0.56621074 0.57254397
## 5 0.1642756 0.8587275 -0.59146649 1.38555106 -0.03437488 -0.06243289
## 6 0.6526887 -1.3680237 0.85019464 2.08402484 0.40117207 1.15490232
The accessor-function return as expected matrices of estimates. For the plotting functions an additional index of the ith study or z-score is required.
traces(bc, burnin=FALSE, index=3)
posteriors(bc, index=3)
fit(bc, index=3, n=100)
plot(bc, type="hist")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 240 rows containing non-finite values (stat_bin).
## Warning: Removed 6 rows containing missing values (geom_bar).
## Warning: Removed 43 rows containing missing values (geom_path).
## Warning: Removed 43 rows containing missing values (geom_path).
plot(bc, type="qq")
Session Info
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 LTS
##
## 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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bacon_1.0.5 ellipse_0.3-8 BiocParallel_1.6.6
## [4] ggplot2_2.1.0 BiocStyle_2.0.3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.6 digest_0.6.10 plyr_1.8.4 grid_3.3.1
## [5] gtable_0.2.0 formatR_1.4 magrittr_1.5 evaluate_0.9
## [9] scales_0.4.0 stringi_1.1.1 rmarkdown_1.0 labeling_0.3
## [13] tools_3.3.1 stringr_1.1.0 munsell_0.4.3 parallel_3.3.1
## [17] yaml_2.1.13 colorspace_1.2-6 htmltools_0.3.5 knitr_1.14
References
Iterson, M van, E van Zwet, P Slagboom, and B.T Heijmans. 2016. “Controlling Bias and Inflation in Association Studies Using the Empirical Null Distribution.” Manuscript in Preparation.
