Convert Bayes factors to beta and standard error

statistics
fine mapping
Author

Gibran Hemani

Published

January 26, 2023

Background

Is it possible to convert BF to beta and standard error? According to Giambartolomei et al 2014 -

\[ ABF = \sqrt{1-r} \times exp(rZ^2/2) \]

so

\[ |Z| = \sqrt{\frac{2 * log(ABF) - log(\sqrt{1-r})}{r}} \]

here \(r = W / V\) where V is the variance of the SNP effect estimate

\[ V \approx \frac{1}{2np(1-p)} \]

where n is sample size and p is allele frequency (assumes small amount of variance explained in trait and sd of trait is 1).

Run simulation

  1. Use regional LD matrix to generate summary statistics with a single causal variant
  2. Use SuSiE to perform fine mapping
  3. Convert SuSiE Bayes Factors into Z scores, betas, standard errors
  4. Compare converted Z, beta, se against original simulated Z, beta, SE

Simulation

Libraries

library(simulateGP)
library(susieR)
library(here)
here() starts at /Users/gh13047/repo/lab-book
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Conversion function for logBF to z, beta, se

#' Convert log Bayes Factor to summary stats
#'
#' @param lbf p-vector of log Bayes Factors for each SNP
#' @param n Overall sample size
#' @param af p-vector of allele frequencies for each SNP
#' @param prior_v Variance of prior distribution. SuSiE uses 50
#'
#' @return tibble with lbf, af, beta, se, z
lbf_to_z_cont <- function(lbf, n, af, prior_v=50)
{
  se = sqrt(1 / (2 * n * af * (1-af)))
  r = prior_v / (prior_v + se^2)
  z = sqrt((2 * lbf - log(sqrt(1-r)))/r)
  beta <- z * se
  return(tibble(lbf, af, z, beta, se))
}

Read in example LD matrix from simulateGP repository

map <- readRDS(url("https://github.com/explodecomputer/simulateGP/raw/master/data/ldobj_5_141345062_141478055.rds", "rb"))
glimpse(map)
List of 3
 $ ld  : num [1:501, 1:501] 1 0.565 0.566 0.565 0.565 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:501] "V2000" "V2001" "V2002" "V2003" ...
 $ map : tibble [501 × 6] (S3: tbl_df/tbl/data.frame)
  ..$ chr: int [1:501] 5 5 5 5 5 5 5 5 5 5 ...
  ..$ snp: chr [1:501] "rs252141" "rs252140" "rs252139" "rs187544" ...
  ..$ pos: int [1:501] 141345062 141345192 141345218 141345361 141345678 141345805 141346830 141347360 141347465 141347931 ...
  ..$ alt: chr [1:501] "T" "T" "C" "G" ...
  ..$ ref: chr [1:501] "C" "C" "T" "T" ...
  ..$ af : num [1:501] 0.627 0.831 0.83 0.831 0.831 ...
  ..- attr(*, ".internal.selfref")=<externalptr> 
 $ nref: num 503

Generate summary statistics for a single causal variant and

set.seed(1234)
ss <- map$map %>%
    generate_gwas_params(h2=0.003, Pi=1/nrow(.)) %>%
    generate_gwas_ss(50000, ld=map$ld)
table(ss$beta == 0)

FALSE  TRUE 
    1   500 
plot(-log10(pval) ~ pos, ss)

Run SuSiE

sout <- susie_rss(ss$bhat / ss$se, R = map$ld, n = 50000, bhat = ss$bhat, var_y=1)
WARNING: XtX is not symmetric; forcing XtX to be symmetric by replacing XtX with (XtX + t(XtX))/2
summary(sout)

Variables in credible sets:

 variable variable_prob cs
      286     0.1604616  1
      306     0.1604616  1
      291     0.1604616  1
      300     0.1604616  1
      274     0.1604616  1
      284     0.1604616  1

Credible sets summary:

 cs cs_log10bf cs_avg_r2 cs_min_r2                variable
  1   30.37357         1         1 274,284,286,291,300,306
