Last updated: 2022-06-10

Checks: 7 0

Knit directory: finemap/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it's best to always run the code in an empty environment.

Seed: set.seed(1)

The command set.seed(1) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 38a8ed5

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 38a8ed5. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 


Ignored files:
    Ignored:    data/small_data_11.ld_refout_file.refout.ld
    Ignored:    data/small_data_11_sim_gaussian_pve_n_8_get_sumstats_n_1.ld_sample_n_file.in_n.ld
    Ignored:    output/sim1.config
    Ignored:    output/sim2.config

Unstaged changes:
    Modified:   output/sim1.cred1
    Modified:   output/sim1.cred2
    Modified:   output/sim1.log_sss
    Modified:   output/sim1.snp
    Modified:   output/sim2.bf5
    Modified:   output/sim2.cred5
    Modified:   output/sim2.log_sss
    Modified:   output/sim2.snp

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/large_effect.Rmd) and HTML (docs/large_effect.html) files. If you've configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 38a8ed5 Peter Carbonetto 2022-06-10 workflowr::wflow_publish("large_effect.Rmd", view = FALSE, verbose = TRUE)
html e597078 Peter Carbonetto 2022-06-10 Added finemap results to large_effect example.
Rmd c3f0aa3 Peter Carbonetto 2022-06-10 workflowr::wflow_publish("large_effect.Rmd", view = FALSE, verbose = TRUE)
Rmd 6c863a1 Peter Carbonetto 2022-06-10 Working on adding finemap analysis to large_effect example.
html f40a337 Peter Carbonetto 2022-06-09 Adjusted plot dims slightly.
Rmd 9dc1f67 Peter Carbonetto 2022-06-09 workflowr::wflow_publish("large_effect.Rmd", verbose = TRUE)
html 37627a7 Peter Carbonetto 2022-06-09 Added susie-rss fit using out-of-sample LD.
Rmd 2c327c8 Peter Carbonetto 2022-06-09 workflowr::wflow_publish("large_effect.Rmd", verbose = TRUE)
Rmd 7af3892 Peter Carbonetto 2022-06-08 Working on large_effect demo.
html cd367c3 Peter Carbonetto 2022-06-08 Added step to run susie with in-sample ld matrix.
Rmd 1f607ce Peter Carbonetto 2022-06-08 workflowr::wflow_publish("large_effect.Rmd")
Rmd f3b3f26 Peter Carbonetto 2022-06-08 More improvements to large_effect demo.
Rmd ba1823b Peter Carbonetto 2022-06-08 Working on large_effect demo.
html ba1823b Peter Carbonetto 2022-06-08 Working on large_effect demo.
html e68ebf0 Peter Carbonetto 2022-06-07 Build site.
Rmd 5ba80d9 Peter Carbonetto 2022-06-07 Small edit.
html 5737016 Peter Carbonetto 2022-06-07 First build of large_effect analysis.
Rmd 4b1ca7f Peter Carbonetto 2022-06-07 workflowr::wflow_publish("large_effect.Rmd")
Rmd eae747f Peter Carbonetto 2022-06-07 Revised workflowr files.

In this small example drawn from our simulations, we show that that FINEMAP works well with an "in-sample LD" matrix—that is, a correlation matrix that was estimated using the same sample that was used to compute the single-SNP association statistics—but, can perform surprisingly poorly with an "out-of-sample" LD matrix. We have observed that this degradation in performance is unusual—specifically, when the effects of the causal SNPs are very large. Nonetheless, this example may be instructive. In contrast to FINEMAP, SuSiE performs similarly well in this example with either the in-sample and out-of-sample LD matrix.

First, load some packages, and set the seed for reproducibility.

library(data.table)
library(susieR)
set.seed(1)

Load the summary data: the least-squares effect estimates \(\hat{\beta}_i\) and their standard errors \(\hat{s}_i\) for each SNP \(i\). Here we also compute the z-scores since SuSiE accepts z-scores as input.

dat1 <- readRDS("../data/small_data_11.rds")
dat3 <- readRDS("../data/small_data_11_sim_gaussian_pve_n_8_get_sumstats_n_1.rds")
maf  <- dat1$maf$in_sample
bhat <- dat3$sumstats$bhat
shat <- dat3$sumstats$shat
z    <- bhat/shat

In this simulation, 2 out of 1,001 SNPs affect the phenotype:

dat2 <- readRDS("../data/small_data_11_sim_gaussian_pve_n_8.rds")
b <- drop(dat2$meta$true_coef)
vars <- which(b != 0)
vars
# [1] 305 740

In-sample LD

Run SuSiE with the "in-sample" LD estimate:

