Here I show how to compute polygenic scores using LDpred2, as well as inferring genetic architecture parameters with LDpred2-auto.

Here is a (slightly outdated) video of me going through the tutorial and explaining the different steps:

## Installation

In R, run

# install.packages("remotes")
remotes::install_github("privefl/bigsnpr")

or for the CRAN version

install.packages("bigsnpr")
• If you install {bigsnpr} >= v1.10.4, LDpred2-grid and LDpred2-auto should be much faster for large data.

• if you install {bigsnpr} >= v1.11.4, there is a new version LDpred2-auto that was validated for inferring parameters of genetic architectures (cf. at the end of this tutorial).

## Which set of variants to use?

To run LDpred2, you need

• GWAS summary statistics with marginal effect sizes, their standard errors, and the corresponding sample size(s),
• an LD (correlation) matrix computed from individuals of the same genetic ancestry as individuals used in the GWAS,
• individual-level data for tuning hyper-parameters (when using LDpred2-grid or lassosum2) and for testing the final models.

You need to restrict to genetic variants in common between all these datasets.

We previously recommended to use a set of 1,054,330 HapMap3 variants for LDpred2, because they provide a good coverage of the genome and are generally well imputed and available in most studies. In this preprint, we recently proposed an extension of this set, to provide an even better coverage of the genome by adding 37% more variants, and called it HapMap3+. This is the preferred set of variants to use with LDpred2 now.

You should use these sets of variants only when your data is imputed so that the overlap is good. If you only have access to a chip of genotyped variants, you should use these. In that case, you should compute the LD information yourself (as done for the tutorial data below).

If you use HM3/HM3+ variants with European summary statistics and do not have enough data to use as LD reference (e.g. at least 2000 individuals), we provide LD matrices to be used directly:

Note that forming independent LD blocks in LD matrices can be useful for robustness and extra speed gains (see this paper).

Information about the HapMap3+ variants can be retrieved with

