---

Title: "GxEprs" Authors: "Dovini Jayasinghe, Md Moksedul Momin and Hong Lee"

GxEprs

The 'GxEprs' is an R package to detect and estimate GxE. It uses a novel PRS model that can enhance the prediction accuracy by utilising GxE effects. Firstly it performs Genome Wide Association Studies (GWAS) and Genome Wide Environment Interaction Studies (GWEIS) using the discovery dataset (see functions GWAS_binary(),GWAS_quantitative(), GWEIS_binary(), GWEIS_quantitative()). See the section $\color{red}{IMPORTANT}$ for the discovery models used. Secondly, it uses the GWAS and GWEIS summary statistics generated from the fucntions above to obtain polygenic risk scores (PRSs) (see functions PRS_binary() and PRS_quantitative()) for the target sample. Finally it predicts the risk values of each individual in the target sample (see functions summary_regular_binary() and summary_regular_quantitative()). Note that the users can fit 4 different models when the outcome is a quantitative trait, and 5 different models when the outcome is a binary disease trait. See the section $\color{red} {IMPORTANT}$ for the target models used. Finally, it is recommended to check the p-value from permutations using Model 4 (see function summary_permuted_quantitative()), and Model 5 (see functionsummary_permuted_binary()), to make sure that the significance of GxE is not spurious due to model misspecification (Jayasinghe et. al.).

Package installation

The current GitHub version of GxEprs can be installed via:

library(devtools)
install_github("DoviniJ/GxEprs") 

or the CRAN version via:

install.packages("GxEprs") 

Load the library

library(GxEprs)

Data preparation

File formats

Input files

1. 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.

• family ID (FID)
• individual ID (IID)
• father's ID
• mother's ID
• sex
• phenotype value

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

2. mydata.bim - This is is a file associated with the PLINK binary format file which contains the following columns in order. The example dataset has 1,000 SNPs. Note that the file has no column headings. This follows the PLINK .bim file format.

• chromosome code
• SNP ID
• position of centimorgans
• base-pair coordinate
• minor allele
• reference allele

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

3. mydata.bed - This is the PLINK binary format file which includes genotype information. This follows the PLINK .bed file format.
4. Bpd.txt - This is a .txt file which contains the following columns in order. The discovery dataset has 800 individuals. Note that the file has no column headings.

• FID
• IID
• binary phenotype (1=controls, 2=cases) of the discovery sample

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

5. Bcd.txt - This is a .txt file which contains the following columns in order. The discovery dataset has 800 individuals. Note that the file has no column headings.

• FID
• IID
• standardized covariate
• square of the standardized covariate (Note: This column is required to address model misspecification, if present. See reference for details.)
• 14 confounders of the discovery sample (Note: These columns are optional. Can use any number of columns as confounders to adjust the phenotype upon user requirement.)

ID_1 ID_1  0.8041577 -0.1393189 -0.44222329 -1.4193491  0.1347118 -1.0567922  0.36022401 -0.55028448 -0.7246779 -1.2961161 -1.25407218 -0.7776324 -0.79831383 -0.2090511   0  10
ID_2 ID_2 -0.1119921 -0.3742044 -1.16598893  0.9200682 -1.3906769  1.4219627 -1.40585866  0.81717581  0.9228321 -1.9086084  0.05270867  2.0209031  0.91342992 -1.0045477   1  19
ID_3 ID_3  0.1840184 -0.3681062 -0.04625744  0.3044321 -0.2966651  1.3790970  0.97412531  0.20573056  0.8691321  0.3116082 -0.69680184  0.3584836  0.20014714 -0.7928459   1   7
ID_4 ID_4 -0.6973568 -0.1901045 -0.01990925  1.4125771  1.1931358  0.5793260 -0.18750284 -0.57664383 -0.4667130  0.7389666 -0.01384177  1.2543045 -0.04840392  0.2529956   0  20
ID_5 ID_5  3.7386616  4.9092974 -0.48295819  1.5357043  2.2851890  0.6243784  0.02419407 -2.91672253 -0.7156356  1.3521021  1.95512167  0.7267245  0.16639161 -0.9327111   1  20
ID_6 ID_6  1.7694612  0.7986924 -0.85380248  0.9200682 -0.4687957 -0.6998884  0.73996797  0.03228818 -0.0421946 -0.7071904  0.48762762 -0.2809428  1.38168040 -0.7227888   1  15

6. Bpt.txt - This is a .txt file which contains the following columns in order. The target dataset has 200 individuals who are independent from the discovery dataset. Note that the file has no column headings.

• FID
• IID
• binary phenotype (0=controls, 1=cases) of the target sample

ID_801 ID_801 0
ID_802 ID_802 1
ID_803 ID_803 0
ID_804 ID_804 0
ID_805 ID_805 1

7. Bct.txt - This is a .txt file which contains the following columns in order. The target dataset has 200 individuals who are independent from the discovery dataset. Note that the file has no column headings.

• FID
• IID
• standardized covariate
• square of the standardized covariate (Note: This column is required to address model misspecification, if present. See reference for details.)
• 14 confounders of the target sample (Note: These columns are optional. Can use any number of columns as confounders to adjust the phenotype upon user requirement.)

ID_801 ID_801 -0.5131023 -0.3734833 -0.8635131  0.07458083 -0.7505894  1.13088081 -1.9116502 -0.18663504 -0.56556069 -0.005592953  1.4391434  0.363682416  0.1412291 -1.2498470   1   7
ID_802 ID_802 -0.1599420 -0.4582160  1.7090613 -0.04195172 -0.8165283 -0.46395143  0.5117234 -0.15821914 -0.09286106  0.056003549  2.4639712  1.184874461 -1.3725107  0.7221745   1  19
ID_803 ID_803 -1.4541690  0.4138515 -0.4870243  0.54071104  1.4915497 -0.17995300  0.4740690 -0.10666278 -0.59424960 -0.441967244 -1.1969526 -0.417668942  0.1143860  0.9186369   0  20
ID_804 ID_804  0.6481305 -0.2194557  2.3160535 -0.85767959 -0.1411417  0.17691493 -0.5480497  1.60549905  2.34863885 -0.160569475  1.3116324  0.006854464  1.0720956  0.2352379   0  10
ID_805 ID_805 -0.7592490 -0.2464120 -0.6506579 -0.04195172 -0.2966719 -0.68373479 -0.6326805  0.79193870  0.87519018 -1.530327179 -0.4491451  1.245483066  0.5153288 -1.8453277   0   7
ID_806 ID_806 -0.4762691 -0.3876925 -0.7696400 -1.32380980 -1.2491145  0.01209383 -0.4076871  0.07842574  0.51256343  1.061430666 -0.4944131  0.753688465  0.5268681  0.4881057   1  19

