CPBayes (Bayesian meta analysis for studying cross-phenotype genetic associations) package

Arunabha Majumdar, Tanushree Haldar, John Witte

2017-04-05

Introduction

Simultaneous analysis of genetic associations with multiple phenotypes may reveal shared genetic susceptibility across traits (pleiotropy). CPBayes is a Bayesian meta analysis method for studying cross-phenotype genetic associations. It uses summary-level data across multiple phenotypes to simultaneously measure the evidence of aggregate-level pleiotropic association and estimate an optimal subset of traits associated with the risk locus. CPBayes is based on a spike and slab prior and is implemented by Markov chain Monte Carlo (MCMC) technique Gibbs sampling.

This R-package consists of five main functions:

  1. cpbayes_uncor: It implements CPBayes for uncorrelated summary statistics. The summary statistics across traits/studies are uncorrelated when the studies have no overlapping subject.
  2. cpbayes_cor: It implements CPBayes for correlated summary statistics. The summary statistics across traits/studies are correlated when the studies have overlapping subjects or the phenotypes were measured in a cohort study.
  3. post_summaries: It summarizes the MCMC data produced by cpbayes_uncor or cpbayes_cor. It computes additional summaries to provide a better insight into a pleiotropic signal. It works in the same way for both cpbayes_uncor and cpbayes_cor.
  4. forest_cpbayes: It creates a forest plot presenting the pleiotropy result obtained by cpbayes_uncor or cpbayes_cor. It works in the same way for both cpbayes_uncor and cpbayes_cor.
  5. estimate_corln: It computes an approximate correlation matrix of the beta-hat vector for multiple overlapping case-control studies using the sample-overlap matrices.

Installation

You can install CPBayes from CRAN.

install.packages("CPBayes")
library("CPBayes")

How to run CPBayes for uncorrelated summary statistics.

Get the path to the data.

library("CPBayes")
# Load the beta hat vector
BetaHatfile <- system.file("extdata", "BetaHat.rda", package = "CPBayes")
load(BetaHatfile)
BetaHat
##  [1]  0.03402285 -0.01987781 -0.01088476 -0.00940975  0.02195188
##  [6] -0.01728955  0.14180403  0.06979527 -0.11721340 -0.09231879

BetaHat contains an example data of the main genetic effect (beta/log(odds ratio)) estimates for a single nucleotide polymorphism (SNP) obtained from 10 separate case-control studies for 10 different diseases. In each case-control study comprising a distinct set of 7000 cases and 10000 controls, we fit a logistic regression of the case-control status on the genotype coded as the minor allele count for all the individuals in the sample. One can also include various covariates, such as, age, gender, principal components (PCs) of ancestries in the logistic regression. From each logistic regression for a disease, we obtain the estimate of the main genetic association parameter (beta/log(odds ratio)) along with the corresponding standard error. Since the studies do not have any overlapping subject, beta-hat across the diseases are uncorrelated.

# Load the standard error vector
SEfile <- system.file("extdata", "SE.rda", package = "CPBayes")
load(SEfile)
SE
##  [1] 0.02736746 0.02741876 0.02749038 0.02747637 0.02749800 0.02753650
##  [7] 0.02731022 0.02714103 0.02786046 0.02783711

SE contains the standard errors corresponding to the above beta hat vector across 10 separate case-control studies.

Next we specify the name of the diseases/phenotypes and the genetic variant.

# Specify the name of the traits and the genetic variant.
traitNames <- paste("Disease", 1:10, sep = "")
SNP1 <- "rs1234"
traitNames
##  [1] "Disease1"  "Disease2"  "Disease3"  "Disease4"  "Disease5" 
##  [6] "Disease6"  "Disease7"  "Disease8"  "Disease9"  "Disease10"
SNP1
## [1] "rs1234"

Now we implement CPBayes for this example data. Since the studies are non-overlapping, the summary statistics across traits are uncorrelated. Hence we run the the cpbayes_uncor function.

# Run the uncorrelated version of CPBayes.
result <- cpbayes_uncor(BetaHat, SE, Phenotypes = traitNames, Variant = SNP1)
## log10_BF: 11.39 
## subset:  Disease7 Disease9 Disease10

