ammegandchips/bacon

Mmm (vegan) bacon... controlling bias and inflation in association studies using the empirical null distribution

bacon

Mmm (vegan) bacon... controlling bias and inflation in association studies using the empirical null distribution

bacon can be used to remove inflation and bias often observed in epigenome- and transcriptome-wide association studies (Iterson, Zwet, and Heijmans 2017). You can read more about it here: https://bioconductor.org/packages/release/bioc/vignettes/bacon/inst/doc/bacon.html

Install and load bacon

#try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("bacon")
require(bacon)

When you just have one set of test statistics

This doesn't work if you try to run it on test statistics from METAL. Very small effects and SEs are rounded to 0, so we can't calculate a z score :(

z <- dat$beta / dat$se
bc <- bacon(z)
bc
estimates(bc)
inflation(bc)
bias(bc)

Histogram of Z-scores - with standard normal (black) and estimated empirical null distribution (red)

plot(bc, type="hist")

QQ plots of -log10 transformed P-values - left panel using uncorrected P-values and right panel using bacon bias and inflation corrected P-values

plot(bc, type="qq")

When you have multiple sets of test statistics (e.g. from several cohorts)

Youll probably want to run this in parallel

library(BiocParallel)
register(MulticoreParam(1, log=TRUE))

es is a matrix of effect sizes (beta coefficients) and se is a matrix of the equivalent standard errors. Each row is a probe and each column is a study.

bc <- bacon(NULL, es, se)
bc
estimates(bc)
inflation(bc)
bias(bc)
plot(bc, type="hist")
plot(bc, type="qq")

When you have multiple sets of test statistics and you want to do a fixed-effect meta-analysis

You'll probably want to run this in parallel

library(BiocParallel)
register(MulticoreParam(1, log=TRUE))
bc <- bacon(NULL, es, se)
bc
bcm <- meta(bc)
head(pval(bcm))
print(topTable(bcm))
plot(bcm, type="qq")

An example using the PACE BMI results

load("/Volumes/filestore/Documents/IEU work/PACE/Maternal BMI/Results/Meta-analysis/BMI/BMI_1_random1.Rdata")
es<-na.omit(as.matrix(M_C[,grep(colnames(M_C),pattern="_beta")]))
se<-na.omit(as.matrix(M_C[,grep(colnames(M_C),pattern="_se")]))

bacon doesn't appear to cope with NAs, and since not every cohort supplied results for every probe, there are a lot of NAs. In fact there are only 255,514 probes

How have other people got round this?

bacon's Issues

Twoi questions for paper "Controlling bias and inflation in epigenome-and transcriptome-wide association studies using the empirical null distribution

I still encounter some problems when use bacon packages, and I would like to ask you for advice. Here are two problems:

1. In our study, lambdas (λ) by conventional approach were 1.06( A cohort) to 1.07 (B cohort), and lambdas based on the bacon approach (λ.bacon) were 0.974( A cohort) to 0.982( B cohort).In particular, P value is more significant after performing bacon method while lambda is smaller in A cohort( λ.bacon 0.974), which made me confuse.
2. λ.bacon of A cohort is 0.97, λ.bacon of A cohort is 0.98, while there are "inconsistence" between the "QQ plots" and " λ.bacon " (this QQ plots was generated by bacon packages). Red is A cohort, green is B cohort. Inflation of A cohort seems to be higher than that of cohort B.
   Could you give me some advice about that? Thanks!
   Looking forward to your reply and Wish you everything.

qq plot of meta-analysis using bacon

When I have two sets of test statistics and you want to do a fixed-effect meta-analysis:
library(BiocParallel)
register(MulticoreParam(1, log=TRUE))
bc <- bacon(NULL, es, se)
bc
bcm <- meta(bc)
head(pval(bcm))
print(topTable(bcm))
plot(bcm, type="qq")

When I use above code graph qq plot of meta-analysis, there are two qq line in the plot (following figure). Most importantly, My aim was to graph a qq plot with only one qq line which was used to access the similarity between distribution of bacon adjusted p values of meta analysis and empirical null distribution, but I don't know how to do that using ggplo2 packages in R. Can you give some advice?