8. Qpd.txt - This is a .txt file which contains the following columns in order. The discovery dataset has 800 individuals. Note that the file has no column headings.

• FID
• IID
• quantitative phenotype of the discovery sample

ID_1 ID_1  0.9678664
ID_2 ID_2 -0.4011995
ID_3 ID_3 -0.1261352
ID_4 ID_4 -0.3959457
ID_5 ID_5 -0.1311663
ID_6 ID_6  2.5648687

9. Qcd.txt - This is a .txt file which contains the following columns in order. The discovery dataset has 800 individuals. Note that the file has no column headings.

• FID
• IID
• standardized covariate
• square of the standardized covariate (Note: Although squared covariate term is not involved when analyzing quantitative traits, this file can/may be used as the covariate file when modulating with a binary outcome. In that event, this column is required to address model misspecification, if present. See reference for details.)
• 14 confounders of the discovery sample

ID_1 ID_1 -0.6020207 -0.51317685 -0.71568121  0.8887260 -1.0317054  1.1314438  0.9745861  1.44600211  0.2312134 -1.159085 -0.26023775 -0.9819280  0.23214540 -0.004934219   0   7
ID_2 ID_2  0.0186192 -0.88475796  0.83380501  1.1415191  0.9028625 -1.0928737  0.3828286 -0.55168684 -0.2858747  1.871461  0.06188261  1.4708022 -0.48702755 -0.522925465   1  20
ID_3 ID_3  2.1908589  3.18159877  0.53689379 -0.2488433  1.5688764 -0.3966288 -0.3706038  1.84225474  1.7810244  1.979390 -0.79032089  0.5935680  0.67895479 -2.075748756   1  10
ID_4 ID_4  2.1908589  3.18159877  2.95462493 -1.2600159  1.9837589 -1.8543612 -0.6100434 -0.05907848 -0.3582189  1.472357 -0.21521644  0.5676725  0.03205355  0.617574423   1  20
ID_5 ID_5 -0.9123407 -0.07182008  0.05713021  0.2567431 -0.3646653 -1.6332468 -0.1856945  0.06542749 -0.9031977 -1.577442 -1.72217862 -0.7251734  1.41905618  0.825306062   1  20
ID_6 ID_6  0.0186192 -0.88475796 -0.91012230 -0.6280330  0.2949665 -1.8988156 -0.7582624 -0.04166171 -0.4340410 -0.155675 -1.11987653 -0.1821531 -0.43905130  0.479906553   1  10

10. Qpt.txt - This is a .txt file which contains the following columns in order. The target dataset has 200 individuals who are independent from the discovery dataset. Note that the file has no column headings.

• FID
• IID
• quantitative phenotype of the target sample

ID_801 ID_801 -0.1113010
ID_802 ID_802 -1.4393114
ID_803 ID_803  0.1337529
ID_804 ID_804 -1.7578268
ID_805 ID_805 -0.7995653
ID_806 ID_806 -0.4449718

11. Qct.txt - This is a .txt file which contains the following columns in order. The target dataset has 200 individuals who are independent from the discovery dataset. Note that the file has no column headings.

• FID
• IID
• standardized covariate
• square of the standardized covariate (Note: Although squared covariate term is not involved when analyzing quantitative traits, this file can/may be used as the covariate file when modulating with a binary outcome. In that event, this column is required to address model misspecification, if present. See reference for details.)
• 14 confounders of the target sample

ID_801 ID_801 -0.5248929 -0.44946058 -0.7134238 1.4813617 -0.99625567 0.1148177 -0.003464572 -0.388698341 -0.7729983  2.172350119  0.2971839 -0.4673828  0.7151812  0.5900938   0   7
ID_802 ID_802 -0.8560557  0.02018765  1.6066575 0.3952309  0.02677194 0.1674581  0.533868884  1.581875902  1.4101015 -0.569170793 -0.7407161  0.4997040  1.1072589 -0.1405227   1  19
ID_803 ID_803  0.1374325 -0.84486046  0.4800323 0.6365933  0.86951606 2.0404667 -0.693365818  0.009317369 -0.6600856 -0.004717947 -0.4101216 -0.2819847 -0.3857549  0.7299139   0   7
ID_804 ID_804 -0.5248929 -0.44946058 -0.2386328 1.2399993  1.50104915 0.2116517 -0.348460870 -0.171542519 -0.1042097 -0.815440139  0.3767518 -1.5353941 -1.9075240  1.1727743   0  10
ID_805 ID_805  0.4685953 -0.77061209 -0.5185197 1.2399993  0.26865598 1.3655378  1.738861438  0.505504084 -0.9220277 -1.182118202 -2.2407192 -0.5302259  0.4076562  0.5962252   1  15
ID_806 ID_806  1.1309208 -0.07821876  2.0740335 0.5159121  0.47693972 1.3072682 -0.322090873  1.092847799  0.3954790 -0.865146524 -0.7739106  0.5371777  0.8537627  0.3906600   1   7

Quick start

We will illustrate the usage of GxEprs using a few example datasets. Follow the step wise process:

Step 1: Download and install plink2 in your machine. The package supports both Linux and Windows version.

Link: https://www.cog-genomics.org/plink/2.0/

Step 2: Obtain the path to the executable plink application <plink_path>

Step 3.1: Run the following code (functions should be run in the given order) to obtain the risk scores of individuals in the target dataset:

Step 3.1.1 Give the path where plink executable file is located

plink_path <- "<plink_path>/plink2" 

Step 3.1.2 It is always recommended to check how the files look like before using them in functions, for better understanding. You may directly use the data files embedded in the package as a trial. You can check the top proportion of each data file using the following code:

head(Bphe_discovery) #phenotype file of the discovery sample when the outcome is binary
head(Bcov_discovery) #covariate file of the discovery sample when the outcome is binary
head(Bphe_target) #phenotype file of the target sample when the outcome is binary
head(Bcov_target) #covariate file of the target sample when the outcome is binary
head(Qphe_discovery) #phenotype file of the discovery sample when the outcome is quantitative
head(Qcov_discovery) #covariate file of the discovery sample when the outcome is quantitative
head(Qphe_target) #phenotype file of the target sample when the outcome is quantitative
head(Qcov_target) #covariate file of the target sample when the outcome is quantitative

Step 3.1.3 To call the data files saved in "inst" directory, you can follow the following code to obtain the path of each data file.