ldinfile <- "../data/small_data_11_sim_gaussian_pve_n_8_get_sumstats_n_1.ld_sample_n_file.in_n.ld"
Rin <- as.matrix(fread(ldinfile))
fit1 <- susie_rss(z,Rin,n = 800,min_abs_corr = 0.1,refine = FALSE,
                 verbose = TRUE)
# HINT: If the in-sample LD matrix is available, we recommend calling susie_rss with the in-sample LD matrix, and setting estimate_residual_variance = TRUE
# HINT: For large R or large XtX, consider installing theRfast package for better performance.
# WARNING: XtX is not symmetric; forcing XtX to be symmetric by replacing XtX with (XtX + t(XtX))/2
# [1] "objective:-1022.11016014142"
# [1] "objective:-1022.09485818746"
# [1] "objective:-1022.08938881341"
# [1] "objective:-1022.08631623581"
# [1] "objective:-1022.08453069756"
# [1] "objective:-1022.08341717686"
# [1] "objective:-1022.08267674839"

(Although the recommendation is to estimate the residual variance, to maintain consistency with the analysis using out-of-sample LD, we ignore this advice.)

SuSiE returns a single credible set (CS) containing a large number of strongly correlated SNPs. One of the SNPs in this CS is a (true) causal SNP.

print(fit1$sets[c("cs","purity")])
# $cs
# $cs$L1
#  [1] 195 197 203 213 226 237 238 243 247 248 249 254 255 278 294 296 301 305 319
# [20] 325 351 371 380 381 389 390 393 405 420 421 422 424 427 434 435 436 437 438
# [39] 441 442 443 445 448 450 452 454 456 459 462 464 466 467 468 473 477 478 479
# [58] 483 484 485 486 487 488 489 490 492 493 497 503 504 512 520 535 552 554 555
# [77] 558 571
# 
# 
# $purity
#    min.abs.corr mean.abs.corr median.abs.corr
# L1    0.9827454     0.9993761       0.9999995

Here's a visualization of this result. (In this plot, the CS is depicted by the light blue circles, and the two causal SNPs are drawn as red triangles.)

par(mar = c(4,4,1,1))
cs1 <- fit1$sets$cs$L1
plot(1:1001,fit1$pip,pch = 20,cex = 0.8,ylim = c(0,0.1),
     xlab = "SNP",ylab = "susie PIP")
points(cs1,fit1$pip[cs1],pch = 1,cex = 1,col = "cyan")
points(vars,fit1$pip[vars],pch = 2,cex = 0.8,col = "tomato")

Version Author Date
e597078 Peter Carbonetto 2022-06-10
f40a337 Peter Carbonetto 2022-06-09
37627a7 Peter Carbonetto 2022-06-09

Now let's try running FINEMAP on these same data:

run_finemap <- function (bhat, shat, maf, prefix) {
  p <- length(b)
  dat <- data.frame(rsid       = 1:p,
                    chromosome = rep(1,p),
                    position   = rep(1,p),
                    allele1    = rep("A",p),
                    allele2    = rep("C",p),
                    maf        = round(maf,digits = 6),
                    beta       = round(bhat,digits = 6),
                    se         = round(shat,digits = 6))
  outfile <- paste(prefix,"z",sep = ".")
  masterfile <- paste(prefix,"master",sep = ".")
  write.table(dat,outfile,quote = FALSE,col.names = TRUE,row.names = FALSE)
  out <- system(paste("./finemap_v1.4.1_x86_64 --sss --log --n-causal-snps 5",
                      "--in-files",masterfile),intern = TRUE)
  out <- out[which(!grepl("Reading",out,fixed = TRUE))]
  out <- out[which(!grepl("Computing",out,fixed = TRUE))]
  out <- out[which(!grepl("evaluated",out,fixed = TRUE))]
  cat(out,sep = "\n")
  return(invisible(out))
}
setwd("../output")
run_finemap(bhat,shat,maf,prefix = "sim1")
setwd("../analysis")
# 
# |--------------------------------------|
# | Welcome to FINEMAP v1.4.1            |
# |                                      |
# | (c) 2015-2022 University of Helsinki |
# |                                      |
# | Help :                               |
# | - ./finemap --help                   |
# | - www.finemap.me                     |
# | - www.christianbenner.com            |
# |                                      |
# | Contact :                            |
# | - finemap@christianbenner.com        |
# | - matti.pirinen@helsinki.fi          |
# |--------------------------------------|
# 
# --------
# SETTINGS
# --------
# - dataset            : all
# - corr-config        : 0.95
# - n-causal-snps      : 5
# - n-configs-top      : 50000
# - n-conv-sss         : 100
# - n-iter             : 100000
# - n-threads          : 1
# - prior-k0           : 0
# - prior-std          : 0.05 
# - prob-conv-sss-tol  : 0.001
# - prob-cred-set      : 0.95
# 
# ------------
# FINE-MAPPING (1/1)
# ------------
# - GWAS summary stats               : sim1.z
# - SNP correlations                 : ../data/small_data_11_sim_gaussian_pve_n_8_get_sumstats_n_1.ld_sample_n_file.in_n.ld
# - Causal SNP stats                 : sim1.snp
# - Causal configurations            : sim1.config
# - Credible sets                    : sim1.cred
# - Log file                         : sim1.log_sss
# 
# 
- Updated prior SD of effect sizes : 0.05 0.0902 0.163 0.293 
# 
# - Number of GWAS samples           : 800
# - Number of SNPs                   : 1001
# - Prior-Pr(# of causal SNPs is k)  : 
#   (0 -> 0)
#    1 -> 0.583
#    2 -> 0.291
#    3 -> 0.097
#    4 -> 0.0242
#    5 -> 0.00483
# - Regional SNP heritability        : 0.246 (SD: 0.0278 ; 95% CI: [0.194,0.303])
# - Log10-BF of >= one causal SNP    : 57.1
# - Post-expected # of causal SNPs   : 1.74
# - Post-Pr(# of causal SNPs is k)   : 
#   (0 -> 0)
#    1 -> 0.262
#    2 -> 0.738
#    3 -> 0
#    4 -> 0
#    5 -> 0
# - Run time                         : 0 hours, 0 minutes, 10 seconds