# $pos is in build GRCh37 / hg19, but we provide positions in 2 other builds info <- readRDS(runonce::download_file( "https://figshare.com/ndownloader/files/37802721", dir = "tmp-data", fname = "map_hm3_plus.rds")) str(info) ## Classes 'tbl_df', 'tbl' and 'data.frame': 1444196 obs. of 10 variables: ##$ chr     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $pos : int 68082 77866 87409 87647 88144 89654 565935 568631 720797 741397 ... ##$ a0      : chr  "T" "C" "C" "T" ...
##  $a1 : chr "C" "T" "T" "C" ... ##$ rsid    : chr  "rs367789441" "rs563593912" "rs139490478" "rs146836579" ...
##  $af_UKBB : num 0.0708 0.07691 0.07683 0.01353 0.00857 ... ##$ ld      : num  1.81 1.3 1.57 1.34 1.16 ...
##  $block_id: int 1 1 1 1 1 1 1 1 1 1 ... ##$ pos_hg18: int  57945 67729 77272 77510 78007 79517 555798 558494 710660 731260 ...
##  $pos_hg38: int 68082 77866 87409 87647 88144 89654 630555 633251 785417 806017 ... ## Downloading genotype data and summary statistics for the tutorial This tutorial uses fake data for educational purposes only. Another tutorial using another dataset can be found at https://privefl.github.io/bigsnpr-extdoc/polygenic-scores-pgs.html. You can download the tutorial data and unzip files in R. We store those files in a directory called "tmp-data" here. # install.packages("runonce") zip <- runonce::download_file( "https://github.com/privefl/bigsnpr/raw/master/data-raw/public-data3.zip", dir = "tmp-data") unzip(zip) First, you need to read genotype data from the PLINK files (or BGEN files) as well as the text file containing summary statistics. # Load packages bigsnpr and bigstatsr library(bigsnpr) ## Loading required package: bigstatsr # Read from bed/bim/fam, it generates .bk and .rds files. snp_readBed("tmp-data/public-data3.bed") ## [1] "C:\\Users\\au639593\\Desktop\\bigsnpr\\tmp-data\\public-data3.rds" # Attach the "bigSNP" object in R session obj.bigSNP <- snp_attach("tmp-data/public-data3.rds") # See how the file looks like str(obj.bigSNP, max.level = 2, strict.width = "cut") ## List of 3 ##$ genotypes:Reference class 'FBM.code256' [package "bigstatsr"] with 16 ..
##   ..and 26 methods, of which 12 are  possibly relevant:
##   ..  add_columns, as.FBM, bm, bm.desc, check_dimensions,
##   ..  check_write_permissions, copy#envRefClass, initialize,
##   ..  initialize#FBM, save, show#envRefClass, show#FBM
##  $fam :'data.frame': 503 obs. of 6 variables: ## ..$ family.ID  : int [1:503] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$sample.ID : chr [1:503] "HG00096" "HG00097" "HG00099" "HG00100" ... ## ..$ paternal.ID: int [1:503] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$maternal.ID: int [1:503] 0 0 0 0 0 0 0 0 0 0 ... ## ..$ sex        : int [1:503] 1 2 2 2 1 2 1 1 2 1 ...
##   ..$affection : num [1:503] -0.547 0.188 -0.407 -0.28 0.398 ... ##$ map      :'data.frame':   45337 obs. of  6 variables:
##   ..$chromosome : int [1:45337] 1 1 1 1 1 1 1 1 1 1 ... ## ..$ marker.ID   : chr [1:45337] "rs3934834" "rs12726255" "rs11260549" "..
##   ..$genetic.dist: num [1:45337] 0.359 0.408 0.932 0.986 1.107 ... ## ..$ physical.pos: int [1:45337] 1005806 1049950 1121794 1162435 1314015..
##   ..$allele1 : chr [1:45337] "T" "G" "A" "A" ... ## ..$ allele2     : chr [1:45337] "C" "A" "G" "C" ...
##  - attr(*, "class")= chr "bigSNP"
# Get aliases for useful slots
G   <- obj.bigSNP$genotypes CHR <- obj.bigSNP$map$chromosome POS <- obj.bigSNP$map$physical.pos y <- obj.bigSNP$fam$affection (NCORES <- nb_cores()) ## [1] 4 # Read external summary statistics sumstats <- bigreadr::fread2("tmp-data/public-data3-sumstats.txt") str(sumstats) ## 'data.frame': 50000 obs. of 9 variables: ##$ rsid   : chr  "rs3934834" "rs12726255" "rs11260549" "rs3766186" ...
##  $chr : int 1 1 1 1 1 1 1 1 1 1 ... ##$ pos    : int  995669 1039813 1111657 1152298 1303878 1495118 1833906 2041373 2130121 2201709 ...
##  $a0 : chr "G" "T" "C" "C" ... ##$ a1     : chr  "A" "C" "T" "A" ...
##  $beta : num 0.0125 0.027 0.0171 -0.0195 -0.0057 ... ##$ beta_se: num  0.0157 0.0167 0.0179 0.0199 0.0213 ...
##  $N : int 15155 15155 15155 15155 15155 15155 15155 15155 15155 15155 ... ##$ p      : num  0.426 0.107 0.342 0.328 0.789 ...

We split the individuals from the genotype data into “validation” (to choose best-performing hyper-parameters) and “test” (to evaluate final polygenic scores). Here we consider that there are 350 individuals to be used as validation set to tune hyper-parameters for LDpred2-grid and lassosum2. The other 153 individuals are used as test set to evaluate the final models.

set.seed(1)
ind.val <- sample(nrow(G), 350)
ind.test <- setdiff(rows_along(G), ind.val)

## Matching variants between genotype data and summary statistics

To match variants contained in genotype data and summary statistics, the variables "chr" (chromosome number), "pos" (physical genetic position in bp), "a0" (reference allele) and "a1" (alternative allele) should be available in the summary statistics and in the genotype data. These 4 variables are used to match variants between the two data frames. From the summary statistics, you need to get "beta", "beta_se" (standard errors), and "n_eff" (the effective sample sizes per variant for a GWAS using logistic regression, and simply the sample size for continuous traits).

# sumstats$n_eff <- 4 / (1 / sumstats$n_case + 1 / sumstats$n_control) # sumstats$n_case <- sumstats$n_control <- NULL sumstats$n_eff <- sumstats$N map <- setNames(obj.bigSNP$map[-3], c("chr", "rsid", "pos", "a1", "a0"))
df_beta <- snp_match(sumstats, map)
## 50,000 variants to be matched.
## 0 ambiguous SNPs have been removed.
## 4 variants have been matched; 2 were flipped and 2 were reversed.
## Error: Not enough variants have been matched.