inst_path <- system.file(package = "GxEprs") 
DummyData <- paste0(inst_path, "/DummyData") #this contains all .fam, .bed and .bim files. They can be accessed by a direct call of prefix "DummyData"
Bphe_discovery <- paste0(inst_path, "/Bphe_discovery.txt")
Bcov_discovery <- paste0(inst_path, "/Bcov_discovery.txt")
Bphe_target <- paste0(inst_path, "/Bphe_target.txt")
Bcov_target <- paste0(inst_path, "/Bcov_target.txt")
Qphe_discovery <- paste0(inst_path, "/Qphe_discovery.txt")
Qcov_discovery <- paste0(inst_path, "/Qcov_discovery.txt")
Qphe_target <- paste0(inst_path, "/Qphe_target.txt")
Qcov_target <- paste0(inst_path, "/Qcov_target.txt")

Note that the step 3.1.3 described above is to call the embedded data files in this package itself. However, when users have to call their own data, they can follow the same approach. It is more convenient if the users can store all their data files in the same working directory. For example, assume that the file names are as follows (Refer to 'File formats' section of this document to view the formatting details of each of the following input file):

• binary files: mydata.fam, mydata.bim and mydata.bed
• phenotype file of discovery sample (binary outcome): Bpd.txt
• covariate file of discovery sample (binary outcome): Bcd.txt
• phenotype file of target sample (binary outcome): Bpt.txt
• covariate file of target sample (binary outcome): Bct.txt
• phenotype file of discovery sample (quantitative outcome): Qpd.txt
• covariate file of discovery sample (quantitative outcome): Qcd.txt
• phenotype file of target sample (quantitative outcome): Qpt.txt
• covariate file of target sample (quantitative outcome): Qct.txt

Additional note:

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;

• instead of "mydata" use "<path>/mydata" as DummyData
• instead of "Bpd.txt" use "<path>/Bpd.txt" as Bphe_discovery
• instead of "Bcd.txt" use "<path>/Bcd.txt" as Bcov_discovery
• instead of "Bpt.txt" use "<path>/Bpt.txt" as Bphe_target
• instead of "Bct.txt" use "<path>/Bct.txt" as Bcov_target
• instead of "Qpd.txt" use "<path>/Qpd.txt" as Qphe_discovery
• instead of "Qcd.txt" use "<path>/Qcd.txt" as Qcov_discovery
• instead of "Qpt.txt" use "<path>/Qpt.txt" as Qphe_target
• instead of "Qct.txt" use "<path>/Qct.txt" as Qcov_target

Step 3.1.4 Set the working directory and run the following R functions

setwd("<path to working directory>") #set the working directory

$\color{red}{NOTE:}$ Read manual.pdf document for descriptions of arguments passed for each function.

When the outcome variable is binary

Command

x <- GWAS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt", thread = 20)

As explained above, “mydata” is the prefix of the PLINK format files, “Bpd.txt” is binary phenotype file of the discovery sample, "Bcd.txt" is the covariate file of the discovery sample, thread indicates the number of CPUs used to run the command which can be optionally specified by the user (default is 20). This command performs GWAS using a logistic regression, and outputs GWAS summary statistics of all additive SNP effects.

To perform the same on embedded data

x <- GWAS_binary(plink_path, DummyData, Bphe_discovery, Bcov_discovery, thread = 20)

Output

  CHROM     POS    ID REF ALT A1 OBS_CT        BETA       SE    Z_STAT         P
1     1  768448 SNP_1   G   A  A    800 -0.09124944 0.421016 -0.216738 0.8284130
2     1  853954 SNP_2   A   C  C    800  0.58099141 0.266219  2.182390 0.0290808
3     1  880390 SNP_3   C   A  A    800  0.44618470 0.744786  0.599076 0.5491220
4     1  940203 SNP_4   G   A  A    800  0.47615964 0.456313  1.043500 0.2967180
5     1  987670 SNP_5   G   T  T    798 -0.35811312 0.540045 -0.663118 0.5072550
6     1 1021695 SNP_6   G   A  A    790 -0.63421057 0.331609 -1.912520 0.0558090

x contains GWAS summary statistics of all additive SNP effects, when the outcome is binary.

• x$CHROM : chromosome number
• x$POS : base pair position
• x$ID : SNP ID
• x$REF : reference allele
• x$ALT : alternate allele
• x$A1 : minor allele
• x$OBS_CT : number of allele observations
• x$BETA : SNP effects
• x$SE : standard errors the SNP effects
• x$Z_STAT : test statistics of SNP effects
• x$P : p values of SNP effects

See the topic GWAS_binary (page 6) in manual.pdf for examples.

Command

x <- GWEIS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt", thread = 20)

This performs GWEIS using a logistic regression, and outputs GWEIS summary statistics of both additive and interaction SNP effects.

To perform the same on embedded data

x <- GWEIS_binary(plink_path, DummyData, Bphe_discovery, Bcov_discovery, thread = 20)

Output

  CHROM     POS    ID REF ALT A1 OBS_CT    ADD_BETA   ADD_SE ADD_Z_STAT     ADD_P INTERACTION_BETA INTERACTION_SE INTERACTION_Z_STAT INTERACTION_P
1     1  768448 SNP_1   G   A  A    800 -0.04159104 0.548199  -0.075868 0.9395240       -0.1510322       0.374891          -0.402868      0.687045
2     1  853954 SNP_2   A   C  C    800  0.69288215 0.339653   2.039970 0.0413535       -0.1130265       0.209144          -0.540423      0.588905
3     1  880390 SNP_3   C   A  A    800 -0.02979757 1.283090  -0.023223 0.9814720        0.6069464       0.927640           0.654292      0.512924
4     1  940203 SNP_4   G   A  A    800  0.53499517 0.578765   0.924371 0.3552930       -0.3087518       0.386390          -0.799069      0.424251
5     1  987670 SNP_5   G   T  T    798 -0.43978175 0.755131  -0.582393 0.5603020        0.1012742       0.651603           0.155419      0.876491
6     1 1021695 SNP_6   G   A  A    790 -0.96678961 0.484625  -1.994920 0.0460513        0.1914052       0.288440           0.663593      0.506951

x contains GWEIS summary statistics of all additive and interaction SNP effects, when the outcome is binary.

• x$CHROM : chromosome number
• x$POS : base pair position
• x$ID : SNP ID
• x$REF : reference allele
• x$ALT : alternate allele
• x$A1 : minor allele
• x$OBS_CT : number of allele observations
• x$ADD_BETA : additive SNP effects
• x$ADD_SE : standard errors of the additive SNP effects
• x$ADD_Z_STAT : test statistics of additive SNP effects
• x$ADD_P : p values of additive SNP effects
• x$INTERACTION_BETA : SNP effect of interaction SNP effects
• x$INTERACTION_SE : standard errors of the interaction SNP effects
• x$INTERACTION_Z_STAT : test statistics of interaction SNP effects
• x$INTERACTION_P : p values of interaction SNP effects