There are more options of arguments to pass into the function (see the Arguments section of cpbayes_uncor in the CPBayes manual). After running cpbayes_uncor, it prints the log10(Bayes factor) (denoted as log10_BF) and the subset of associated/non-null traits (denoted as subset) produced by CPBayes. The Bayes factor evaluates the overall pleiotropic association and the subset of non-null traits are the most important phenotypes that underlie the pleiotropic signal. However, the printed outputs are only a part of ‘result’ which is a list that constitutes of various components. An overall summary of ‘result’ can be seen by using the str() function (as shown below).

# Overall summary of the primary results produced by cpbayes_uncor.
str(result)
## List of 7
##  $ variantName     : chr "rs1234"
##  $ log10_BF        : num 11.4
##  $ PPNA            : num 3.99e-15
##  $ subset          : chr [1:3] "Disease7" "Disease9" "Disease10"
##  $ important_traits:'data.frame':    3 obs. of  2 variables:
##   ..$ traits: chr [1:3] "Disease7" "Disease9" "Disease10"
##   ..$ PPAj  : num [1:3] 1 0.962 0.6
##  $ auxi_data       :List of 8
##   ..$ traitNames     : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ K              : int 10
##   ..$ mcmc.samplesize: num 10000
##   ..$ PPAj           : num [1:10] 0.029 0.0191 0.0214 0.0159 0.0232 ...
##   ..$ Z.data         : num [1:10000, 1:10] 0 0 0 0 0 0 0 0 0 1 ...
##   ..$ sim.beta       : num [1:10000, 1:10] -0.000319 0.010188 0.012115 -0.003118 0.003075 ...
##   ..$ betahat        : num [1:10] 0.03402 -0.01988 -0.01088 -0.00941 0.02195 ...
##   ..$ se             : num [1:10] 0.0274 0.0274 0.0275 0.0275 0.0275 ...
##  $ runtime         :Class 'proc_time'  Named num [1:5] 1.804 0.014 1.822 0 0
##   .. ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ...

result$variantName returns the name of the genetic variant specified by the user. Here, it is ‘rs1234’. result$log10_BF provides the log10(Bayes factor) (discussed above). result$PPNA provides the posterior probability of null association which is another measure of the aggregate-level pleiotropic association. result$subset provides the optimal subset of associated/non-null traits selected by CPBayes. Thus, for the example data, CPBayes produced a log10(Bayes factor) of 11.4, PPNA of 3.99*10^(-15), and selected Disease7, Disease9, and Disease10 as associated/non-null. A detailed interpretation of all the outputs are described in the Value section of cpbayes_uncor in the CPBayes manual.

The post_summaries function provides important insights into an obseved pleiotropic signal, e.g., the direction of associations, trait-specific posterior probability of associations (PPAj), posterior mean/median and 95% credible interval (Bayesian analog of the confidence interval) of the unknown true genetic effect (beta/odds ratio) on each trait, etc.

# Post summary of the MCMC data produced by cpbayes_uncor.
PleioSumm <- post_summaries(result, level = 0.05)  
str(PleioSumm)
## List of 9
##  $ variantName       : chr "rs1234"
##  $ log10_BF          : num 11.4
##  $ PPNA              : num 3.99e-15
##  $ subset            :'data.frame':  3 obs. of  3 variables:
##   ..$ traits   : chr [1:3] "Disease7" "Disease9" "Disease10"
##   ..$ PPAj     : num [1:3] 1 0.962 0.6
##   ..$ direction: chr [1:3] "positive" "negative" "negative"
##  $ important_traits  :'data.frame':  3 obs. of  3 variables:
##   ..$ traits   : chr [1:3] "Disease7" "Disease9" "Disease10"
##   ..$ PPAj     : num [1:3] 1 0.962 0.6
##   ..$ direction: chr [1:3] "positive" "negative" "negative"
##  $ traitNames        : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##  $ PPAj              :'data.frame':  10 obs. of  2 variables:
##   ..$ traits: chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ PPAj  : num [1:10] 0.029 0.0191 0.0214 0.0159 0.0232 ...
##  $ poste_summary_beta:'data.frame':  10 obs. of  6 variables:
##   ..$ traits      : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ poste_mean  : num [1:10] 0.00498 -0.00271 -0.00136 -0.00128 0.00278 ...
##   ..$ poste_median: num [1:10] 0.0043 -0.00237 -0.0012 -0.00124 0.00259 ...
##   ..$ poste_se    : num [1:10] 0.0115 0.0103 0.0104 0.0101 0.0104 ...
##   ..$ lCl         : num [1:10] -0.014 -0.0225 -0.0212 -0.0205 -0.0161 ...
##   ..$ uCl         : num [1:10] 0.0266 0.0162 0.0177 0.0178 0.0225 ...
##  $ poste_summary_OR  :'data.frame':  10 obs. of  5 variables:
##   ..$ traits      : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ poste_mean  : num [1:10] 1.005 0.997 0.999 0.999 1.003 ...
##   ..$ poste_median: num [1:10] 1.004 0.998 0.999 0.999 1.003 ...
##   ..$ lCl         : num [1:10] 0.986 0.978 0.979 0.98 0.984 ...
##   ..$ uCl         : num [1:10] 1.03 1.02 1.02 1.02 1.02 ...