Here, there is problem with the matching; this is due to having different genome builds. You can either convert between builds with snp_modifyBuild() (or directly use the converted positions in info), or match by rsIDs instead.

df_beta <- snp_match(sumstats, map, join_by_pos = FALSE)  # use rsid instead of pos
## 50,000 variants to be matched.
## 0 ambiguous SNPs have been removed.
## 45,337 variants have been matched; 22,758 were flipped and 15,092 were reversed.

If no or few variants are actually flipped, you might want to disable the strand flipping option (strand_flip = FALSE) and maybe remove the few that were flipped (errors?).

## Quality control of GWAS summary statistics

Some quality control on summary statistics is highly recommended. A refined QC is described in this new paper. See e.g. the code that was used to prepare the sumstats there.

Most of the QC recommended for LDpred2 consists in computing $$$\text{sd}(G_j) \approx \dfrac{\text{sd}(y)}{\sqrt{n_j ~ \text{se}(\hat{\gamma}_j)^2 + \hat{\gamma}_j^2}} ~,$$$ where $$\hat{\gamma}_j$$ is the marginal GWAS effect size of variant $$j$$, $$n_j$$ is the GWAS sample size associated with variant $$j$$, $$y$$ is the vector of phenotypes and $${G}_j$$ is the vector of genotypes for variant $$j$$.

Note that, for a binary trait for which logistic regression is used, we have instead $$$\label{eq:approx-sd-log} \text{sd}(G_j) \approx \dfrac{2}{\sqrt{n_j^\text{eff} ~ \text{se}(\hat{\gamma}_j)^2 + \hat{\gamma}_j^2}} ~,$$$ where $$n_j^\text{eff} = \frac{4}{1 / n_j^\text{cases} + 1 / n_j^\text{controls}}$$.

Then, these standard deviations of genotypes inferred from GWAS summary statistics can be compared to $$\sqrt{2 \cdot f_j \cdot (1 - f_j) \cdot \text{INFO}_j}$$, where $$f_j$$ is the allele frequency of variant $$j$$ and $$\text{INFO}_j$$ is used to correct for the reduced variance of imputed dosages.

In the equation for continuous traits, we can estimate $$\text{sd}(y)$$ by the first percentile of $$\sqrt{0.5 \left(n_j ~ \text{se}(\hat{\gamma}_j)^2 + \hat{\gamma}_j^2\right)}$$, where the first percentile approximates the minimum and is robust to outliers. In practice, a consistent deviation between these two estimates of standard deviations (from summary statistics and from allele frequencies) can also be explained by using wrong estimates for $$n_j$$ or $$n_j^\text{eff}$$.

## Computing LDpred2 scores genome-wide

Note that you should run LDpred2 genome-wide. Just build the SFBM (the sparse LD matrix on disk) so that it contains selected variants for all chromosomes at once (see the for-loop below).

### Correlation

First, you need to compute correlations between variants. We recommend to use a window size of 3 cM (see the LDpred2 paper).

# To convert physical positions (in bp) to genetic positions (in cM), use
# POS2 <- snp_asGeneticPos(CHR, POS, dir = "tmp-data", ncores = NCORES)
# To avoid downloading "large" files, POS2 has been precomputed here
POS2 <- obj.bigSNP$map$genetic.dist

Let us create the on-disk sparse genome-wide correlation matrix on-the-fly:

tmp <- tempfile(tmpdir = "tmp-data")

for (chr in 1:22) {

# print(chr)

## indices in 'df_beta'
ind.chr <- which(df_beta$chr == chr) ## indices in 'G' ind.chr2 <- df_beta$_NUM_ID_[ind.chr]

corr0 <- snp_cor(G, ind.col = ind.chr2, size = 3 / 1000,
infos.pos = POS2[ind.chr2], ncores = NCORES)

if (chr == 1) {
ld <- Matrix::colSums(corr0^2)
corr <- as_SFBM(corr0, tmp, compact = TRUE)
} else {
ld <- c(ld, Matrix::colSums(corr0^2))
corr$add_columns(corr0, nrow(corr)) } } To use the “compact” format for SFBMs, you need packageVersion("bigsparser") >= package_version("0.5"). Make sure to reinstall {bigsnpr} after updating {bigsparser} to this new version (to avoid crashes). file.size(corr$sbk) / 1024^3  # file size in GB
## [1] 0.03205058

