Version: 1.0.0
Authors: Dovini Jayasinghe and S. Hong Lee
The iPRSue is a novel method for estimating polygenic risk scores (PRSs) and their uncertainty (variance) using Genome Wide Association Studies (GWAS) summary statistics which can be constructed by GWAS_QT() or GWAS_BT() functions (see the section $\color{red} {IMPORTANT}$ for details). This method can be applied to both quantitative and binary traits through the functions iPRSue_estimates_QT() and iPRSue_estimates_BT(), providing unbiased and precise estimates. In addition, we have incorporated traditionally used Best linear unbiased estimator (BLUE) method for PRS and variance estimation, available for both trait types via BLUE_estimates_QT() and BLUE_estimates_BT(). While the BLUE method offers slightly but significantly higher prediction accuracy compared to iPRSue, it is limited to cases where the number of SNP genotypes is smaller than the discovery sample size. In contrast, iPRSue does not have this limitation and can be applied to any data dimension.
install.packages("iPRSue")
Install devtools:
install.packages("devtools")
Install iPRSue:
devtools::install_github("DoviniJ/iPRSue")
library(iPRSue)
The package iPRSue depends on the following R packages. Please install them and call the respective libraries.
install.packages(c("data.table", "bigstatsr", "logistf"))
library(data.table)
library(bigstatsr)
library(logistf)
ID_1 ID_1 0.967866365493657
ID_2 ID_2 -0.401199467463008
ID_3 ID_3 -0.1261351992308
ID_4 ID_4 -0.395945731035134
ID_5 ID_5 -0.131166319708339
ID_801 ID_801 -0.111301037986145
ID_802 ID_802 -1.43931139178665
ID_803 ID_803 0.133752876818589
ID_804 ID_804 -1.75782675330745
ID_805 ID_805 -0.799565331602807
ID_1 ID_1 1
ID_2 ID_2 2
ID_3 ID_3 1
ID_4 ID_4 1
ID_5 ID_5 1
ID_1 ID_1 0
ID_2 ID_2 1
ID_3 ID_3 0
ID_4 ID_4 0
ID_5 ID_5 0
ID_801 ID_801 1
ID_802 ID_802 2
ID_803 ID_803 1
ID_804 ID_804 1
ID_805 ID_805 2
ID_1 ID_1 -0.715681209642834 0.888725964414309 -1.03170535509313 1.13144382806602 0.974586116983568 1.44600211115777 0.231213420238433 -1.15908518979155 -0.260237749851626 -0.981928033025326 0.232145395698731 -0.00493421851594367 -1.02467943705236 -1.38309046960338
ID_2 ID_2 0.833805005322291 1.14151912762549 0.902862549695279 -1.09287371734997 0.382828629498391 -0.551686839383708 -0.285874699239703 1.87146117182544 0.0618826065454539 1.47080223642316 -0.487027551176616 -0.522925465187803 0.974695074269322 1.10982822777533
ID_3 ID_3 0.536893785388661 -0.248843270036006 1.568876367973 -0.396628752413828 -0.370603810588219 1.84225473719443 1.78102440765412 1.97939010218386 -0.790320886262911 0.593568016724195 0.678954789531668 -2.07574875557482 0.974695074269322 -0.807801539439061
ID_4 ID_4 2.95462492768725 -1.26001592288073 1.98375887378537 -1.85436124493203 -0.610043439191099 -0.0590784765825099 -0.358218902818925 1.47235715582949 -0.21521644298211 0.56767247164066 0.0320535549435434 0.617574422540434 0.974695074269322 1.10982822777533
ID_5 ID_5 0.0571302091333667 0.256743056386356 -0.364665321378934 -1.63324675483113 -0.185694458411861 0.0654274876199731 -0.903197651489459 -1.5774423323997 -1.72217861845979 -0.725173364377149 1.41905617836277 0.825306062484362 0.974695074269322 1.10982822777533
Note: If you have no confounder, then leave the 3rd column constant (e.g. vector of 1s)
cd1.txt
ID_1 ID_1 1
ID_2 ID_2 1
ID_3 ID_3 1
ID_4 ID_4 1
ID_5 ID_5 1
Gd.txt - This is a .txt file which contains scaled (column standardized) genotype matrix with the dimension discovery_sample_size x number_of_snps (e.g. 800 x 1,000), of the discovery individuals. Note that the file has neither row nor column headings. (Note: Gd_reduced.txt file contains discovery genotype information for 5 SNPs and is generated to support BLUE estimates where number of individuals cannot be less than number of SNPs)
Gt.txt - This is a .txt file which contains scaled (column standardized) genotype matrix with the dimension target_sample_size x number_of_snps (e.g. 200 x 1,000), of the target individuals. Note that the file has neither row nor column headings. (Note: Gt_reduced.txt file contains target genotype information for 5 SNPs and is generated to support BLUE estimates where number of individuals cannot be less than number of SNPs)
mydata.fam - This is a file associated with the PLINK binary format file which contains the following columns in order. The example dataset has 1,000 individuals. Note that the file has no column headings. This follows the PLINK .fam file format.
ID_1 ID_1 0 0 1 -9
ID_2 ID_2 0 0 2 -9
ID_3 ID_3 0 0 2 -9
ID_4 ID_4 0 0 2 -9
ID_5 ID_5 0 0 1 -9
1 SNP_1 0 768448 A G
1 SNP_2 0 853954 C A
1 SNP_3 0 880390 A C
1 SNP_4 0 940203 A G
1 SNP_5 0 987670 T G
Download all the above example data files from inst folder to your working directory. You can copy the path to inst folder https://github.com/DoviniJ/iPRSue/tree/main/inst and download the directory using https://download-directory.github.io/ page, for instant downloading. Follow the commands below to obtain the expected outputs. Once practiced, you may use your own data files to generate individual PRS confidence intervals by using either the novel method iPRSue or the BLUE method.
Link: https://www.cog-genomics.org/plink/2.0/
Give the path where plink executable file is located
plink_path <- "<plink_path>/plink2"
Note that, all these files (user's own data files) can be placed in a separate location. It is always upto the user's choice. In that case remember to give the full path to the file location since R identifies files by name, only when they are in the same directory. For example, if your data files are in a different directory;
"mydata" use "<path>/mydata"
Commands
x <- GWAS_QT(plink_path = "plink2", b_file = "mydata", discovery_pheno = "Qpd.txt", discovery_cov = "cd.txt", thread = 48)
y <- iPRSue_estimates_QT(gwas = x, target_pheno = "Qpt.txt", target_geno_mat = "Gt.txt", no_of_PRSs = 500, significance_level = 0.05, seed = set.seed(1))
In GWAS_QT(), plink_path is the path or executable name for PLINK2 (e.g., “plink2”); b_file is the prefix of the PLINK binary set (.bed/.bim/.fam) that defines the SNPs and discovery sample IDs; discovery_pheno provides the scaled quantitative phenotype file for the discovery sample, and discovery_cov supplies the matching covariate file (for example, sex and principal components), with IDs consistent with those in b_file; thread specifies the number of CPU threads to use (optional, default 20). The function runs linear-regression GWAS and returns additive SNP effect estimates and their standard errors, ordered exactly as in the corresponding .bim file.
In iPRSue_estimates_QT(), gwas is the GWAS results object produced by GWAS_QT(). The target_pheno argument points to the target sample’s quantitative phenotype file, while target_geno_mat is a text file containing the scaled target genotype matrix aligned to the GWAS SNP order. The argument no_of_PRSs sets the number of PRS draws/replicates used internally for uncertainty estimation, significance_level is the alpha level for confidence intervals (e.g., 0.05 for 95% CIs), and seed fixes the random seed for reproducibility. Using these inputs, the function computes iPRSue point estimates together with uncertainty metrics for each target individual.
Outputs
beta se
1 0.0182285503 0.03568526
2 -0.0517348657 0.03601479
3 0.0067835942 0.03327651
4 -0.0659989174 0.03386520
5 -0.0004306981 0.03517784
6 0.0284129553 0.03487432
Returns the dataframe x with following columns:
IID PRS Variance Lower_Limit Upper_Limit
1 ID_801 -0.8710867 1.326542 -3.1149828 1.2669814
2 ID_802 -2.9520295 1.194713 -5.2039386 -0.9847323
3 ID_803 0.2631768 1.382384 -2.1256907 2.5201000
4 ID_804 2.7146966 1.469916 0.2658072 5.1874619
5 ID_805 2.3592619 1.364946 0.1601901 4.6957551
6 ID_806 -1.9697418 1.210845 -4.0841447 0.2151955
Returns the dataframe y with following columns:
Command
z <- BLUE_estimates_QT(discovery_pheno = "Qpd.txt", discovery_geno_mat = "Gd_reduced.txt", target_pheno = "Qpt.txt", target_geno_mat = "Gt_reduced.txt", significance_level = 0.05)
In BLUE_estimates_QT(), discovery_pheno is the file containing scaled quantitative phenotypes for the discovery sample, and discovery_geno_mat is a text file with the corresponding scaled discovery genotype matrix. The target_pheno argument provides the target sample’s quantitative phenotypes, while target_geno_mat supplies the scaled target genotype matrix. The discovery and target genotype matrices must refer to the same SNP subset in the same column order, and sample IDs must align with their respective phenotype files. The significance_level sets the alpha for confidence intervals (e.g., 0.05 for 95% CIs). The function fits a BLUE multiple linear regression on the discovery data to estimate SNP effects, then applies those effects to the target genotypes to compute per-individual PRS along with uncertainty metrics (standard errors and confidence intervals) under the usual linear-model assumptions.
Output
IID PRS Variance Lower_Limit Upper_Limit
1 ID_801 0.08620601 0.002480396 -0.01140720 0.1838192
2 ID_802 0.03582957 0.029809734 -0.30256793 0.3742271
3 ID_803 0.06101779 0.007735024 -0.11135906 0.2333946
4 ID_804 -0.08975104 0.129695600 -0.79559826 0.6160962
5 ID_805 0.08620601 0.002480396 -0.01140720 0.1838192
6 ID_806 0.04970188 0.004847807 -0.08676296 0.1861667
Returns the dataframe z with following columns:
Commands
x <- GWAS_BT(plink_path = "plink2", b_file = "mydata", discovery_pheno = "Bpd.txt", discovery_cov = "cd.txt", thread = 48)
y <- iPRSue_estimates_BT(gwas = x, target_pheno = "Bpt.txt", target_geno_mat = "Gt.txt", no_of_PRSs = 500, significance_level = 0.05, seed = set.seed(1))
The arguments of GWAS_BT() function, namely, discovery_pheno and discovery_geno_mat should be the text file names of phenotype and scaled genotype matrix for the individuals in the discovery dataset. The function iPRSue_estimates() use the estimated values from GWAS_BT() function as an input to the argument gwas. Moreover, target_pheno and target_geno_mat should be the text file names of phenotype and scaled genotype matrix for the individuals in the target dataset. The no_of_PRSs specifies the size of the PRS distribution per target individual (default is 500), significance_level sets the level of significance for constructing PRS confidence intervals (default is 0.05), and seed is used to control random number generation for reproducibility (default is NULL).
Outputs
beta se
1 0.01128249 0.1828907
2 0.36766352 0.1820150
3 0.03322998 0.1599603
4 0.13476046 0.1486067
5 -0.08688650 0.1981246
6 -0.38928581 0.2070150
Returns the dataframe x with following columns:
IID PRS Variance Lower_Limit Upper_Limit
1 ID_801 4.146077 36.03383 -7.8230214 15.0955206
2 ID_802 -4.176091 35.71329 -16.4947165 6.9894897
3 ID_803 11.533248 39.03824 -0.3957665 23.5524960
4 ID_804 -3.503107 43.23646 -16.4752543 9.6227589
5 ID_805 -13.343637 40.25580 -25.1031448 -0.6833604
6 ID_806 6.039555 36.87239 -5.7316201 18.1581323
Returns the dataframe y with following columns:
Command
z <- BLUE_estimates_BT(discovery_pheno = "Bpd_0_1.txt", discovery_geno_mat = "Gd_reduced.txt", target_pheno = "Bpt.txt", target_geno_mat = "Gt_reduced.txt", significance_level = 0.05, max_iterations = 100)
In BLUE_estimates_BT(), discovery_pheno is the discovery sample’s binary phenotype file coded 0/1, and discovery_geno_mat is a text file containing the corresponding scaled discovery genotype matrix. The target_pheno argument supplies the target sample’s binary phenotypes, while target_geno_mat provides the scaled target genotype matrix. Discovery and target genotype matrices must reference the same SNP subset in identical column order, and sample IDs must align with their respective phenotype files. The significance_level sets the alpha for confidence intervals (e.g., 0.05 for 95% CIs), and max_iterations controls the cap on optimization steps for the Firth-corrected logistic fit. The function fits a BLUE multiple logistic regression with Firth’s bias correction on the discovery data to obtain SNP effects, then applies those effects to the target genotypes to compute per-individual PRS along with uncertainty metrics.
Output
IID PRS Variance Lower_Limit Upper_Limit
1 ID_801 -0.46466860 0.05424898 -0.9211719 -0.008165357
2 ID_802 -0.67096196 0.78880018 -2.4116927 1.069768757
3 ID_803 -0.56781528 0.19833426 -1.4406800 0.305049490
4 ID_804 1.06695650 2.35839577 -1.9429738 4.076886769
5 ID_805 -0.46466860 0.05424898 -0.9211719 -0.008165357
6 ID_806 0.04296879 0.10720391 -0.5987628 0.684700358
Returns the dataframe z with following columns:
When individual level data to use as a discovery sample are unavailable, users can often chose to utilize publicly available summary statistics data. However, these publicly available GWAS summary statistics often include only the (p)-value. In such cases, one can derive the Z-statistic ($Z_i$) assuming a standard normal distribution, and then compute the standard error as $\hat\sigma_i^{(\mathrm{scaled})} = 1/\sqrt{n}$ and effect size as $\hat\beta_i^{(\mathrm{scaled})} = Z_i/\sqrt{n}$, where $n$ corresponds to the discovery sample size. Once both are obtained, the iPRSue method can be applied using either iPRSue_estimates_QT() or iPRSue_estimates_BT() functions accordingly.