So we have to pass the list ‘result’ returned by cpbayes_uncor as the first argument and the ‘level’ as the second argument into the post_summaries function. If ‘level’ is not specified, the default value is 0.05. Note that post_summaries computes (1-level)% credible interval of the unknown true genetic effect (beta/odds ratio) on each trait. It estimates the direction of association with the important traits, the vector of trait-specific posterior probability of association (PPAj), etc. For detailed description of different outputs provided by this function, please see the Value section of post_summaries in the CPBayes manual.

Next we run the forest_cpbayes function to create a forest plot that presents the pleiotropy result produced by cpbayes_uncor.

# Forest plot for the pleiotropy result obtained by cpbayes_uncor.
forest_cpbayes(result, level = 0.05)
## quartz_off_screen 
##                 2

Similarly as for the post_summaries function, we need to pass the same list `result’ returned by cpbayes_uncor as the first argument into the function. Second argument is the level whose default value is 0.05. In the forest plot, (1-level)% confidence interval of the beta/log(odds ratio) parameter is plotted for each trait. For more details, please see the section of forest_cpbayes function in the CPBayes manual.

How to run CPBayes for correlated summary statistics.

Next we demonstrate how to run CPBayes for correlated summary statistics. Get the path to the data.

# Load the beta-hat vector
datafile <- system.file("extdata", "cBetaHat.rda", package = "CPBayes")
load(datafile)
cBetaHat
##  [1] -0.042216290 -0.035624411 -0.058063641 -0.026037581 -0.014434543
##  [6] -0.004717336 -0.043771014 -0.025293739  0.059845798 -0.177518725

Here ‘c’ in cBetaHat stands for correlated case. cBetaHat contains an example data of the main genetic association parameter (beta/log(odds ratio)) estimates for a SNP across 10 overlapping case-control studies for 10 different diseases. Each of the 10 studies has a distinct set of 7000 cases and a common set of 10000 controls shared across all the studies. In each case-control study, we fit a logistic regression of the case-control status on the genotype coded as the minor allele count for all the individuals in the sample. One can also include various covariates, such as, age, gender, principal components (PCs) of ancestries in the logistic regression. From each logistic regression for a disease, we obtain the estimate of the main genetic effect (beta/log(odds ratio)) along with the corresponding standard error. Since the studies have overlapping subjects, beta-hat across the diseases are correlated.

# Load the standard error vector
datafile <- system.file("extdata", "cSE.rda", package = "CPBayes")
load(datafile)
cSE
##  [1] 0.02386863 0.02389768 0.02404311 0.02395724 0.02383974 0.02376651
##  [7] 0.02402079 0.02392416 0.02377087 0.02446383

cSE contains the standard errors corresponding to the above beta hat vector across 10 overlapping case-control studies.

# Load the correlation matrix of the beta-hat vector (cBetaHat)
datafile <- system.file("extdata", "cor.rda", package = "CPBayes")
load(datafile)
cor
##            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
##  [1,] 1.0000000 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647
##  [2,] 0.4117647 1.0000000 0.4117647 0.4117647 0.4117647 0.4117647
##  [3,] 0.4117647 0.4117647 1.0000000 0.4117647 0.4117647 0.4117647
##  [4,] 0.4117647 0.4117647 0.4117647 1.0000000 0.4117647 0.4117647
##  [5,] 0.4117647 0.4117647 0.4117647 0.4117647 1.0000000 0.4117647
##  [6,] 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647 1.0000000
##  [7,] 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647
##  [8,] 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647
##  [9,] 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647
## [10,] 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647 0.4117647
##            [,7]      [,8]      [,9]     [,10]
##  [1,] 0.4117647 0.4117647 0.4117647 0.4117647
##  [2,] 0.4117647 0.4117647 0.4117647 0.4117647
##  [3,] 0.4117647 0.4117647 0.4117647 0.4117647
##  [4,] 0.4117647 0.4117647 0.4117647 0.4117647
##  [5,] 0.4117647 0.4117647 0.4117647 0.4117647
##  [6,] 0.4117647 0.4117647 0.4117647 0.4117647
##  [7,] 1.0000000 0.4117647 0.4117647 0.4117647
##  [8,] 0.4117647 1.0000000 0.4117647 0.4117647
##  [9,] 0.4117647 0.4117647 1.0000000 0.4117647
## [10,] 0.4117647 0.4117647 0.4117647 1.0000000

The correlation matrix of the beta-hat vector (cBetaHat) is given by ‘cor’ which we estimated by employing the estimate_corln function (demonstrated later in this tutorial) using the sample-overlap matrices (explained later in this tutorial). Next we run the correlated version of CPBayes for this example data.

# Run the correlated version of CPBayes.
result <- cpbayes_cor(cBetaHat, cSE, cor, Phenotypes = traitNames, Variant = SNP1)
## log10_BF: 18.87 
## subset:  Disease9 Disease10

There are more options of arguments to pass into the function (see the Arguments section of cpbayes_cor in the CPBayes manual). After running cpbayes_cor, it prints the log10(Bayes factor) (denoted as log10_BF) and the subset of non-null/associated traits (denoted as subset) produced by CPBayes. The Bayes factor evaluates the overall pleiotropic association and the subset of non-null traits are the most important phenotypes that underlie the pleiotropic signal. However, the printed outputs are only a part of ‘result’ which is a list that constitutes of various components. An overall summary of ‘result’ can be seen by using the str() function (as shown below).

# Overall summary of the primary results produced by cpbayes_cor.
str(result)
## List of 8
##  $ variantName     : chr "rs1234"
##  $ log10_BF        : num 18.9
##  $ PPNA            : num 1.31e-22
##  $ subset          : chr [1:2] "Disease9" "Disease10"
##  $ important_traits:'data.frame':    2 obs. of  2 variables:
##   ..$ traits: chr [1:2] "Disease9" "Disease10"
##   ..$ PPAj  : num [1:2] 0.91 1
##  $ auxi_data       :List of 8
##   ..$ traitNames     : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ K              : int 10
##   ..$ mcmc.samplesize: num 10000
##   ..$ PPAj           : num [1:10] 0.0109 0.0105 0.0369 0.0069 0.0119 ...
##   ..$ Z.data         : num [1:10000, 1:10] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$ sim.beta       : num [1:10000, 1:10] -0.00604 -0.01855 0.00852 0.00693 -0.00578 ...
##   ..$ betahat        : num [1:10] -0.0422 -0.0356 -0.0581 -0.026 -0.0144 ...
##   ..$ se             : num [1:10] 0.0239 0.0239 0.024 0.024 0.0238 ...
##  $ uncor_use       : chr "No"
##  $ runtime         :Class 'proc_time'  Named num [1:5] 4.53 0.047 4.581 0 0
##   .. ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ...

result$variantName returns the name of the genetic variant specified by the user. Here, it is ‘rs1234’. result$log10_BF provides the log10(Bayes factor) (discussed above). result$PPNA provides the posterior probability of null association which is another measure of the aggregate-level pleiotropic association. result$subset provides the optimal subset of associated traits selected by CPBayes. Thus, for the example data, CPBayes produced a log10(Bayes factor) of 18.9, PPNA of 1.31*10^(-22), and selected Disease9 and Disease10 as associated/non-null. A detailed interpretation of all the outputs are described in the Value section of cpbayes_cor in the CPBayes manual.

The post_summaries function provides important insights into an observed pleiotropic signal, e.g., the direction of associations, trait-specific posterior probability of associations (PPAj), posterior mean/median and 95% credible interval (Bayesian analog of the confidence interval) of the unknown true genetic effect (beta/odds ratio) on each trait, etc.

# Post summary of the MCMC data produced by cpbayes_cor.
PleioSumm <- post_summaries(result, level = 0.05)  
str(PleioSumm)
## List of 9
##  $ variantName       : chr "rs1234"
##  $ log10_BF          : num 18.9
##  $ PPNA              : num 1.31e-22
##  $ subset            :'data.frame':  2 obs. of  3 variables:
##   ..$ traits   : chr [1:2] "Disease9" "Disease10"
##   ..$ PPAj     : num [1:2] 0.91 1
##   ..$ direction: chr [1:2] "positive" "negative"
##  $ important_traits  :'data.frame':  2 obs. of  3 variables:
##   ..$ traits   : chr [1:2] "Disease9" "Disease10"
##   ..$ PPAj     : num [1:2] 0.91 1
##   ..$ direction: chr [1:2] "positive" "negative"
##  $ traitNames        : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##  $ PPAj              :'data.frame':  10 obs. of  2 variables:
##   ..$ traits: chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ PPAj  : num [1:10] 0.0109 0.0105 0.0369 0.0069 0.0119 ...
##  $ poste_summary_beta:'data.frame':  10 obs. of  6 variables:
##   ..$ traits      : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ poste_mean  : num [1:10] -0.004227 -0.002783 -0.008707 -0.000548 0.002431 ...
##   ..$ poste_median: num [1:10] -0.004274 -0.002593 -0.007977 -0.000616 0.002421 ...
##   ..$ poste_se    : num [1:10] 0.00932 0.00931 0.01096 0.00902 0.00935 ...
##   ..$ lCl         : num [1:10] -0.0222 -0.021 -0.0326 -0.0182 -0.0157 ...
##   ..$ uCl         : num [1:10] 0.01351 0.01485 0.00943 0.01721 0.02038 ...
##  $ poste_summary_OR  :'data.frame':  10 obs. of  5 variables:
##   ..$ traits      : chr [1:10] "Disease1" "Disease2" "Disease3" "Disease4" ...
##   ..$ poste_mean  : num [1:10] 0.996 0.997 0.991 0.999 1.002 ...
##   ..$ poste_median: num [1:10] 0.996 0.997 0.992 0.999 1.002 ...
##   ..$ lCl         : num [1:10] 0.978 0.979 0.968 0.982 0.984 ...
##   ..$ uCl         : num [1:10] 1.01 1.01 1.01 1.02 1.02 ...

Note that, post_summaries works exactly in the same way for both cpbayes_cor and cpbayes_uncor. Thus, we need to pass the list ‘result’ returned by cpbayes_cor as the first argument and ‘level’ as the second argument. If ‘level’ is not specified, the default value is 0.05. Note that post_summaries computes (1-level)% credible interval of the unknown true genetic effect (beta/odds ratio) on each trait. It provides the direction of association with the important traits, the vector of trait-specific posterior probability of association (PPAj), etc. For detailed description of different outputs provided by post_summaries, please see the Value section of post_summaries in the CPBayes manual.

Next we run the forest_cpbayes function to create a forest plot that presents the pleiotropy result produced by cpbayes_cor.

# Forest plot for the pleiotropy result obtained by cpbayes_cor.
forest_cpbayes(result, level = 0.05)
## quartz_off_screen 
##                 2

Note that, forest_cpbayes works exactly in the same way for both cpbayes_cor and cpbayes_uncor. Thus, we need to pass the list `result’ returned by cpbayes_cor as the first argument into the function. Second argument is the level for which the default value is 0.05. In the forest plot, (1-level)% confidence interval of the beta/log(odds ratio) parameter is plotted for each trait. For more details, please see the section of forest_cpbayes function in the CPBayes manual.