See the topic GWEIS_binary (page 9) in manual.pdf for examples.

Command

a <- GWAS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]

x <- PRS_binary(plink_path, "mydata", summary_input = trd)
y <- PRS_binary(plink_path, "mydata", summary_input = add)
z <- PRS_binary(plink_path, "mydata", summary_input = gxe)

As explained above, “mydata” is the prefix of the PLINK format files, , add and gxe are summary statistics generated from previous functions and used as an input for this function to construct PRS. These commands compute polygenic risk scores for each individual in the target dataset and outputs the PRSs of all individuals.

To perform the same on embedded data

a <- GWAS_binary(plink_path, DummyData, Bphe_discovery, Bcov_discovery)
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_binary(plink_path, DummyData, Bphe_discovery, Bcov_discovery)
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]

x <- PRS_binary(plink_path, DummyData, summary_input = trd)
y <- PRS_binary(plink_path, DummyData, summary_input = add)
z <- PRS_binary(plink_path, DummyData, summary_input = gxe)

Output

   FID  IID         PRS
1 ID_1 ID_1 -0.01121010
2 ID_2 ID_2  0.00675935
3 ID_3 ID_3 -0.00961251
4 ID_4 ID_4 -0.00707935
5 ID_5 ID_5 -0.01488360
6 ID_6 ID_6  0.00223896

x contains the following columns in order.

• x$FID : family IDs of the full dataset
• x$IID : individual IDs of the full dataset
• x$PRS : polygenic risk scores (PRSs), computed from the additive effects of GWAS summary statistics, of the full dataset

   FID  IID          PRS
1 ID_1 ID_1 -7.22524e-03
2 ID_2 ID_2  3.00896e-02
3 ID_3 ID_3 -3.66145e-05
4 ID_4 ID_4 -7.42899e-03
5 ID_5 ID_5  2.75918e-03
6 ID_6 ID_6  1.10225e-02

y contains the the following columns in order.

• y$FID : family IDs of the full dataset
• y$IID : individual IDs of the full dataset
• y$PRS : polygenic risk scores (PRSs), computed from the additive effects of GWEIS summary statistics), of the full dataset

   FID  IID         PRS
1 ID_1 ID_1 -0.00567046
2 ID_2 ID_2 -0.02604500
3 ID_3 ID_3 -0.01254020
4 ID_4 ID_4 -0.00127688
5 ID_5 ID_5 -0.03479880
6 ID_6 ID_6 -0.01316010

z contains the the following columns in order.

• z$FID : family IDs of the full dataset
• z$IID : individual IDs of the full dataset
• z$PRS : polygenic risk scores (PRSs), computed from the interaction effects of GWEIS summary statistics), of the full dataset

See the topic PRS_binary (page 12) in manual.pdf for examples.

Command

a <- GWAS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_binary(plink_path, "mydata", summary_input = trd)
q <- PRS_binary(plink_path, "mydata", summary_input = add)
r <- PRS_binary(plink_path, "mydata", summary_input = gxe)

u <- summary_regular_binary("Bpt.txt", "Bct.txt", add_score = p, gxe_score = r, Model = 0) #see section IMPORTANT for details
v <- summary_regular_binary("Bpt.txt", "Bct.txt", add_score = p, Model = 1)
w <- summary_regular_binary("Bpt.txt", "Bct.txt", add_score = q, Model = 2)
x <- summary_regular_binary("Bpt.txt", "Bct.txt", add_score = q, gxe_score = r, Model = 3)
y <- summary_regular_binary("Bpt.txt", "Bct.txt", add_score = q, gxe_score = r, Model = 4)
z <- summary_regular_binary("Bpt.txt", "Bct.txt", add_score = q, gxe_score = r, Model = 5)

“Bpt.txt” is binary phenotype file of the target sample, "Bct.txt" is the covariate file of the target sample, p, q and r are the PRSs generated. If users choose to use their own PRSs stored previously, then they should create their own files (having column names as FID, IID and PRS) to use as inputs to summary_regular_binary() function (see examples.R script file attached to this repository). Depending on the model used, the input should be varied. (See section $\color{red}{IMPORTANT}$ for the target models.) This function outputs both the summary of the fitted regular model for binary outcome and all the calculated individual risk scores of the fitted regular model for binary outcome. Note that if the input for add_score is provided the same, output of fitting Model 1 and Model 2 will be equivalent. For example, summary_regular_binary("Qpt.txt", "Qct.txt", add_score = p, Model = 1) is equivalent to summary_regular_binary("Qpt.txt", "Qct.txt", add_score = p, Model = 2) where both functions will output the results after fitting Model 1 (as the input for add_score is p). Similarly, summary_regular_binary("Qpt.txt", "Qct.txt", add_score = q, Model = 1) is equivalent to summary_regular_binary("Qpt.txt", "Qct.txt", add_score = q, Model = 2) where both functions will output the results after fitting Model 2 (as the input for add_score is q).

Refer to the section $\color{red}{IMPORTANT}$ at the end of this document for details about models fitted at this step.

To perform the same on embedded data:

a <- GWAS_binary(plink_path, DummyData, Bphe_discovery, Bcov_discovery)
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_binary(plink_path, DummyData, Bphe_discovery, Bcov_discovery)
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_binary(plink_path, DummyData, summary_input = trd)
q <- PRS_binary(plink_path, DummyData, summary_input = add)
r <- PRS_binary(plink_path, DummyData, summary_input = gxe)

u <- summary_regular_binary(Bphe_target, Bcov_target, add_score = p, gxe_score = r, Model = 0) #see section IMPORTANT for details
v <- summary_regular_binary(Bphe_target, Bcov_target, add_score = p, Model = 1)
w <- summary_regular_binary(Bphe_target, Bcov_target, add_score = q, Model = 2)
x <- summary_regular_binary(Bphe_target, Bcov_target, add_score = q, gxe_score = r, Model = 3)
y <- summary_regular_binary(Bphe_target, Bcov_target, add_score = q, gxe_score = r, Model = 4)
z <- summary_regular_binary(Bphe_target, Bcov_target, add_score = q, gxe_score = r, Model = 5)

Note: For demonstration, we used Model = 5:

Output

                Coefficient  Std.Error Test.Statistic     pvalue