glimpse(sout)
List of 18
 $ alpha                 : num [1:10, 1:501] 6.9e-35 2.0e-03 2.0e-03 2.0e-03 2.0e-03 ...
 $ mu                    : num [1:10, 1:501] 0.000761 0 0 0 0 ...
 $ mu2                   : num [1:10, 1:501] 2.04e-05 0.00 0.00 0.00 0.00 ...
 $ KL                    : num [1:10] 6.75 -1.24e-14 -1.24e-14 -1.24e-14 -1.24e-14 ...
 $ lbf                   : num [1:10] 6.99e+01 1.24e-14 1.24e-14 1.24e-14 1.24e-14 ...
 $ lbf_variable          : num [1:10, 1:501] -2.51 0 0 0 0 ...
 $ sigma2                : num 1
 $ V                     : num [1:10] 0.00307 0 0 0 0 ...
 $ pi                    : num [1:501] 0.002 0.002 0.002 0.002 0.002 ...
 $ null_index            : num 0
 $ XtXr                  : num [1:501, 1] -0.328 70.558 72.085 70.558 70.558 ...
 $ converged             : logi TRUE
 $ elbo                  : num [1:2] -70876 -70876
 $ niter                 : int 2
 $ X_column_scale_factors: num [1:501] 1 1 1 1 1 1 1 1 1 1 ...
 $ intercept             : num NA
 $ sets                  :List of 5
  ..$ cs                :List of 1
  .. ..$ L1: int [1:6] 274 284 286 291 300 306
  ..$ purity            :'data.frame':  1 obs. of  3 variables:
  .. ..$ min.abs.corr   : num 1
  .. ..$ mean.abs.corr  : num 1
  .. ..$ median.abs.corr: num 1
  ..$ cs_index          : int 1
  ..$ coverage          : num 0.963
  ..$ requested_coverage: num 0.95
 $ pip                   : num [1:501] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "class")= chr "susie"

Get Z scores from lbf

a <- lbf_to_z_cont(sout$lbf_variable[1,], 50000, ss$af, prior_v = 50)
a
# A tibble: 501 × 5
     lbf     af     z    beta      se
   <dbl>  <dbl> <dbl>   <dbl>   <dbl>
 1 -2.51 0.373   1.41 0.00919 0.00654
 2 -2.43 0.169   1.37 0.0115  0.00844
 3 -2.37 0.17    1.42 0.0119  0.00842
 4 -2.43 0.169   1.37 0.0115  0.00844
 5 -2.43 0.169   1.37 0.0115  0.00844
 6 -2.52 0.191   1.32 0.0106  0.00805
 7 -2.43 0.169   1.37 0.0115  0.00844
 8 -2.44 0.17    1.36 0.0115  0.00842
 9  2.23 0.0139  3.17 0.0855  0.0270 
10 -2.43 0.169   1.37 0.0115  0.00844
# … with 491 more rows

Relationship between lbf and re-estimated z

plot(z ~ lbf, a)

New Z vs original Z

plot(a$z^2 ~ ss$fval)

lm(a$z^2 ~ ss$fval)

Call:
lm(formula = a$z^2 ~ ss$fval)

Coefficients:
(Intercept)      ss$fval  
     1.5141       0.9834

New beta vs original beta

plot(a$beta ~ ss$bhat)

Two causal variants

Set two causal variants at either end of the region

set.seed(12)
param <- map$map
param$beta <- 0
param$beta[c(10, 490)] <- 0.3
ss <- generate_gwas_ss(param, 50000, ld=map$ld)
plot(-log10(pval) ~ pos, ss)

First variant

sout <- susie_rss(ss$bhat / ss$se, R = map$ld, n = 50000, bhat = ss$bhat, var_y=1)
WARNING: XtX is not symmetric; forcing XtX to be symmetric by replacing XtX with (XtX + t(XtX))/2
a1 <- lbf_to_z_cont(sout$lbf_variable[1,], 50000, ss$af, prior_v = 50)
Warning in sqrt((2 * lbf - log(sqrt(1 - r)))/r): NaNs produced
plot(a1$beta ~ ss$bhat)

a2 <- lbf_to_z_cont(sout$lbf_variable[2,], 50000, ss$af, prior_v = 50)
plot(a2$beta ~ ss$bhat)

This looks good - it’s setting different values to 0 in the two lbf vectors that correspond to two causal variants

sessionInfo()
R version 4.2.1 Patched (2022-09-06 r82817)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.6.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dplyr_1.0.10     here_1.0.1       susieR_0.12.27   simulateGP_0.1.2

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9         plyr_1.8.7         compiler_4.2.1     pillar_1.8.1      
 [5] tools_4.2.1        digest_0.6.31      jsonlite_1.8.4     evaluate_0.19     
 [9] lifecycle_1.0.3    tibble_3.1.8       gtable_0.3.1       lattice_0.20-45   
[13] pkgconfig_2.0.3    rlang_1.0.6        Matrix_1.4-1       DBI_1.1.3         
[17] cli_3.5.0          yaml_2.3.6         xfun_0.36          fastmap_1.1.0     
[21] stringr_1.5.0      knitr_1.41         generics_0.1.3     vctrs_0.5.1       
[25] htmlwidgets_1.5.4  rprojroot_2.0.3    tidyselect_1.2.0   grid_4.2.1        
[29] reshape_0.8.9      glue_1.6.2         R6_2.5.1           fansi_1.0.3       
[33] rmarkdown_2.16     mixsqp_0.3-48      irlba_2.3.5.1      ggplot2_3.4.0     
[37] magrittr_2.0.3     MASS_7.3-58.1      matrixStats_0.63.0 scales_1.2.1      
[41] htmltools_0.5.4    assertthat_0.2.1   colorspace_2.0-3   utf8_1.2.2        
[45] stringi_1.7.8      munsell_0.5.0      crayon_1.5.2