How to run estimate_corln.

The function estimate_corln estimates the correlation matrix of the beta-hat vector for multiple overlapping case-control studies using the sample-overlap matrices which describe the number of cases or controls shared between studies/traits, and the number of subjects who are case for one study/trait but control for another study/trait. For a cohort study, the phenotypic correlation matrix should be a reasonable substitute of this correlation matrix.

# Example data of sample-overlap matrices
SampleOverlapMatrixFile <- system.file("extdata", "SampleOverlapMatrix.rda", package = "CPBayes")
load(SampleOverlapMatrixFile)
SampleOverlapMatrix
## $n11
##        trait1 trait2 trait3 trait4 trait5
## trait1   9048   4647   2985   2835   1812
## trait2   4647  13565   3873   4245   2419
## trait3   2985   3873  14681   6285   2044
## trait4   2835   4245   6285  16697   2059
## trait5   1812   2419   2044   2059   7121
## 
## $n00
##        trait1 trait2 trait3 trait4 trait5
## trait1  44683  35765  32987  30821  39374
## trait2  35765  40166  29358  27714  35464
## trait3  32987  29358  39050  28638  33973
## trait4  30821  27714  28638  37034  31972
## trait5  39374  35464  33973  31972  46610
## 
## $n10
##        trait1 trait2 trait3 trait4 trait5
## trait1      0   4401   6063   6213   7236
## trait2   8918      0   9692   9320  11146
## trait3  11696  10808      0   8396  12637
## trait4  13862  12452  10412      0  14638
## trait5   5309   4702   5077   5062      0