(Intercept)     -1.02602149 1.00854092    -1.01733254 0.30899526
E                0.10336521 0.43227954     0.23911659 0.81101518
`E squared`      0.04255889 0.44343401     0.09597569 0.92353988
PRS_add         -1.33954526 0.92174446    -1.45327183 0.14614830
PRS_gxe         -1.80677738 1.20370678    -1.50101121 0.13335266
`PRS_gxe x E`    0.31122769 0.59556423     0.52257619 0.60126921
`Confounder 1`   0.21383426 0.26910611     0.79460948 0.42684070
`Confounder 2`  -0.34560722 0.34255458    -1.00891140 0.31301713
`Confounder 3`   0.39600696 0.37262842     1.06273955 0.28790008
`Confounder 4`   0.18857671 0.32469004     0.58078993 0.56138204
`Confounder 5`   0.16670465 0.32441357     0.51386459 0.60734667
`Confounder 6`  -0.26352081 0.48246734    -0.54619408 0.58493254
`Confounder 7`   0.38994824 0.52679490     0.74022782 0.45916177
`Confounder 8`   0.14785002 0.34954545     0.42297796 0.67231134
`Confounder 9`   0.46205534 0.32386129     1.42670750 0.15366422
`Confounder 10`  0.10277880 0.36177449     0.28409631 0.77633658
`Confounder 11`  0.02771563 0.37190774     0.07452286 0.94059435
`Confounder 12`  0.17135167 0.37399582     0.45816466 0.64683415
`Confounder 13`  1.15950406 0.68751357     1.68651808 0.09169606
`Confounder 14` -0.20112706 0.08172877    -2.46090890 0.01385856

z$summary contains the target regular model summary output, when the outcome is binary.

• z$summary[,"Coefficient"] : regression coefficients of all model components
• z$summary[,"Std.Error"] : standard errors of all regression coefficients
• z$summary[,"Test.Statistic"] : test statistic value of all regression coefficients
• z$summary[,"pvalue"] : p-value of all regression coefficients

  FID      IID      Risk.Values          
1 "ID_801" "ID_801" "0.32863562946532"   
2 "ID_802" "ID_802" "0.178604626660152"  
3 "ID_803" "ID_803" "0.00732496962395319"
4 "ID_804" "ID_804" "0.285968921409738"  
5 "ID_805" "ID_805" "0.0699590517867513" 
6 "ID_806" "ID_806" "0.0420642156888637" 

z$risk.values contains all the calculated individual risk scores using the target dataset (e.g. Model 5), when the outcome is binary. The columns denote the following in order.

• z$risk.values[,"FID"] : family IDs of the target dataset
• z$risk.values[,"IID"] : individual IDs of the target dataset
• z$risk.values[,"Risk.Values] : estimated risk values of the target dataset

Note: It is recommended to fit both regular and permuted models and obtain the summary of both fitted models (using summary_regular_binary("Bpt.txt", "Bct.txt", add_score = q, gxe_score = r, Model = 5) and summary_permuted_binary("Bpt.txt", "Bct.txt", iterations = 1000, add_score = q, gxe_score = r)), if you choose to fit 'PRS_gxe x E' interaction component (i.e. novel proposed model, Model 5) when generating risk scores. If the 'PRS_gxe x E' term is significant in Model 5, and insignificant in Model 5* (permuted p value), consider that the 'PRS_gxe x E' interaction component is actually insignificant (always give priority to the p value obtained from the permuted model). See Figure 5 in Jayasinghe et.al. for simulation verification.

See the topic summary_regular_binary (page 19) in manual.pdf for examples.

Command

a <- GWEIS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
add <- a[c("ID", "A1", "ADD_BETA")]
gxe <- a[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_binary(plink_path, "mydata", summary_input = add)
q <- PRS_binary(plink_path, "mydata", summary_input = gxe)
x <- summary_permuted_binary("Bpt.txt", "Bct.txt", iterations = 1000, add_score = p, gxe_score = q)

To perform the same on embedded data:

a <- GWEIS_binary(plink_path, DummyData, Bphe_discovery, Bcov_discovery)
add <- a[c("ID", "A1", "ADD_BETA")]
gxe <- a[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_binary(plink_path, DummyData, summary_input = add)
q <- PRS_binary(plink_path, DummyData, summary_input = gxe)
x <- summary_permuted_binary(Bphe_target, Bcov_target, iterations = 1000, add_score = p, gxe_score = q)

x outputs the p value of the fitted permuted model for binary outcome.

When the outcome variable is quantitative

Command

x <- GWAS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt", thread = 20)

As explained above, “mydata” is the prefix of the PLINK format files, “Qpd.txt” is quantitative phenotype file of the discovery sample, "Qcd.txt" is the covariate file of the discovery sample, thread indicates the number of CPUs used to run the command which can be optionally specified by the user (default is 20). This command performs GWAS using a linear regression, and outputs GWAS summary statistics of all additive SNP effects.

To perform the same on embedded data:

x <- GWAS_quantitative(plink_path, DummyData, Qphe_discovery, Qcov_discovery, thread = 20)

Output

  CHROM     POS    ID REF ALT A1 OBS_CT        BETA        SE     T_STAT         P
1     1  768448 SNP_1   G   A  A    800  0.04087160 0.0800126  0.5108150 0.6096240
2     1  853954 SNP_2   A   C  C    800 -0.07500170 0.0522118 -1.4364900 0.1512620
3     1  880390 SNP_3   C   A  A    800  0.03102760 0.1522040  0.2038550 0.8385200
4     1  940203 SNP_4   G   A  A    800 -0.19965600 0.1024470 -1.9488800 0.0516666
5     1  987670 SNP_5   G   T  T    798 -0.00115168 0.0940650 -0.0122435 0.9902350
6     1 1021695 SNP_6   G   A  A    790  0.04444630 0.0545538  0.8147240 0.4154810

x contains GWAS summary statistics of all additive SNP effects, when the outcome is quantitative. V1 to V13 denote the following columns in order. Note that all summary statistics follow the same structure.

• x$CHROM : chromosome number
• x$POS : base pair position
• x$ID : SNP ID
• x$REF : reference allele
• x$ALT : alternate allele
• x$A1 : minor allele
• x$OBS_CT : number of allele observations
• x$BETA : SNP effects
• x$SE : standard errors of SNP effects
• x$T_STAT : test statistics of SNP effects
• x$P : p values of SNP effects

See the topic GWAS_quantitative (page 7) in manual.pdf for examples.

Command

x <- GWEIS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt", thread = 20)

This performs GWEIS using a linear regression, and outputs GWEIS summary statistics of all additive and interaction SNP effects.

To perform the same on embedded data

x <- GWEIS_quantitative(plink_path, DummyData, Qphe_discovery, Qcov_discovery, thread = 20)

Output

  CHROM     POS    ID REF ALT A1 OBS_CT    ADD_BETA    ADD_SE ADD_T_STAT     ADD_P INTERACTION_BETA INTERACTION_SE INTERACTION_T_STAT INTERACTION_P
1     1  768448 SNP_1   G   A  A    800  0.03435200 0.0805291  0.4265790 0.6698030        0.0771699      0.0818123           0.943256     0.3458410
2     1  853954 SNP_2   A   C  C    800 -0.07503330 0.0523249 -1.4339900 0.1519750        0.0145521      0.0507679           0.286640     0.7744640
3     1  880390 SNP_3   C   A  A    800  0.05319270 0.1622820  0.3277790 0.7431670        0.0804831      0.1953130           0.412073     0.6803990
4     1  940203 SNP_4   G   A  A    800 -0.20973600 0.1040870 -2.0150000 0.0442465       -0.0579509      0.1064690          -0.544299     0.5863910
5     1  987670 SNP_5   G   T  T    798 -0.00105564 0.0939844 -0.0112321 0.9910410       -0.1747820      0.0934860          -1.869610     0.0619128
6     1 1021695 SNP_6   G   A  A    790  0.04103200 0.0546241  0.7511700 0.4527790       -0.0805943      0.0565680          -1.424730     0.1546380

x contains GWEIS summary statistics of all additive and interaction SNP effects, when the outcome is quantitative.

• x$CHROM : chromosome number
• x$POS : base pair position
• x$ID : SNP ID
• x$REF : reference allele
• x$ALT : alternate allele
• x$A1 : minor allele
• x$OBS_CT : number of allele observations
• x$ADD_BETA : effects of additive SNP effects
• x$ADD_SE : standard errors of the effects of additive SNP effects
• x$ADD_T_STAT : test statistics of additive SNP effects
• x$ADD_P : p values of additive SNP effects
• x$INTERACTION_BETA : interaction SNP effects
• x$INTERACTION_SE : standard errors of interaction SNP effects
• x$INTERACTION_T_STAT : test statistics of interaction SNP effects
• x$INTERACTION_P : p values of interaction SNP effects

See the topic GWEIS_quantitative (page 10) in manual.pdf for examples.

Command

a <- GWAS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt")
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt")
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]

x <- PRS_quantitative(plink_path, "mydata", summary_input = trd)
y <- PRS_quantitative(plink_path, "mydata", summary_input = add)
z <- PRS_quantitative(plink_path, "mydata", summary_input = gxe)

As explained above, “mydata” is the prefix of the PLINK format files, , add and gxe are summary statistics generated from previous functions and used as an input for this function to construct PRS. These commands compute polygenic risk scores for each individual in the target dataset and outputs the PRSs of all individuals.

To perform the same on embedded data

a <- GWAS_quantitative(plink_path, DummyData, Qphe_discovery, Qcov_discovery)
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_quantitative(plink_path, DummyData, Qphe_discovery, Qcov_discovery)
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]

x <- PRS_quantitative(plink_path, DummyData, summary_input = trd)
y <- PRS_quantitative(plink_path, DummyData, summary_input = add)
z <- PRS_quantitative(plink_path, DummyData, summary_input = gxe)

Output

   FID  IID          PRS
1 ID_1 ID_1 -9.72106e-04
2 ID_2 ID_2  2.69304e-04
3 ID_3 ID_3  3.52687e-04
4 ID_4 ID_4 -9.61218e-05
5 ID_5 ID_5  1.17225e-03
6 ID_6 ID_6  1.33232e-03

x contains the following columns in order.

• x$FID : family IDs of the full dataset
• x$IID : individual IDs of the full dataset
• x$PRS : polygenic risk scores (PRSs), computed from the additive effects of GWAS summary statistics, of the full dataset

   FID  IID          PRS
1 ID_1 ID_1 -9.56472e-04
2 ID_2 ID_2  3.24770e-04
3 ID_3 ID_3  3.55841e-04
4 ID_4 ID_4 -9.98136e-06
5 ID_5 ID_5  1.28597e-03
6 ID_6 ID_6  1.38409e-03

y contains the the following columns in order.

• y$FID : family IDs of the full dataset
• y$IID : individual IDs of the full dataset
• y$PRS : polygenic risk scores (PRSs), computed from the additive effects of GWEIS summary statistics), of the full dataset

   FID  IID          PRS
