R/sherlock1.R

In jean997/sherlockAsh: CAUSE: Causal Analysis Using Summary Effect Estimates

#'@title Calculate Bayes Factor for original Sherlock
#'@param eqtl_pval Vector of p-values for gene expression (length p)
#'@param gwas_pval Vector of p-values for gwas trait (length p)
#'@param maf Minor allele frequences (length p)
#'@param n1,n2 sample sizes for eqtl and gwas studies
#'@param alpah,beta,sigma_a_eqtl,sigma_a_pheno Tuning parameters (see Xin et al 2011).
#'Default values are alpha=5e-5, beta=5e-4, sigma_a_eqtl = 0.5, sigam_a_pheno=0.5
#'@return log Bayes factor
#'@export
sherlock1 <- function(eqtl_pval, gwas_pval, maf, n1, n2,
                      alpha=5e-5, beta=5e-4, sigma_a_eqtl=0.5, sigma_a_pheno=0.5){
  p <- length(eqtl_pval)
  stopifnot(length(gwas_pval)==p & length(maf)==p)
  stopifnot(length(alpha) %in% c(1, p))

  #Calculate B_ix
  abs_zscore_eqtl <- qnorm(eqtl_pval/2)
  var_u1 <- 1 + 2*n1*maf*(1-maf)*sigma_a_eqtl^2
  B_ix <- dnorm(abs_zscore_eqtl, mean=0, sd=sqrt(var_u1))/dnorm(abs_zscore_eqtl)

  #Calculate B_iy
  abs_zscore_gwas <- qnorm(gwas_pval/2)
  var_v1 <- 1 + 2*n2*maf*(1-maf)*sigma_a_pheno^2
  B_iy <- dnorm(abs_zscore_gwas, mean=0, sd=sqrt(var_v1))/dnorm(abs_zscore_gwas)

  #Calculate Bayes factor
  pi_x <- (alpha*B_ix)/(1 - alpha + alpha*B_ix)

  B_i <- (1-pi_x) + (pi_x)*(B_iy/(1-beta+beta*B_iy))

  return(sum(log(B_i)))
}