SampleOverlapMatrix is a list that contains an example of the sample overlap matrices for five different diseases in the Kaiser GERA cohort (a real data). The list constitutes of three matrices as follows. SampleOverlapMatrix$n11 provides the number of cases shared between all possible pairs of studies/traits. SampleOverlapMatrix$n00 provides the number of controls shared between all possible pairs of studies/traits. SampleOverlapMatrix$n10 provides the number of subjects who are case for one study/trait and control for another study/trait. For more detailed explanation, see the Arguments section of estimate_corln in the CPBayes manual.

# Estimate the correlation matrix of correlated beta-hat vector
n11 <- SampleOverlapMatrix$n11
n00 <- SampleOverlapMatrix$n00
n10 <- SampleOverlapMatrix$n10
cor <- estimate_corln(n11, n00, n10)
cor
##             trait1      trait2     trait3       trait4      trait5
## trait1 1.000000000 0.270490702 0.05723195  0.002505875  0.08989408
## trait2 0.270490702 1.000000000 0.01601813  0.002744961  0.07849158
## trait3 0.057231953 0.016018131 1.00000000  0.155476792  0.01211052
## trait4 0.002505875 0.002744961 0.15547679  1.000000000 -0.01824859
## trait5 0.089894085 0.078491585 0.01211052 -0.018248589  1.00000000

