coloc (version 3.2-1)

coloc.bma: Wrapper to use colocalization testing within a Bayesian model averaging structure.

Description

Performs the colocalisation tests described in Plagnol et al (2009) and Wallace et al (2012).

Usage

coloc.bma(df1, df2, snps = intersect(setdiff(colnames(df1), c(response1,
  stratum1)), setdiff(colnames(df2), c(response2, stratum2))),
  response1 = "Y", response2 = "Y", stratum1 = NULL,
  stratum2 = NULL, family1 = "binomial", family2 = "binomial",
  bayes = !is.null(bayes.factor), thr = 0.01, nsnps = 2,
  n.approx = 1001, bayes.factor = NULL, plot.coeff = FALSE,
  r2.trim = 0.95, quiet = FALSE, bma = FALSE, ...)

Arguments

df1, df2

Each is a dataframe, containing response and potential explanatory variables for two independent datasets.

snps

The SNPs to consider as potential explanatory variables

response1, response2

The names of the response variables in df1 and df2 respectively

stratum1

optional column name of df1 that gives stratum information

stratum2

optional column name of df2 that gives stratum information

family1, family2

the error family for use in glm

bayes

Logical, indicating whether to perform Bayesian inference for the coefficient of proportionality, eta. If bayes.factor is supplied, Bayes factors are additionally computed for the specificed values. This can add a little time as it requires numerical integration, so can be set to FALSE to save time in simulations, for example.

thr

posterior probability threshold used to trim SNP list. Only SNPs with a marginal posterior probability of inclusion greater than this with one or other trait will be included in the full BMA analysis

nsnps

number of SNPs required to model both traits. The BMA analysis will average over all possible nsnp SNP models, subject to thr above.

n.approx

number of values at which to numerically approximate the posterior

bayes.factor

if true, compare specific models

plot.coeff

deprecated

r2.trim

for pairs SNPs with r2>r2.trim, only one SNP will be retained. This avoids numerical instability problems caused by including two highly correlated SNPs in the model.

quiet

suppress messages about how the model spaced is trimmed for BMA

bma

if true (default), average over models

...

other parameters passed to coloc.test

Value

a coloc or colocBayes object

Details

This is a test for proportionality of regression coefficients from two independent regressions. Analysis can either be based on a profile likelihood approach, where the proportionality coefficient, eta, is replaced by its maximum likelihood value, and inference is based on a chisquare test (p.value), or taking a hybrid-Bayesian approach and integrating the p value over the posterior distribution of eta, which gives a posterior predictive p value. The Bayesian approach can also be used to give a credible interval for eta. See the references below for further details.

References

Wallace et al (2012). Statistical colocalisation of monocyte gene expression and genetic risk variants for type 1 diabetes. Hum Mol Genet 21:2815-2824. http://europepmc.org/abstract/MED/22403184

Plagnol et al (2009). Statistical independence of the colocalized association signals for type 1 diabetes and RPS26 gene expression on chromosome 12q13. Biostatistics 10:327-34. http://www.ncbi.nlm.nih.gov/pubmed/19039033

Examples

# NOT RUN {
 ## simulate covariate matrix (X) and continuous response vector (Y)
 ## for two populations/triats Y1 and Y2 depend equally on f1 and f2
 ## within each population, although their distributions differ between
 ## populations.  They are compatible with a null hypothesis that they
 ## share a common causal variant
set.seed(1)
 X1 <- matrix(rbinom(2000,1,0.4),ncol=4)
 Y1 <- rnorm(500,rowSums(X1[,1:2]),2)
 X2 <- matrix(rbinom(2000,1,0.6),ncol=4)
 Y2 <- rnorm(500,rowSums(X2[,1:2]),5)
 
 boxplot(list(Y1,Y2),names=c("Y1","Y2"))
 
 ## fit and store linear model objects
 colnames(X1) <- colnames(X2) <- sprintf("f%s",1:ncol(X1))
 summary(lm1 <- lm(Y1~f1+f2+f3+f4,data=as.data.frame(X1)))
 summary(lm2 <- lm(Y2~f1+f2+f3+f4,data=as.data.frame(X2)))
 

## test colocalisation using bma
df1=as.data.frame(cbind(Y1=Y1,X1))
df2=as.data.frame(cbind(Y2=Y2,X2))

result <- coloc.bma( df1, df2, snps=colnames(X1), response1="Y1", response2="Y2",
family1="gaussian", family2="gaussian",
nsnps=2,bayes.factor=c(1,2,3))
result
plot(result)

## test colocalisation when one dataset contains a stratifying factor in column named "s"
df1$s <- rbinom(500,1,0.5)
result <- coloc.bma( df1, df2, snps=colnames(X1), response1="Y1", response2="Y2",
stratum1="s",
family1="gaussian", family2="gaussian",
nsnps=2,bayes.factor=c(1,2,3))
result
plot(result)

# }