/prs-uncertainty

Primary LanguageR

Estimate PRS Uncertainty for Biobank

Here we demonstrate how to estimate PRS uncertainty (full posterior distribution of genetic value) as described in Large uncertainty in individual PRS estimation impacts PRS-based risk stratification.

Data

Prepare the following data:

  • train_bfile: Genotype file for training individuals.
  • train_sumstats: GWAS association summary statistics output from plink.
  • val_bfile: Genotype file for validation individuals.
  • val_pheno: Phenotype file for validation individuals.
  • test_bfile: Genotype file for testing individuals

Estimate PRS Uncertainty

Install the following packages:

install.packages('optparse', 'bigsnpr', 'tibble', 'readr', 'dplyr')

Download the prs_uncertainty.R and summary.R scripts:

Compute the PRS uncertainty:

#!/bin/sh
#$ -cwd
#$ -j y
#$ -l h_data=20G,h_rt=10:00:00 -pe shared 6
#$ -o ./job_out
#$ -t 1-22

chr_i=$SGE_TASK_ID

Rscript prs_uncertainty.R
    --chr_i=$chr_i
    --train_bfile=train_bfile \
    --train_sumstats=all.assoc.linear \
    --val_bfile=val_bfile \
    --val_pheno=height.val.regressed_pheno \
    --test_bfile=test_bfile \
    --out_dir=out_dir  #store outputfiles \
    --cache_dir=cache_file #store LD cache\ 
    --n_cores=6 \
    --num_burn_in=100 \
    --num_iter=500 \

To reduce runtime and memory required for training PRS models, we train PRS models separtely for each of the twentytwo chromosomes and then merge them to obtain a final PRS model for the whole genome using the script summary.R

#!/bin/sh
#$ -cwd
#$ -j y
#$ -l h_data=10G,h_rt=1:00:00 
#$ -o ./job_out

Rscript summary.R --out_dir=out_dir