The function estimate_corln computes an approximate correlation matrix of the correlated beta-hat vector obtained from multiple overlapping case-control studies using the sample-overlap matrices. Note that for a cohort study, the phenotypic correlation matrix should be a reasonable substitute of this correlation matrix. These approximations of the correlation structure are accurate when none of the diseases/traits is associated with the environmental covariates and genetic variant. While demonstrating cpbayes_cor, we used simulated data for 10 overlapping case-control studies with each study having a distinct set of 7000 cases and a common set of 10000 controls shared across all the studies. We used the estimate_corln function to estimate the correlation matrix of the correlated beta-hat vector using the sample-overlap matrices.

Important note on the estimation of correlation structure of correlated beta-hat vector: In general, environmental covariates are expected to be present in a study and associated with the phenotypes of interest. Also, a small proportion of genome-wide genetic variants are expected to be associated. Hence the above approximations of the correlation matrix may not be accurate. So in general, we recommend an alternative strategy to estimate the correlation matrix using the genome-wide summary statistics data across traits as follows. First, extract all the SNPs for each of which the trait-specific univariate association p-value across all the traits are > 0.1. The trait-specific univariate association p-values are obtained using the beta-hat and standard error for each trait. Each of the SNPs selected in this way is either weakly or not associated with any of the phenotypes (null SNP). Next, select a set of independent null SNPs from the initial set of null SNPs by using a threshold of r^2 < 0.01 (r: the correlation between the genotypes at a pair of SNPs). In the absence of in-sample linkage disequilibrium (LD) information, one can use the reference panel LD information for this screening. Finally, compute the correlation matrix of the effect estimates (beta-hat vector) as the sample correlation matrix of the beta-hat vector across all the selected independent null SNPs. This strategy is more general and applicable to a cohort study or multiple overlapping studies for binary or quantitative traits with arbitrary distributions. It is also useful when the beta-hat vector for multiple non-overlapping studies become correlated due to genetically related individuals across studies. Misspecification of the correlation structure can affect the results produced by CPBayes to some extent. Hence, if genome-wide summary statistics data across traits is available, we highly recommend to use this alternative strategy to estimate the correlation matrix of the beta-hat vector.

See our paper for more details: Arunabha Majumdar, Tanushree Haldar, Sourabh Bhattacharya, John Witte. An efficient Bayesian meta-analysis approach for studying cross-phenotype genetic associations (submitted), available at: http://biorxiv.org/content/early/2017/01/18/101543.