Note that you will need at least the same memory as this file size (to keep it cached for faster processing) + some other memory for all the results returned. If you do not have enough memory, processing will be very slow (because you would read the data from disk all the time). If using the one million HapMap3 variants, having 60 GB of memory should be enough.

### LDpred2-inf: infinitesimal model

# Estimate of h2 from LD Score regression
(ldsc <- with(df_beta, snp_ldsc(ld, length(ld), chi2 = (beta / beta_se)^2,
sample_size = n_eff, blocks = NULL)))
##       int        h2
## 0.9496371 0.2981213
h2_est <- ldsc[["h2"]]
beta_inf <- snp_ldpred2_inf(corr, df_beta, h2 = h2_est)
pred_inf <- big_prodVec(G, beta_inf, ind.row = ind.test, ind.col = df_beta[["_NUM_ID_"]])
pcor(pred_inf, y[ind.test], NULL)
## [1] 0.3171149 0.1668316 0.4529985

LDpred2-inf would very likely perform worse than the other models presented hereinafter. We actually recommend not to use it anymore.

### LDpred2(-grid): grid of models

In practice, we recommend to test multiple values for h2 and p.

(h2_seq <- round(h2_est * c(0.3, 0.7, 1, 1.4), 4))
## [1] 0.0894 0.2087 0.2981 0.4174
(p_seq <- signif(seq_log(1e-5, 1, length.out = 21), 2))
##  [1] 1.0e-05 1.8e-05 3.2e-05 5.6e-05 1.0e-04 1.8e-04 3.2e-04 5.6e-04
##  [9] 1.0e-03 1.8e-03 3.2e-03 5.6e-03 1.0e-02 1.8e-02 3.2e-02 5.6e-02
## [17] 1.0e-01 1.8e-01 3.2e-01 5.6e-01 1.0e+00
(params <- expand.grid(p = p_seq, h2 = h2_seq, sparse = c(FALSE, TRUE)))
##          p     h2 sparse
## 1  1.0e-05 0.0894  FALSE
## 2  1.8e-05 0.0894  FALSE
## 3  3.2e-05 0.0894  FALSE
## 4  5.6e-05 0.0894  FALSE
## 5  1.0e-04 0.0894  FALSE
## 6  1.8e-04 0.0894  FALSE
## 7  3.2e-04 0.0894  FALSE
## 8  5.6e-04 0.0894  FALSE
## 9  1.0e-03 0.0894  FALSE
## 10 1.8e-03 0.0894  FALSE
##  [ reached 'max' / getOption("max.print") -- omitted 158 rows ]
set.seed(1)  # to get the same result every time
# takes less than 2 min with 4 cores
beta_grid <- snp_ldpred2_grid(corr, df_beta, params, ncores = NCORES)
pred_grid <- big_prodMat(G, beta_grid, ind.col = df_beta[["_NUM_ID_"]])
params$score <- apply(pred_grid[ind.val, ], 2, function(x) { if (all(is.na(x))) return(NA) summary(lm(y[ind.val] ~ x))$coef["x", 3]
# summary(glm(y[ind.val] ~ x, family = "binomial"))$coef["x", 3] }) Note that missing values represent models that diverged substantially. library(ggplot2) ggplot(params, aes(x = p, y = score, color = as.factor(h2))) + theme_bigstatsr() + geom_point() + geom_line() + scale_x_log10(breaks = 10^(-5:0), minor_breaks = params$p) +
facet_wrap(~ sparse, labeller = label_both) +
labs(y = "GLM Z-Score", color = "h2") +
theme(legend.position = "top", panel.spacing = unit(1, "lines"))

library(dplyr)
params %>%
mutate(sparsity = colMeans(beta_grid == 0), id = row_number()) %>%
arrange(desc(score)) %>%
mutate_at(c("score", "sparsity"), round, digits = 3) %>%
slice(1:10)
##       p     h2 sparse  score sparsity  id
## 1 0.018 0.2981  FALSE 10.274    0.000  56
## 2 0.010 0.4174  FALSE 10.256    0.000  76
## 3 0.018 0.4174   TRUE 10.217    0.775 161
## 4 0.010 0.2981  FALSE 10.212    0.000  55
## 5 0.010 0.2087  FALSE 10.208    0.000  34
##  [ reached 'max' / getOption("max.print") -- omitted 5 rows ]

You can then choose the best model according to your preferred criterion (e.g. max AUC or $$r^2$$). Here, we use the Z-Score from the (linear or logistic) regression of the phenotype by the PRS since we have found it more robust than using the correlation or the AUC. It also enables adjusting for covariates in this step.

Also note that we separate both sparse and non-sparse models to show that their predictive performance are similar (in the original LDpred2 paper). In practice, if you do not really care about sparsity, you could choose the best LDpred2-grid model among all sparse and non-sparse models. If you do, choose the best sparse one (if it is close enough to the best one).

best_beta_grid <- params %>%
mutate(id = row_number()) %>%
# filter(sparse) %>%
arrange(desc(score)) %>%
slice(1) %>%
print() %>%
pull(id) %>%
beta_grid[, .]
##       p     h2 sparse    score id
## 1 0.018 0.2981  FALSE 10.27429 56
pred <- big_prodVec(G, best_beta_grid, ind.row = ind.test,
ind.col = df_beta[["_NUM_ID_"]])
pcor(pred, y[ind.test], NULL)
## [1] 0.4978420 0.3682548 0.6084543

### LDpred2-auto: automatic model

Contrary to LDpred2-grid, LDpred2-auto does not need a validation set because it can directly infer values for hyper-parameters $$h^2$$ and $$p$$. We recommend to run many chains in parallel with different initial values for p, which will be used for QC afterwards. In this paper, we have also introduced two new parameters in LDpred2-auto for improving its robustness, allow_jump_sign and shrink_corr, and recommend to use them.

set.seed(1)  # to get the same result every time
# takes less than 2 min with 4 cores
multi_auto <- snp_ldpred2_auto(
corr, df_beta, h2_init = h2_est,
vec_p_init = seq_log(1e-4, 0.2, length.out = 30), ncores = NCORES,
# use_MLE = FALSE,  # uncomment if you have convergence issues (need v1.11.9)
allow_jump_sign = FALSE, shrink_corr = 0.95)
str(multi_auto, max.level = 1)
## List of 30
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
##  $:List of 12 ##$ :List of 12
str(multi_auto[[1]], max.level = 1)
## List of 12
##  $beta_est : num [1:45337] 5.92e-05 2.94e-04 9.30e-05 -1.10e-04 2.28e-05 ... ##$ postp_est     : num [1:45337] 0.00564 0.01305 0.00661 0.00699 0.00503 ...
##  $corr_est : num [1:45337] 3.60e-05 1.87e-04 1.24e-05 -6.52e-05 -3.23e-05 ... ##$ sample_beta   :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##  $path_p_est : num [1:700] 0.000567 0.001407 0.001638 0.002563 0.003236 ... ##$ path_h2_est   : num [1:700] 0.084 0.113 0.132 0.127 0.143 ...
##  $path_alpha_est: num [1:700] 0.5 0.5 0.5 0.104 -0.328 ... ##$ h2_est        : num 0.209
##  $p_est : num 0.011 ##$ alpha_est     : num -0.523
##  $h2_init : num 0.298 ##$ p_init        : num 1e-04

You can verify whether the chains “converged” by looking at the path of the chains:

library(ggplot2)
auto <- multi_auto[[1]]  # first chain
plot_grid(
qplot(y = auto$path_p_est) + theme_bigstatsr() + geom_hline(yintercept = auto$p_est, col = "blue") +
scale_y_log10() +
labs(y = "p"),
qplot(y = auto$path_h2_est) + theme_bigstatsr() + geom_hline(yintercept = auto$h2_est, col = "blue") +
labs(y = "h2"),
ncol = 1, align = "hv"
)
## Warning: qplot() was deprecated in ggplot2 3.4.0.

In the LDpred2 paper, we proposed an automatic way of filtering bad chains by comparing the scale of the resulting predictions. We have tested a somewhat equivalent and simpler alternative since, which we recommend here:

# range should be between 0 and 2
(range <- sapply(multi_auto, function(auto) diff(range(auto$corr_est)))) ## [1] 0.1210443 0.1216797 0.1209822 0.1197119 0.1187931 0.1199233 0.1202114 ## [8] 0.1195135 0.1214119 0.1206811 0.1189464 0.1204088 0.1196642 0.1195328 ## [15] 0.1198751 0.1225441 0.1210127 0.1210234 0.1196245 0.1194824 0.1213715 ## [22] 0.1188433 0.1203081 0.1196867 0.1210735 0.1201303 0.1209825 0.1195834 ## [29] 0.1200190 0.1197259 (keep <- (range > (0.95 * quantile(range, 0.95)))) ## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [29] TRUE TRUE To get the final effects / predictions, you should only use chains that pass this filtering: beta_auto <- rowMeans(sapply(multi_auto[keep], function(auto) auto$beta_est))
pred_auto <- big_prodVec(G, beta_auto, ind.row = ind.test, ind.col = df_beta[["_NUM_ID_"]])
pcor(pred_auto, y[ind.test], NULL)
## [1] 0.4967147 0.3669594 0.6075102

### lassosum2: grid of models

lassosum2 is a re-implementation of the lassosum model that now uses the exact same input parameters as LDpred2 (corr and df_beta). It should be fast to run. It can be run next to LDpred2 and the best model can be chosen using the validation set. Note that parameter ‘s’ from lassosum has been replaced by a new parameter ‘delta’ in lassosum2, in order to better reflect that the lassosum model also uses L2-regularization (therefore, elastic-net regularization).

beta_lassosum2 <- snp_lassosum2(corr, df_beta, ncores = NCORES)
(params2 <- attr(beta_lassosum2, "grid_param"))
##       lambda delta num_iter time  sparsity
## 1 0.06111003 0.001        2 0.00 0.9998897
## 2 0.05241382 0.001        3 0.01 0.9997133
## 3 0.04495512 0.001        4 0.00 0.9995368
## 4 0.03855782 0.001        7 0.01 0.9992501
## 5 0.03307088 0.001        8 0.00 0.9984119
## 6 0.02836476 0.001       13 0.02 0.9971105
##  [ reached 'max' / getOption("max.print") -- omitted 114 rows ]
pred_grid2 <- big_prodMat(G, beta_lassosum2, ind.col = df_beta[["_NUM_ID_"]])
params2$score <- apply(pred_grid2[ind.val, ], 2, function(x) { if (all(is.na(x))) return(NA) summary(lm(y[ind.val] ~ x))$coef["x", 3]
# summary(glm(y[ind.val] ~ x, family = "binomial"))$coef["x", 3] }) library(ggplot2) ggplot(params2, aes(x = lambda, y = score, color = as.factor(delta))) + theme_bigstatsr() + geom_point() + geom_line() + scale_x_log10(breaks = 10^(-5:0)) + labs(y = "GLM Z-Score", color = "delta") + theme(legend.position = "top") + guides(colour = guide_legend(nrow = 1)) ## Warning: Removed 10 rows containing missing values (geom_point()). ## Warning: Removed 10 rows containing missing values (geom_line()). library(dplyr) best_grid_lassosum2 <- params2 %>% mutate(id = row_number()) %>% arrange(desc(score)) %>% print() %>% slice(1) %>% pull(id) %>% beta_lassosum2[, .] ## lambda delta num_iter time sparsity score id ## 1 0.02086632 0.001 21 0.01 0.9864570 9.715332 8 ## 2 0.02086632 0.010 20 0.01 0.9864790 9.709724 38 ## 3 0.02086632 0.100 16 0.01 0.9862364 9.667190 68 ## 4 0.02432834 0.001 19 0.00 0.9942652 9.404453 7 ## [ reached 'max' / getOption("max.print") -- omitted 116 rows ] # Choose the best among all LDpred2-grid and lassosum2 models best_grid_overall <- if(max(params2$score, na.rm = TRUE) > max(params$score, na.rm = TRUE), best_grid_lassosum2, best_beta_grid) ## Inference with LDpred2-auto Here, you need packageVersion("bigsnpr") >= package_version("1.11.4"). Here is the accompanying preprint. LDpred2-auto has been very recently extended and improved to allow for estimating $$h^2$$, $$p$$, and $$\alpha$$, a third parameter that controls how expected effect sizes relate to minor allele frequencies. The new model assumed by LDpred2-auto is $$$\beta_j = S_j \gamma_j \sim \left\{ \begin{array}{ll} \mathcal N \big( 0,~\sigma_\beta^2 \cdot (S_j^2)^{(\alpha + 1)} \big) & \mbox{with probability p,} \\ 0 & \mbox{otherwise,} \end{array} \right.$$$ where $$p$$ is the proportion of causal variants, $$h^2$$ the (SNP) heritability, $$\boldsymbol{\gamma}$$ the effect sizes on the allele scale, $$\boldsymbol{S}$$ the standard deviations of the genotypes, and $$\boldsymbol{\beta}$$ the effects of the scaled genotypes. Parameters $$h^2$$, $$p$$, and $$\alpha$$ (and 95% CIs) can be estimated using: coef_shrink <- 0.95 multi_auto <- snp_ldpred2_auto(corr, df_beta, h2_init = h2_est, vec_p_init = seq_log(1e-4, 0.2, length.out = 50), burn_in = 500, num_iter = 500, report_step = 20, allow_jump_sign = FALSE, shrink_corr = coef_shrink, ncores = NCORES) (range <- sapply(multi_auto, function(auto) diff(range(auto$corr_est))))
##  [1] 0.1203440 0.1201514 0.1201017 0.1199062 0.1200198 0.1217100 0.1192894
##  [8] 0.1206322 0.1194016 0.1205341 0.1209561 0.1209005 0.1202831 0.1203131
## [15] 0.1199622 0.1207713 0.1205029 0.1200927 0.1213750 0.1194640 0.1204270
## [22] 0.1200989 0.1205159 0.1199131 0.1200470 0.1212829 0.1207075 0.1203341
## [29] 0.1196193 0.1208523
##  [ reached getOption("max.print") -- omitted 20 entries ]
(keep <- which(range > (0.95 * quantile(range, 0.95, na.rm = TRUE)) & range < 3))
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30
##  [ reached getOption("max.print") -- omitted 20 entries ]
all_h2 <- sapply(multi_auto[keep], function(auto) tail(auto$path_h2_est, 500)) quantile(all_h2, c(0.5, 0.025, 0.975)) ## 50% 2.5% 97.5% ## 0.2178006 0.1894810 0.2483018 all_p <- sapply(multi_auto[keep], function(auto) tail(auto$path_p_est, 500))
quantile(all_p, c(0.5, 0.025, 0.975))
##         50%        2.5%       97.5%
## 0.012294107 0.008627975 0.017112622
all_alpha <- sapply(multi_auto[keep], function(auto) tail(auto$path_alpha_est, 500)) quantile(all_alpha, c(0.5, 0.025, 0.975)) ## 50% 2.5% 97.5% ## -0.4925350 -0.9062350 -0.0419629 Predictive performance $$r^2$$ can also be inferred from the Gibbs sampler: bsamp <- lapply(multi_auto[keep], function(auto) auto$sample_beta)
all_r2 <- do.call("cbind", lapply(seq_along(bsamp), function(ic) {
b1 <- bsamp[[ic]]
Rb1 <- apply(b1, 2, function(x)
coef_shrink * bigsparser::sp_prodVec(corr, x) + (1 - coef_shrink) * x)
b2 <- do.call("cbind", bsamp[-ic])
b2Rb1 <- as.matrix(Matrix::crossprod(b2, Rb1))
}))
quantile(all_r2, c(0.5, 0.025, 0.975))
##       50%      2.5%     97.5%
## 0.1150057 0.1044237 0.1264264

and compared to

beta_auto <- rowMeans(sapply(multi_auto[keep], function(auto) auto$beta_est)) pred_auto <- big_prodVec(G, beta_auto, ind.col = df_beta[["_NUM_ID_"]]) pcor(pred_auto, y, NULL)^2 ## [1] 0.2377858 0.1747458 0.3042043 These are not exactly the same, which we attribute to the small number of variants used in this tutorial data. You can also get per-variant probabilities of being causal (for fine-mapping purposes) postp <- rowMeans(sapply(multi_auto[keep], function(auto) auto$postp_est))
qplot(y = postp, alpha = I(0.2)) + theme_bigstatsr()

You can also infer local per-block or per-annotation $$h^2$$ following the code here and here.

# Some cleaning
rm(corr); gc(); file.remove(paste0(tmp, ".sbk"))
##            used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells  3922430 209.5    6973624 372.5  6973624 372.5
## Vcells 43354351 330.8   75127815 573.2 75127815 573.2
## [1] TRUE