1 ID_1 ID_1  3.99456e-04
2 ID_2 ID_2 -3.80288e-04
3 ID_3 ID_3  5.48983e-05
4 ID_4 ID_4 -3.20377e-04
5 ID_5 ID_5  1.90810e-03
6 ID_6 ID_6 -1.97507e-05

z contains the the following columns in order.

• z$FID : family IDs of the full dataset
• z$IID : individual IDs of the full dataset
• z$PRS : polygenic risk scores (PRSs), computed from the interaction effects of GWEIS summary statistics), of the full dataset

See the topic PRS_quantitative (page 13) in manual.pdf for examples.

Command

a <- GWAS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt")
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt")
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_quantitative(plink_path, "mydata", summary_input = trd)
q <- PRS_quantitative(plink_path, "mydata", summary_input = add)
r <- PRS_quantitative(plink_path, "mydata", summary_input = gxe)

u <- summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = p, gxe_score = r, Model = 0) #see section IMPORTANT for details
v <- summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = p, Model = 1)
w <- summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = q, Model = 2)
x <- summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = q, gxe_score = r, Model = 3)
y <- summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = q, gxe_score = r, Model = 4)

“Qpt.txt” is quantitative phenotype file of the target sample, "Qct.txt" is the covariate file of the target sample, p, q and r are the PRSs generated. If users choose to use their own PRSs stored previously, then they should create their own files (having column names as FID, IID and PRS) to use as inputs to summary_regular_quantitative() function (see examples.R script file attached to this repository). Depending on the model used, the input should be varied. (See section $\color{red}{IMPORTANT}$ for the target models.) This function outputs both the summary of the fitted regular model for quantitative outcome and all the calculated individual risk scores of the fitted regular model for quantitative outcome. Note that if the input for add_score is provided the same, output of fitting Model 1 and Model 2 will be equivalent. For example, summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = p, Model = 1) is equivalent to summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = p, Model = 2) where both functions will output the results after fitting Model 1 (as the input for add_score is p). Similarly, summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = q, Model = 1) is equivalent to summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = q, Model = 2) where both functions will output the results after fitting Model 2 (as the input for add_score is q).

Refer to the section $\color{red}{IMPORTANT}$ at the end of this document for details about models fitted at this step.

To perform the same on embedded data:

a <- GWAS_quantitative(plink_path, DummyData, Qphe_discovery, Qcov_discovery)
trd <- a[c("ID", "A1", "BETA")]
b <- GWEIS_quantitative(plink_path, DummyData, Qphe_discovery, Qcov_discovery)
add <- b[c("ID", "A1", "ADD_BETA")]
gxe <- b[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_quantitative(plink_path, DummyData, summary_input = trd)
q <- PRS_quantitative(plink_path, DummyData, summary_input = add)
r <- PRS_quantitative(plink_path, DummyData, summary_input = gxe)

u <- summary_regular_quantitative(Qphe_target, Qcov_target, add_score = p, gxe_score = r, Model = 0) #see section IMPORTANT for details
v <- summary_regular_quantitative(Qphe_target, Qcov_target, add_score = p, Model = 1)
w <- summary_regular_quantitative(Qphe_target, Qcov_target, add_score = q, Model = 2)
x <- summary_regular_quantitative(Qphe_target, Qcov_target, add_score = q, gxe_score = r, Model = 3)
y <- summary_regular_quantitative(Qphe_target, Qcov_target, add_score = q, gxe_score = r, Model = 4)

Note: For demonstration, we used Model = 4:

Output

                Coefficient  Std.Error Test.Statistic      pvalue
(Intercept)      0.13646629 0.23224944      0.5875850 0.557543186
E               -0.01197672 0.07255329     -0.1650747 0.869069369
PRS_add         -0.04723272 0.07310681     -0.6460783 0.519046949
PRS_gxe         -0.02954316 0.07596885     -0.3888852 0.697818196
`PRS_gxe x E`   -0.01731704 0.07376726     -0.2347525 0.814666223
`Confounder 1`  -0.17016325 0.07040773     -2.4168261 0.016649362
`Confounder 2`  -0.05331027 0.07294422     -0.7308362 0.465824222
`Confounder 3`   0.07103913 0.07787459      0.9122247 0.362863822
`Confounder 4`   0.12420393 0.07342810      1.6915039 0.092461587
`Confounder 5`  -0.12654205 0.07743131     -1.6342491 0.103944656
`Confounder 6`  -0.08223335 0.09958450     -0.8257646 0.410024755
`Confounder 7`   0.13187809 0.09925324      1.3287032 0.185618601
`Confounder 8`  -0.07076306 0.07489718     -0.9448027 0.346019099
`Confounder 9`  -0.09654830 0.07574628     -1.2746276 0.204074978
`Confounder 10` -0.08690212 0.08328035     -1.0434888 0.298113621
`Confounder 11` -0.02116981 0.07955715     -0.2660956 0.790468452
`Confounder 12`  0.02804831 0.08793378      0.3189708 0.750116468
`Confounder 13`  0.44386988 0.14215406      3.1224566 0.002088447
`Confounder 14` -0.02322767 0.01469411     -1.5807462 0.115681987

y$summary contains the target regular model summary output, when the outcome is quantitative.

• y$summary[,"Coefficient"] : regression coefficients of all model components
• y$summary[,"Std.Error"] : standard errors of all regression coefficients
• y$summary[,"Test.Statistic"] : test statistic value of all regression coefficients
• y$summary[,"pvalue"] : p-value of all regression coefficients

  FID      IID      Risk.Values          
1 "ID_801" "ID_801" "-0.250074533388036" 
2 "ID_802" "ID_802" "-0.0219035934167711"
3 "ID_803" "ID_803" "0.18809813264688"   
4 "ID_804" "ID_804" "0.214732793779811"  
5 "ID_805" "ID_805" "0.314792745880588"  
6 "ID_806" "ID_806" "0.335673297263822"  

y$risk.values contains all the calculated individual risk scores using the target dataset (e.g. Model 4), when the outcome is quantitative. The columns denote the following in order.

• y$risk.values[,"FID"] : family IDs of the target dataset
• y$risk.values[,"IID"] : individual IDs of the target dataset
• y$risk.values[,"Risk.Values] : estimated risk values of the target dataset

Note: It is recommended to fit both regular and permuted models and obtain the summary of both fitted models (using summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = q, gxe_score = r, Model = 4) and summary_permuted_quantitative("Qpt.txt", "Qct.txt", iterations = 1000, add_score = q, gxe_score = r)), if you choose to fit 'PRS_gxe x E' interaction component (i.e. novel proposed model, Model 4) when generating risk scores. If the 'PRS_gxe x E' term is significant in Model 4, and insignificant in Model 4* (permuted p value), consider that the 'PRS_gxe x E' interaction component is actually insignificant (always give priority to the p value obtained from the permuted model). See Figure 3 in Jayasinghe et. al. for simulation verification.

See the topic summary_regular_quantitative (page 21) in manual.pdf for examples.

Command

a <- GWEIS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt")
add <- a[c("ID", "A1", "ADD_BETA")]
gxe <- a[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_quantitative(plink_path, "mydata", summary_input = add)
q <- PRS_quantitative(plink_path, "mydata", summary_input = gxe)
x <- summary_permuted_quantitative("Qpt.txt", "Qct.txt", iterations = 1000, add_score = p, gxe_score = q)

To perform the same on embedded data:

a <- GWEIS_quantitative(plink_path, DummyData, Qphe_discovery, Qcov_discovery)
add <- a[c("ID", "A1", "ADD_BETA")]
gxe <- a[c("ID", "A1", "INTERACTION_BETA")]
p <- PRS_quantitative(plink_path, DummyData, summary_input = add)
q <- PRS_quantitative(plink_path, DummyData, summary_input = gxe)
x <- summary_permuted_quantitative(Qphe_target, Qcov_target, iterations = 1000, add_score = p, gxe_score = q)

x outputs the p value of the fitted permuted model for quantitative outcome.

$$\color{red}{IMPORTANT}$$

The discovery model used in GWAS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt", thread = 20) or GWAS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt", thread = 20) is as follows:

• y = b_trd.W + error where y is the outcome variable, b_trd is the estimated SNP effect and W is the SNP genotype.

The discovery model used in GWEIS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt", thread = 20) is as follows:

• y = b_add.W + b_cov.E + b_cov2.E^2 + b_gxe.(WxE) + error where y is the outcome variable, b_add is the estimated additive SNP effect, E is the covariate, W is the SNP genotype, b_cov is the estimated effect of the covariate, b_cov2 is the estimated effect of the squared covariate and b_gxe is the estimated effect of the Hadamard product of WxE.

The discovery model used in GWEIS_quantitative(plink_path, "mydata", "Qpd.txt", "Qcd.txt", thread = 20) is as follows:

• y = b_add.W + b_cov.E + b_gxe.(WxE) + error where y is the outcome variable, b_add is the estimated additive SNP effect, E is the covariate, W is the SNP genotype, b_cov is the estimated effect of the covariate and b_gxe is the estimated effect of the Hadamard product of WxE.

The fitted (target) models in summary_regular_binary("Bpt.txt", "Bct.txt", add_score = NULL, gxe_score = NULL, Model) or summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = NULL, gxe_score = NULL, Model) are as follows:

