From 0e50174c045dcbc82acfa98f3c6670d5990f895b Mon Sep 17 00:00:00 2001 From: chr1swallace Date: Sun, 30 Mar 2014 20:24:39 +0100 Subject: [PATCH 1/6] use new plot() method in coloc package --- R/coloc.mdf.R | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/R/coloc.mdf.R b/R/coloc.mdf.R index e5d2daf..90b5415 100755 --- a/R/coloc.mdf.R +++ b/R/coloc.mdf.R @@ -204,7 +204,8 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response), var.2 <- lapply(var.2, diag) if(plot.coeff) { - coloc:::coeff.plot(unlist(coef.1),unlist(coef.2), + message("Plots generated directly are deprecated. Please plot returned object instead.") + coloc:::coeff.plot(unlist(coef.1),unlist(coef.2), unlist(var.1),unlist(var.2), eta=stats["eta.hat"], main="Coefficients", @@ -214,11 +215,13 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response), if(!bayes) { return(new("coloc", result=c(stats["eta.hat"],chisquare=NA,stats["n"],stats["p"]), - method="BMA")) + method="BMA", + plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]))) } else { return(new("colocBayes", result=c(stats["eta.hat"],chisquare=NA,stats["n"],stats["p"]), method="BMA", + plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]), ppp=stats["ppp"], credible.interval=ci, bayes.factor=bf)) @@ -412,7 +415,7 @@ coloc.var.test <- function(X,vars.drop=NULL, ...) { ##' @param k Theta has a Cauchy(0,k) prior. The default, k=1, is equivalent to a ##' uniform (uninformative) prior. We have found varying k to have little effect ##' on the results. -##' @param plot.coeff \code{TRUE} if you want to generate a plot showing the +##' @param plot.coeff DEPRECATED. Please \code{plot()} returned object instead. \code{TRUE} if you want to generate a plot showing the ##' coefficients from the two regressions together with confidence regions. ##' @param bma parameter set to \code{TRUE} when \code{coloc.test} is called by \code{coloc.bma}. DO NOT SET THIS WHEN CALLING \code{coloc.test} DIRECTLY! ##' @param plots.extra list with 2 named elements, x and y, equal length @@ -435,7 +438,7 @@ coloc.var.test <- function(X,vars.drop=NULL, ...) { ##' approximating the posterior distribution at \code{n.approx} distinct values. ##' @param nsnps The SNPs to consider as potential explanatory variables ##' @author Mary Fortune -coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=TRUE,plots.extra=NULL,bayes=!is.null(bayes.factor), +coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=FALSE,plots.extra=NULL,bayes=!is.null(bayes.factor), n.approx=1001, level.ci=0.95, bayes.factor=NULL, bma=FALSE, nsnps=2) { snpnum <- length(b1) @@ -541,7 +544,8 @@ coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=TRUE,plots.extra=NULL, ## plots if(plot.coeff) { - coloc:::coeff.plot(b1,b2,diag(V11),diag(V22),eta=eta.hat, + message("Plots generated directly are deprecated. Please plot returned object instead.") + coloc:::coeff.plot(b1,b2,diag(V11),diag(V22),eta=eta.hat, main="Coefficients", # sub=paste("ppp =",format.pval(ppp$value,digits=2),"p =",format.pval(pchisq(X2,df=nsnps-1,lower.tail=FALSE),digits=2)), xlab=expression(b[1]),ylab=expression(b[2])) @@ -564,12 +568,14 @@ coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=TRUE,plots.extra=NULL, if(!bayes) { return(new("coloc", result=c(eta.hat=eta.hat,chisquare=X2,n=nsnps), - method="single")) + method="single", + plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]))) } else { if(!bma) { return(new("colocBayes", result=c(eta.hat=eta.hat,chisquare=X2,n=nsnps), method="single", + plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]), ppp=ppp$value, credible.interval=cred.int, bayes.factor=post.bf)) @@ -577,6 +583,7 @@ coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=TRUE,plots.extra=NULL, return(new("colocBayesBMA", result=c(eta.hat=eta.hat,chisquare=X2,n=nsnps), method="single", + plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]), ppp=ppp$value, bma=post.bma, bayes.factor=post.bf)) From 998cc1f0c2d162903c6bf0fa8dcd76c552b42295 Mon Sep 17 00:00:00 2001 From: Chris Wallace Date: Sun, 30 Mar 2014 20:30:32 +0100 Subject: [PATCH 2/6] Clarify final output from RA/T1D -> trait1/2 --- R/coloc.mdf.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/coloc.mdf.R b/R/coloc.mdf.R index 90b5415..4744237 100755 --- a/R/coloc.mdf.R +++ b/R/coloc.mdf.R @@ -668,8 +668,8 @@ singlesnp.twotrait<-function(df1,response="Y",snps=setdiff(colnames(df1),respons results.t2<-single.snp.tests(df.t2[,response], snp.data=snp.data.t2) capture.output(show.results.t2<-show(results.t2)) - cat("The minimum single SNP p value for RA is: ", min(show.results.t1[,4]),"\n") - cat("The minimum single SNP p value for T1D is: ", min(show.results.t2[,4]),"\n") + message("The minimum single SNP p value for trait 1 is: ", min(show.results.t1[,4])) + message("The minimum single SNP p value for trait 2 is: ", min(show.results.t2[,4])) return(c(min(show.results.t1[,4]),min(show.results.t2[,4]))) } From 96379807874617058f276123212c46d9ab700180 Mon Sep 17 00:00:00 2001 From: chr1swallace Date: Mon, 31 Mar 2014 11:37:05 +0100 Subject: [PATCH 3/6] use response instead of Y --- R/coloc.mdf.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/coloc.mdf.R b/R/coloc.mdf.R index 4744237..bd05b13 100755 --- a/R/coloc.mdf.R +++ b/R/coloc.mdf.R @@ -173,7 +173,7 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response), snps2<-paste("2:",unlist(modelsnps),sep="") if(!quiet) cat(".") - capture.output(multiglm <- multinom(as.formula(paste(response"Y~",paste(modelsnps,collapse="+"))), data=df1,maxit=1000)) + capture.output(multiglm <- multinom(as.formula(paste(response,"~",paste(modelsnps,collapse="+"))), data=df1,maxit=1000)) coef.1[[i]] <- coefficients(multiglm)[1,modelsnps] coef.2[[i]] <- coefficients(multiglm)[2,modelsnps] var[[i]] <- vcov(multiglm)[c(snps1,snps2),c(snps1,snps2)] From c27bfe43058dc3af869890ec1d629c2932adf7b2 Mon Sep 17 00:00:00 2001 From: chr1swallace Date: Mon, 31 Mar 2014 11:44:42 +0100 Subject: [PATCH 4/6] added dependency on mlogitBMA --- DESCRIPTION | 1 + 1 file changed, 1 insertion(+) diff --git a/DESCRIPTION b/DESCRIPTION index 8d53703..dfc97b7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,6 +4,7 @@ Depends: methods, MASS, BMA, + mlogitBMA, coloc Suggests: ggplot2, From 1921a159cbf2377e3bc089972d4a78bea4fa5c14 Mon Sep 17 00:00:00 2001 From: chr1swallace Date: Mon, 31 Mar 2014 16:06:37 +0100 Subject: [PATCH 5/6] minor changes + plot --- DESCRIPTION | 1 + R/coloc.mdf.R | 8 +++++--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index dfc97b7..a8623e0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -5,6 +5,7 @@ Depends: MASS, BMA, mlogitBMA, + nnet, coloc Suggests: ggplot2, diff --git a/R/coloc.mdf.R b/R/coloc.mdf.R index bd05b13..c871d57 100755 --- a/R/coloc.mdf.R +++ b/R/coloc.mdf.R @@ -31,8 +31,8 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response), n.orig <- length(snps) if(n.orig<3) return(1) - prep <- coloc:::prepare.df(df1, snps, r2.trim=r2.trim, dataset=1, quiet=quiet) - df1 <- prep$df + prep <- coloc:::prepare.df(df1=df1[,c(response,snps)], r2.trim=r2.trim, dataset=1, quiet=quiet) + df1 <- prep$df1[,c(response,prep$snps)] snps <- prep$snps if(!quiet) @@ -216,7 +216,9 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response), return(new("coloc", result=c(stats["eta.hat"],chisquare=NA,stats["n"],stats["p"]), method="BMA", - plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]))) + plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2), + var1=unlist(var.1),var2=unlist(var.2), + model.prob=probs[wh]))) } else { return(new("colocBayes", result=c(stats["eta.hat"],chisquare=NA,stats["n"],stats["p"]), From bdae3dd22cddea62c92738fc7aeae287cfe6756e Mon Sep 17 00:00:00 2001 From: Chris Wallace Date: Wed, 12 Aug 2015 11:52:40 +0100 Subject: [PATCH 6/6] Update README.md Updated Readme.md to point to paper --- README.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index dd4a97f..188539a 100755 --- a/README.md +++ b/README.md @@ -1,4 +1,7 @@ colocCommonControl ===== -Repo for the R package colocCommonControl \ No newline at end of file +Repo for the R package colocCommonControl + +For details on the method, and for citations, please see Fortune et al (2015) [http://www.nature.com/ng/journal/v47/n7/full/ng.3330.html]. Preprint is at bioRxiv [http://biorxiv.org/content/early/2015/06/08/020651]. +