Wrapper to use colocalization testing within a Bayesian model averaging structure.
Performs the colocalisation tests described in Plagnol et al (2009) and Wallace et al (2012).
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, ...)
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
|
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 |
bayes |
Logical, indicating whether to perform Bayesian
inference for the coefficient of proportionality, eta. If
|
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 |
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> |
quiet |
suppress messages about how the model spaced is trimmed for BMA |
bma |
if true (default), average over models |
... |
other parameters passed to |
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.
a coloc or colocBayes object
Chris Wallace
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
## 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)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.