• Model 1: y = PRS_trd + E + PRS_trd x E + confounders + error
• Model 2: y = PRS_add + E + PRS_add x E + confounders + error
• Model 3: y = PRS_add + E + PRS_gxe x E + confounders + error
• Model 4: y = PRS_add + E + PRS_gxe + PRS_gxe x E + confounders + error
• Model 4*: permuted Model 4
• Model 5: y = PRS_add + E + E^2 + PRS_gxe + PRS_gxe x E + confounders + error
• Model 5*: permuted Model 5

where y is the outcome variable, E is the covariate of interest, PRS_trd and PRS_add are the polygenic risk scores computed using additive SNP effects of GWAS and GWEIS summary statistics respectively, and PRS_gxe is the polygenic risk scores computed using GxE interaction SNP effects of GWEIS summary statistics.

• Model 0: We denote the reduced model of each of the above models as Model 0 in general. For example, the reduced model of Model 1 will be y = PRS_trd + E + confounders + error. In case if the users are not interested to fit the PRSxE component they can use Model 0. By this approach, they can use this option to make comparisons across each full model and reduced model, to evaluate the significance of the PRSxE component of the model of interest (given that they are nested models). Additionally, if the users are interested to fit only the additive model (for example y = PRS_trd + confounders + error) they can use Model 0 option available in summary_regular_binary("Bpt.txt", "Bct.txt", add_score = NULL, gxe_score = NULL, Model = 0) or summary_regular_quantitative("Qpt.txt", "Qct.txt", add_score = NULL, gxe_score = NULL, Model = 0). However, in this instance, users are required to provide the input files appropriately.

Example 1:

#Model 0: y = PRS_trd + confounders + error and assume y is a binary trait. Note that the data columns in covariate file for confounders are optional. See 7. in "Input Files".

a <- GWAS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
trd <- a[c("ID", "A1", "BETA")]

p <- PRS_binary(plink_path, "mydata", summary_input = trd)

r <- as.data.frame(cbind(p$FID, p$IID, p$PRS/p$PRS)) #this will create columns of FID, IID and 1's. This is passed to the argument gxe_score in summary_regular_binary() function, to make the unwanted PRSxE variable neutral in the model
colnames(r) <- c("FID", "IID", "PRS") #to make the correct file format


Bct <- read.table("Bct.txt", header = FALSE) #read the covariate file
Bct$V3 <- rep(1, length(Bct$V3)) #the 3rd column of Bct.txt file is the covariate term. Then assign a vector of same values (e.g. 1) to make the unwanted E variable neutral in the model
Bct$V4 <- rep(1, length(Bct$V4)) #the 4th column of Bct.txt file is the covariate term. Then assign a vector of same values (e.g. 1) to make the unwanted E^2 variable neutral in the model
sink("Bct_new.txt") #create a new covariate file where E and E^2 have all 1's
write.table(Bct, row.names=F, col.names=F, quote=F)
sink()

y0 <- summary_regular_binary("Bpt.txt", "Bct_new.txt", add_score = p, gxe_score = r, Model = 0) #finally fit Model 0

Example 2:

#Model 0: y = PRS_add + PRS_gxe + confounders + error and assume y is a binary trait. Note that the data columns in covariate file for confounders are optional. See 7. in "Input Files".

a <- GWEIS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
add <- a[c("ID", "A1", "ADD_BETA")]
gxe <- a[c("ID", "A1", "INTERACTION_BETA")]

p <- PRS_binary(plink_path, "mydata", summary_input = add)
r <- PRS_binary(plink_path, "mydata", summary_input = gxe)

Bct <- read.table("Bct.txt", header = FALSE) #read the covariate file
Bct$V3 <- rep(1, length(Bct$V3)) #the 3rd column of Bct.txt file is the covariate term. Then assign a vector of same values (e.g. 1) to make the unwanted E variable neutral in the model
Bct$V4 <- rep(1, length(Bct$V4)) #the 4th column of Bct.txt file is the covariate term. Then assign a vector of same values (e.g. 1) to make the unwanted E^2 variable neutral in the model
sink("Bct_new.txt") #create a new covariate file where E and E^2 have all 1's
write.table(Bct, row.names=F, col.names=F, quote=F)
sink()

y0 <- summary_regular_binary("Bpt.txt", "Bct_new.txt", add_score = p, gxe_score = r, Model = 0) #finally fit Model 0

Example 3:

#Model 0: y = PRS_add + PRS_gxe + E + E^2 confounders + error and assume y is a binary trait. Note that the data columns in covariate file for confounders are optional. See 7. in "Input Files".

a <- GWEIS_binary(plink_path, "mydata", "Bpd.txt", "Bcd.txt")
add <- a[c("ID", "A1", "ADD_BETA")]
gxe <- a[c("ID", "A1", "INTERACTION_BETA")]

p <- PRS_binary(plink_path, "mydata", summary_input = add)
r <- PRS_binary(plink_path, "mydata", summary_input = gxe)

y0 <- summary_regular_binary("Bpt.txt", "Bct.txt", add_score = p, gxe_score = r, Model = 0) #finally fit Model 0

To address the potential issue of spurious GxE signals that may arise from unknown relationships between the quantitative outcome variable (y) and covariate (E), we implemented a permutation procedure on the term PRS_gxe in the interaction component in Model 4 and denoted it as Model 4*. The purpose of this permutation was to maintain the correlation structure between the outcome variable and other model components while specifically permuting the interaction component. Similarly, in Model 5 for binary outcome, we performed the same permutation on the PRS_gxe term in the interaction component and denoted it as Model 5*. It’s important to note that the number of permutations performed in Models 4* or 5* was determined based on the p value obtained in Models 4 or 5 to ensure an adequate number of permutations (with a minimum of 1000). By selectively permuting the interaction term while preserving the correlation structure, we aimed to assess the robustness of the GxE effects in our models and mitigate the risk of spurious signals. This approach allowed us to obtain more reliable results and evaluate the significance of the observed GxE interactions.

Contact

Please contact Dovini Jayasinghe (dovini.jayasinghe@mymail.unisa.edu.au), Md Moksedul Momin (md_moksedul.momin@mymail.unisa.edu.au) or Hong Lee (hong.lee@unisa.edu.au) for queries.

Reference

Jayasinghe, D., Momin, M. M., Beckmann, K., Hyppönen, E., Benyamin, B., & Lee, S. H. (2024). Mitigating type 1 error inflation and power loss in GxE PRS: Genotype–environment interaction in polygenic risk score models. Genetic Epidemiology.