FINEMAP predicts 2 causal SNPs (with 74% probability). The second FINEMAP credible set is a diffuse set with lots of uncorrelated SNPs, so we visualize the first CS only:

par(mar = c(4,4,1,1))
finemap <- read.table("../output/sim1.cred2",header = TRUE)
pip   <- rep(0,1001)
cs1   <- finemap$cred1
rows1 <- which(!is.na(cs1))
cs1   <- cs1[rows1]
pip[cs1] <- pip[cs1] + finemap$prob1[rows1]
plot(1:1001,pip,pch = 20,cex = 0.8,ylim = c(0,0.1),
     xlab = "SNP",ylab = "finemap PIP")
points(cs1,pip[cs1],pch = 1,cex = 1,col = "cyan")
points(vars,pip[vars],pch = 2,cex = 0.8,col = "tomato")

Version Author Date
e597078 Peter Carbonetto 2022-06-10

This looks very similar to the SuSiE credible set, which is reassuring.

Out-of-sample LD

Let's try running SuSiE and FINEMAP as before, except using an out-of-sample LD estimate.

ldoutfile <- "../data/small_data_11.ld_refout_file.refout.ld"
Rout <- as.matrix(fread(ldoutfile))
fit2 <- susie_rss(z,Rout,n = 800,min_abs_corr = 0.1,refine = FALSE,
                 verbose = TRUE)
# HINT: If the in-sample LD matrix is available, we recommend calling susie_rss with the in-sample LD matrix, and setting estimate_residual_variance = TRUE
# HINT: For large R or large XtX, consider installing theRfast package for better performance.
# WARNING: XtX is not symmetric; forcing XtX to be symmetric by replacing XtX with (XtX + t(XtX))/2
# [1] "objective:-1022.15334027678"
# [1] "objective:-1022.1533319747"

As above, SuSiE returns a single CS containing one of the two causal SNPs.

print(fit2$sets[c("cs","purity")])
# $cs
# $cs$L1
#  [1] 195 197 203 213 226 237 238 243 247 248 249 254 255 278 294 296 301 305 319
# [20] 325 351 371 380 381 389 390 393 405 420 421 422 424 427 434 435 436 437 438
# [39] 441 442 443 445 448 450 452 454 456 459 462 464 466 467 468 473 477 478 479
# [58] 483 484 485 486 487 488 489 490 492 493 497 503 504 512 520 535 552 554 555
# [77] 558 571
# 
# 
# $purity
#    min.abs.corr mean.abs.corr median.abs.corr
# L1    0.9759149     0.9986827        0.999996

In fact, in this particular example the CS is exactly the same as the CS above obtained using the in-sample LD:

all(fit1$sets$cs$L1 == fit2$sets$cs$L1)
# [1] TRUE

Now let's run FINEMAP:

setwd("../output")
run_finemap(bhat,shat,maf,prefix = "sim2")
setwd("../analysis")
# 
# |--------------------------------------|
# | Welcome to FINEMAP v1.4.1            |
# |                                      |
# | (c) 2015-2022 University of Helsinki |
# |                                      |
# | Help :                               |
# | - ./finemap --help                   |
# | - www.finemap.me                     |
# | - www.christianbenner.com            |
# |                                      |
# | Contact :                            |
# | - finemap@christianbenner.com        |
# | - matti.pirinen@helsinki.fi          |
# |--------------------------------------|
# 
# --------
# SETTINGS
# --------
# - dataset            : all
# - corr-config        : 0.95
# - n-causal-snps      : 5
# - n-configs-top      : 50000
# - n-conv-sss         : 100
# - n-iter             : 100000
# - n-threads          : 1
# - prior-k0           : 0
# - prior-std          : 0.05 
# - prob-conv-sss-tol  : 0.001
# - prob-cred-set      : 0.95
# 
# ------------
# FINE-MAPPING (1/1)
# ------------
# - GWAS summary stats               : sim2.z
# - SNP correlations                 : ../data/small_data_11.ld_refout_file.refout.ld
# - Causal SNP stats                 : sim2.snp
# - Causal configurations            : sim2.config
# - Credible sets                    : sim2.cred
# - Log file                         : sim2.log_sss
# 
# 
- Updated prior SD of effect sizes : 0.05 0.0906 0.164 0.297 
# 
# - Number of GWAS samples           : 800
# - Number of SNPs                   : 1001
# - Prior-Pr(# of causal SNPs is k)  : 
#   (0 -> 0)
#    1 -> 0.583
#    2 -> 0.291
#    3 -> 0.097
#    4 -> 0.0242
#    5 -> 0.00483
# - Regional SNP heritability        : 0.986 (SD: 0.000581 ; 95% CI: [0.984,0.987])
# - Log10-BF of >= one causal SNP    : 185
# - Post-expected # of causal SNPs   : 5
# - Post-Pr(# of causal SNPs is k)   : 
#   (0 -> 0)
#    1 -> 0
#    2 -> 0
#    3 -> 4.93e-126
#    4 -> 4.06e-113
#    5 -> 1
# - Run time                         : 0 hours, 0 minutes, 23 seconds

With the out-of-sample LD estimate, FINEMAP gives a very different result: it predicts 5 causal SNPs (which is the maximum we allowed), all the CSs contain a single SNP (except the last CS, which has 2 SNPs), and none of the CSs contain a true causal SNP.

cat(readLines("../output/sim2.cred5"),sep = "\n")
# # Post-Pr(# of causal SNPs is 5) = 1
# #log10bf 600.715 NA 531.267 NA 135.593 NA 133.711 NA 117.309 NA
# #min(|ld|) 1 NA 1 NA 1 NA 1 NA 1 NA
# #mean(|ld|) 1 NA 1 NA 1 NA 1 NA 1 NA
# #median(|ld|) 1 NA 1 NA 1 NA 1 NA 1 NA
# index cred1 prob1 cred2 prob2 cred3 prob3 cred4 prob4 cred5 prob5
# 1 765 1 949 1 592 1 583 1 389 0.5
# 2 NA NA NA NA NA NA NA NA 405 0.5

sessionInfo()
# R version 3.5.1 (2018-07-02)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Scientific Linux 7.4 (Nitrogen)
# 
# Matrix products: default
# BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so
# 
# locale:
#  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
# [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] susieR_0.12.07    data.table_1.12.0
# 
# loaded via a namespace (and not attached):
#  [1] tidyselect_1.1.0   xfun_0.31          bslib_0.3.1        purrr_0.3.4       
#  [5] lattice_0.20-38    colorspace_1.3-2   vctrs_0.3.6        generics_0.0.2    
#  [9] htmltools_0.5.2    yaml_2.2.0         utf8_1.1.4         rlang_0.4.10      
# [13] mixsqp_0.3-46      jquerylib_0.1.4    later_0.7.5        pillar_1.5.0      
# [17] glue_1.4.2         DBI_1.0.0          matrixStats_0.54.0 lifecycle_1.0.0   
# [21] plyr_1.8.4         stringr_1.3.1      munsell_0.5.0      gtable_0.2.0      
# [25] workflowr_1.7.0    evaluate_0.15      knitr_1.39         fastmap_1.1.0     
# [29] httpuv_1.4.5       irlba_2.3.3        fansi_0.4.0        highr_0.7         
# [33] Rcpp_1.0.7         promises_1.0.1     backports_1.1.2    scales_1.0.0      
# [37] jsonlite_1.6       fs_1.5.0           ggplot2_3.3.3      digest_0.6.18     
# [41] stringi_1.2.4      dplyr_1.0.5        rprojroot_1.3-2    grid_3.5.1        
# [45] tools_3.5.1        magrittr_1.5       sass_0.4.1         tibble_3.1.0      
# [49] crayon_1.3.4       whisker_0.3-2      pkgconfig_2.0.2    ellipsis_0.3.2    
# [53] Matrix_1.2-15      rmarkdown_2.14     reshape_0.8.8      R6_2.3.0          
# [57] git2r_0.26.1       compiler_3.5.1