Answer sheet for the tutorial
Setting up environment
mkdir -p day7
cd day7
cp -r /faculty/sarah/2022/Day7/* .PRS calculation using PRSice
# PRSice code - run this in the ssh terminal window
Rscript /usr/local/bin/PRSice.R --dir . --prsice /usr/local/bin/PRSice --base weights.prs --target ozbmi --thread 1 --snp SNP --chr CHR --bp BP --A1 A1 --A2 A2 --stat BETA --pvalue P --beta --pheno ozht.pheno --score sum --all-score --bar-levels 0.00000005,0.000001,0.0001,0.01,1 --fastscore --quantile 5 --quant-ref 1 --binary-target F --no-regressOpen tables and merge
pheno <- read.table("ozht.pheno", header = T)
covar <- read.table("ozht.cov", header = T)
PRS <- read.table("PRSice.all_score", header = T)
head(PRS)names(PRS) <- c("FID", "IID", "PRS1", "PRS2", "PRS3", "PRS4",
"PRS5")
head(PRS)pheno$ht <- pheno$ht * 100
# Best to use names instead of column number
temp = merge(pheno, covar[, 2:8], by = "IID")
longData = merge(temp, PRS[, 2:7], by = "IID")
# Run models
nullModel <- lm(ht ~ age + sex + PC1 + PC2 + PC3 + PC4, data = longData)
summary(nullModel)##
## Call:
## lm(formula = ht ~ age + sex + PC1 + PC2 + PC3 + PC4, data = longData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.7268 -6.0023 -0.2827 5.3050 30.4879
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 173.50935 6.65163 26.085 < 2e-16 ***
## age -0.10116 0.01697 -5.961 2.99e-09 ***
## sex -9.80000 0.43327 -22.619 < 2e-16 ***
## PC1 744.83641 373.94861 1.992 0.0465 *
## PC2 240.08216 230.12319 1.043 0.2970
## PC3 -97.14299 151.21233 -0.642 0.5207
## PC4 54.75009 145.82996 0.375 0.7074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.412 on 1887 degrees of freedom
## (38 observations deleted due to missingness)
## Multiple R-squared: 0.2296, Adjusted R-squared: 0.2271
## F-statistic: 93.71 on 6 and 1887 DF, p-value: < 2.2e-16# Then we are going to run a model that includes the
# covariates and the first PRS
thresh1 <- lm(ht ~ age + sex + PC1 + PC2 + PC3 + PC4 + PRS1,
data = longData)
summary(thresh1)##
## Call:
## lm(formula = ht ~ age + sex + PC1 + PC2 + PC3 + PC4 + PRS1, data = longData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.6273 -5.9917 -0.2311 5.5705 29.9557
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 175.28775 6.60888 26.523 < 2e-16 ***
## age -0.10084 0.01684 -5.988 2.54e-09 ***
## sex -9.76070 0.43002 -22.698 < 2e-16 ***
## PC1 667.06438 371.36819 1.796 0.0726 .
## PC2 266.81843 228.42069 1.168 0.2429
## PC3 -96.05118 150.05960 -0.640 0.5222
## PC4 35.77279 144.75946 0.247 0.8048
## PRS1 3.11956 0.56855 5.487 4.64e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.348 on 1886 degrees of freedom
## (38 observations deleted due to missingness)
## Multiple R-squared: 0.2417, Adjusted R-squared: 0.2389
## F-statistic: 85.86 on 7 and 1886 DF, p-value: < 2.2e-16# We can compute the difference in R2 between the 2 models
# using the following code
summary(thresh1)$r.squared - summary(nullModel)$r.squared## [1] 0.01210506Format data - to one row per family (long -> wide)
zyg <- read.table("ozht.zyg", header = T)
# We will then merge on the zyg data using the individual
# ID 'IID' as the key.
longData1 = merge(longData, zyg[, 2:5], by = "IID")
# Take a look at the result to check that it worked
head(longData1)# For openMx we need our data to be in a wide format with
# one line per family so we will need to reshape the data.
# In this code the variables listed in v.names are present
# for both twins The variable Twin specifies if the person
# is twin 1 or 2
wideData = reshape(data = longData1, direction = "wide", v.names = c("IID",
"ht", "age", "sex", "PC1", "PC2", "PC3", "PC4", "PRS1", "PRS2",
"PRS3", "PRS4", "PRS5"), idvar = "FID", timevar = "Twin",
sep = "_")
# We will replace all missingness in the definition
# variables with 0 the people with missing definition
# variables don't have phenotypic data so this will not
# bais estimates
wideData[, c(6:12, 19:25)][is.na(wideData[, c(6:12, 19:25)])] <- 0OpenMx model
library(OpenMx)
source("miFunctions5272022.R")
# Select Variables for Analysis
vars <- "ht" # list of variables names
nv <- 1 # number of variables
ntv <- nv * 2 # number of total variables
depVar <- c("ht_1", "ht_2") # dependent variables
indVar <- c("age_1", "sex_1", "PC1_1", "PC2_1", "PC3_1", "PC4_1",
"PRS1_1", "age_2", "sex_2", "PC1_2", "PC2_2", "PC3_2", "PC4_2",
"PRS1_2") # independent variables or predictors
selVars <- c(depVar, indVar)
# Select Data for Analysis
mzData <- subset(wideData, MZDZ == 1, selVars)
dzData <- subset(wideData, MZDZ == 2, selVars)
# Set Starting Values
svMe <- 170 # start value for means
svVa <- 100 # start value for variance
lbVa <- 1e-05 # lower bound for variance
svBe <- 1e-04 # start value for means
# Create Data Objects for Multiple Groups
dataMZ <- mxData(observed = mzData, type = "raw")
dataDZ <- mxData(observed = dzData, type = "raw")
# Objects to hold definition variables for the regression
dPC1 <- mxMatrix(type = "Full", nrow = 1, ncol = 2, free = FALSE,
labels = c("data.PC1_1", "data.PC1_2"), name = "dPC1")
dPC2 <- mxMatrix(type = "Full", nrow = 1, ncol = 2, free = FALSE,
labels = c("data.PC2_1", "data.PC2_2"), name = "dPC2")
dPC3 <- mxMatrix(type = "Full", nrow = 1, ncol = 2, free = FALSE,
labels = c("data.PC3_1", "data.PC3_2"), name = "dPC3")
dPC4 <- mxMatrix(type = "Full", nrow = 1, ncol = 2, free = FALSE,
labels = c("data.PC4_1", "data.PC4_2"), name = "dPC4")
dsex <- mxMatrix(type = "Full", nrow = 1, ncol = 2, free = FALSE,
labels = c("data.sex_1", "data.sex_2"), name = "dsex")
dage <- mxMatrix(type = "Full", nrow = 1, ncol = 2, free = FALSE,
labels = c("data.age_1", "data.age_2"), name = "dage")
dPRS <- mxMatrix(type = "Full", nrow = 1, ncol = 2, free = FALSE,
labels = c("data.PRS1_1", "data.PRS1_2"), name = "dPRS")
# Objects to hold betas for the covariates
bPC1 <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svBe, label = "b1", name = "bPC1")
bPC2 <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svBe, label = "b2", name = "bPC2")
bPC3 <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svBe, label = "b3", name = "bPC3")
bPC4 <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svBe, label = "b4", name = "bPC4")
bsex <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svBe, label = "b5", name = "bsex")
bage <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svBe, label = "b6", name = "bage")
bPRS <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svBe, label = "b7", name = "bPRS")
# Then we will build the model for the means that includes
# the regression tems (we are using Kronecker products here
# to mutliply the values of the independent variables by
# the betas)
meanG <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = TRUE,
values = svMe, labels = labVars("mean", vars), name = "meanG")
ExpMean <- mxAlgebra(expression = (meanG + bPC1 %x% dPC1 + bPC2 %x%
dPC2 + bPC3 %x% dPC3 + bPC4 %x% dPC4 + bsex %x% dsex + bage %x%
dage + bPRS %x% dPRS), name = "expMean")
# and do some house keeping (making a list of the objects
# so we can pull these into the model later)
pars <- c(bPC1, bPC2, bPC3, bPC4, bsex, bage, bPRS, meanG)
defs <- c(dPC1, dPC2, dPC3, dPC4, dsex, dage, dPRS)
# Create Algebra for expected Variance/Covariance Matrices
covMZ <- mxMatrix(type = "Symm", nrow = ntv, ncol = ntv, free = TRUE,
values = c(0.1, 0.05, 0.1), lbound = valDiag(lbVa, ntv),
labels = c("var", "cMZ21", "var"), name = "covMZ")
covDZ <- mxMatrix(type = "Symm", nrow = ntv, ncol = ntv, free = TRUE,
values = c(0.1, 0.05, 0.1), lbound = valDiag(lbVa, ntv),
labels = c("var", "cDZ21", "var"), name = "covDZ")
# Next we build the expection and model objects Create
# Expectation Objects for Multiple Groups
expMZ <- mxExpectationNormal(covariance = "covMZ", means = "expMean",
dimnames = depVar)
expDZ <- mxExpectationNormal(covariance = "covDZ", means = "expMean",
dimnames = depVar)
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
modelMZ <- mxModel(pars, defs, ExpMean, covMZ, dataMZ, expMZ,
funML, name = "MZ")
modelDZ <- mxModel(pars, defs, ExpMean, covDZ, dataDZ, expDZ,
funML, name = "DZ")
multi <- mxFitFunctionMultigroup(c("MZ", "DZ"))
# Build Saturated Model
model_wPRS <- mxModel("wPRS", pars, modelMZ, modelDZ, multi)
# Run model
fit_wPRS <- mxRun(model_wPRS)
(sum_wPRS <- summary(fit_wPRS))## Summary of wPRS
##
## free parameters:
## name matrix row col Estimate Std.Error A lbound ubound
## 1 b1 bPC1 1 1 494.07780228 460.780896671
## 2 b2 bPC2 1 1 376.45026775 275.668856990
## 3 b3 bPC3 1 1 -52.34236364 170.976227797
## 4 b4 bPC4 1 1 52.11647033 164.470117427
## 5 b5 bsex 1 1 -8.91504139 0.557884837
## 6 b6 bage 1 1 -0.11091751 0.021763799
## 7 b7 bPRS 1 1 2.58520600 0.690176704
## 8 meanht meanG 1 1 173.54694174 8.320062467
## 9 var MZ.covMZ ht_1 ht_1 74.34378646 2.844644713 1e-05
## 10 cMZ21 MZ.covMZ ht_1 ht_2 69.45798062 2.857977923 0
## 11 cDZ21 DZ.covDZ ht_1 ht_2 16.93968677 3.188544188 0
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 11 1883 12435.322
## Saturated: NA NA NA
## Independence: NA NA NA
## Number of observations/statistics: 969/1894
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: 8669.32237 12457.322 12457.598
## BIC: -512.68389 12510.961 12476.025
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2022-06-15 00:48:26
## Wall clock time: 4.534725 secs
## optimizer: NPSOL
## OpenMx version number: 2.19.5.1
## Need help? See help(mxSummary)# And then run a submodel whic drops the PRS effect from
# the model Drop the PRS from the model (running the null
# model)
model_Null <- mxModel(fit_wPRS, name = "noPRS")
model_Null <- omxSetParameters(model_Null, labels = "b7", free = FALSE,
values = 0)
fit_Null <- mxRun(model_Null)
(sum_Null <- summary(fit_Null))## Summary of noPRS
##
## The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. Optimizer was terminated because no further improvement could be made in the merit function (Mx status GREEN).
##
## free parameters:
## name matrix row col Estimate Std.Error A lbound ubound
## 1 b1 bPC1 1 1 596.9028672 444.656189378
## 2 b2 bPC2 1 1 341.5478088 279.361150612
## 3 b3 bPC3 1 1 -60.4220028 141.517064862
## 4 b4 bPC4 1 1 83.9433795 191.642823752
## 5 b5 bsex 1 1 -8.9563185 0.561653963
## 6 b6 bage 1 1 -0.1107354 0.021912748
## 7 meanht meanG 1 1 172.2274408 8.538805730
## 8 var MZ.covMZ ht_1 ht_1 75.2064382 2.879449459 1e-05
## 9 cMZ21 MZ.covMZ ht_1 ht_2 70.3274383 2.892338544 0
## 10 cDZ21 DZ.covDZ ht_1 ht_2 17.6870557 3.240011731 ! 0
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 10 1884 12449.195
## Saturated: NA NA NA
## Independence: NA NA NA
## Number of observations/statistics: 969/1894
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: 8681.1947 12469.195 12469.424
## BIC: -505.6878 12517.957 12486.197
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2022-06-15 00:48:29
## Wall clock time: 2.568562 secs
## optimizer: NPSOL
## OpenMx version number: 2.19.5.1
## Need help? See help(mxSummary)# We can work out the significance of the PRS effect by
# comparing the fit of the models
mxCompare(fit_wPRS, fit_Null)## base comparison ep minus2LL df AIC diffLL diffdf p
## 1 wPRS <NA> 11 12435.322 1883 12457.322 NA NA NA
## 2 wPRS noPRS 10 12449.195 1884 12469.195 13.872356 1 0.00019565525# and calculate r2 as the difference in variance between
# the models gives variance explained by PRS, divided by
# total variance
(fit_Null$MZ$covMZ$values[1, 2] - fit_wPRS$MZ$covMZ$values[1,
2])/fit_Null$MZ$covMZ$values[1, 2]## [1] 0.012362994Calculate GRM using GCTA
./gcta64 --bfile ozbmi --make-grm --out ozGRM## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:48:29 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Options:
##
## --bfile ozbmi
## --make-grm
## --out ozGRM
##
## Note: GRM is computed using the SNPs on the autosomes.
## Reading PLINK FAM file from [ozbmi.fam]...
## 1956 individuals to be included from FAM file.
## 1956 individuals to be included. 651 males, 1305 females, 0 unknown.
## Reading PLINK BIM file from [ozbmi.bim]...
## 9758 SNPs to be included from BIM file(s).
## Computing the genetic relationship matrix (GRM) v2 ...
## Subset 1/1, no. subject 1-1956
## 1956 samples, 9758 markers, 1913946 GRM elements
## IDs for the GRM file have been saved in the file [ozGRM.grm.id]
## Computing GRM...
## 100% finished in 2.1 sec
## 9758 SNPs have been processed.
## Used 9718 valid SNPs.
## The GRM computation is completed.
## Saving GRM...
## GRM has been saved in the file [ozGRM.grm.bin]
## Number of SNPs in each pair of individuals has been saved in the file [ozGRM.grm.N.bin]
##
## Analysis finished at 00:48:31 CEST on Wed Jun 15 2022
## Overall computational time: 2.22 sec.Open GRM and draw some plots
# Open GRM and visualise off-diagonal distribution
ReadGRMBin = function(prefix, AllN = F, size = 4) {
sum_i = function(i) {
return(sum(1:i))
}
BinFileName = paste(prefix, ".grm.bin", sep = "")
NFileName = paste(prefix, ".grm.N.bin", sep = "")
IDFileName = paste(prefix, ".grm.id", sep = "")
id = read.table(IDFileName)
n = dim(id)[1]
BinFile = file(BinFileName, "rb")
grm = readBin(BinFile, n = n * (n + 1)/2, what = numeric(0),
size = size)
NFile = file(NFileName, "rb")
if (AllN == T) {
N = readBin(NFile, n = n * (n + 1)/2, what = numeric(0),
size = size)
} else N = readBin(NFile, n = 1, what = numeric(0), size = size)
i = sapply(1:n, sum_i)
return(list(diag = grm[i], off = grm[-i], id = id, N = N))
}
GRM = ReadGRMBin(prefix = "ozGRM")
hist(GRM$diag, breaks = 100) # Large values can indicate data problems or outstanding homozygosity rate 
#
hist(GRM$off, breaks = 100)
# This plot isn't very informative as most relationship
# coefficents are close to 0. We can update our plot to
# only show relationship coefficents greater than 0.05
# using the following code
hist(GRM$off[which(GRM$off > 0.05)], breaks = 100)
# This plot clearly shows we have a lot of siblings and MZ
# twin pairs. If we zoom in even further we might spot
# those cousins
hist(GRM$off[which(GRM$off > 0.05)], breaks = 100, xlim = c(0.05,
0.3))
Write files for GCTA analysis
prs = read.table("PRSice.all_score", header = T)
covar <- read.table("ozht.cov", header = T)
pheno <- read.table("ozht.pheno", header = T)
# Then rescaling and renaming
pheno$ht <- pheno$ht * 100
names(prs) <- c("FID", "IID", "PRS1", "PRS2", "PRS3", "PRS4",
"PRS5")
# Then merging the data into 1 dataframe
all = merge(covar, prs, by = c("IID", "FID"))
all = merge(all, pheno, by = c("IID", "FID"))
# Then we'll write out the covariate files for the
# continuous covariates that we want to include in the GCTA
# model It's important that the PRS is in the last column
# of this file
write.table(all[, c("FID", "IID", "age", "PC1", "PC2", "PC3",
"PC4", "PRS1")], "ozht.prs1.qcovar", quote = F, row.names = F)
write.table(all[, c("FID", "IID", "age", "PC1", "PC2", "PC3",
"PC4", "PRS2")], "ozht.prs2.qcovar", quote = F, row.names = F)
write.table(all[, c("FID", "IID", "age", "PC1", "PC2", "PC3",
"PC4", "PRS3")], "ozht.prs3.qcovar", quote = F, row.names = F)
write.table(all[, c("FID", "IID", "age", "PC1", "PC2", "PC3",
"PC4", "PRS4")], "ozht.prs4.qcovar", quote = F, row.names = F)
write.table(all[, c("FID", "IID", "age", "PC1", "PC2", "PC3",
"PC4", "PRS5")], "ozht.prs5.qcovar", quote = F, row.names = F)
# Then we'll write out the covariate files for the binary
# covariates that we want to include in the GCTA model
write.table(all[, c("FID", "IID", "sex")], "ozht.covar", quote = F,
row.names = F)
# Then we'll write out the phenotype files for the GCTA
# model
write.table(all[, c("FID", "IID", "ht")], "ozht.phen", quote = F,
row.names = F)Run GCTA with fixed and random effect (AE model)
./gcta64 --pheno ozht.phen --qcovar ozht.prs1.qcovar --covar ozht.covar --grm ozGRM --out ozhtPRS1 --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs2.qcovar --covar ozht.covar --grm ozGRM --out ozhtPRS2 --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs3.qcovar --covar ozht.covar --grm ozGRM --out ozhtPRS3 --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs4.qcovar --covar ozht.covar --grm ozGRM --out ozhtPRS4 --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs5.qcovar --covar ozht.covar --grm ozGRM --out ozhtPRS5 --reml --reml-est-fix
#GCTA will make a .log file and an .hsq file for each model.
#The beta and SE of the PRS effect can be found in the line labeled X_7 in the log files.
#We can pull these out of the log files and save them to a file with the following code
grep X_7 *.log > gcta_results.txt## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:48:33 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs1.qcovar
## --covar ozht.covar
## --grm ozGRM
## --out ozhtPRS1
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs1.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs1.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 2 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 43.9528 34.2896
## logL: -4810.94
## Running AI-REML algorithm ...
## Iter. logL V(G) V(e)
## 1 -4751.85 59.84011 17.40093
## 2 -4647.03 70.13324 10.91462
## 3 -4590.77 76.39449 8.49173
## 4 -4567.32 80.61353 7.36394
## 5 -4557.13 83.63826 6.74011
## 6 -4552.33 85.84919 6.35988
## 7 -4549.96 87.47110 6.11437
## 8 -4548.77 88.65891 5.94998
## 9 -4548.16 89.52614 5.83726
## 10 -4547.85 91.52381 5.58866
## 11 -4547.51 91.78332 5.56727
## 12 -4547.51 91.80849 5.56429
## 13 -4547.51 91.81150 5.56392
## Log-likelihood ratio converged.
##
## Calculating the logLikelihood for the reduced model ...
## (variance component 1 is dropped from the model)
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 69.78139
## logL: -4980.18681
## Running AI-REML algorithm ...
## Iter. logL V(e)
## 1 -4948.05 69.75216
## 2 -4948.05 69.73219
## 3 -4948.05 69.68899
## 4 -4948.05 69.68902
## 5 -4948.05 69.68902
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G) 91.811500 3.895050
## V(e) 5.563916 0.357498
## Vp 97.375416 3.852776
## V(G)/Vp 0.942861 0.004449
##
## Sampling variance/covariance of the estimates of variance components:
## 1.517142e+01 -2.276706e-01
## -2.276706e-01 1.278046e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 169.095957 8.764314
## X_2 -0.113977 0.022456
## X_3 -17.551326 484.849288
## X_4 410.861460 292.698935
## X_5 51.660560 189.243241
## X_6 -41.123137 179.306181
## X_7 3.182480 0.808135
## X_8 -9.535405 0.566582
##
## Summary result of REML analysis has been saved in the file [ozhtPRS1.hsq].
##
## Analysis finished at 00:48:44 CEST on Wed Jun 15 2022
## Overall computational time: 11.16 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:48:44 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs2.qcovar
## --covar ozht.covar
## --grm ozGRM
## --out ozhtPRS2
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs2.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs2.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 2 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 43.828 34.2361
## logL: -4807.36
## Running AI-REML algorithm ...
## Iter. logL V(G) V(e)
## 1 -4747.62 59.53572 17.48921
## 2 -4644.39 69.76967 10.98529
## 3 -4588.49 76.01234 8.53453
## 4 -4564.96 80.21423 7.39286
## 5 -4554.71 83.22495 6.76154
## 6 -4549.88 85.42585 6.37683
## 7 -4547.49 87.04118 6.12843
## 8 -4546.29 88.22501 5.96209
## 9 -4545.67 89.09008 5.84799
## 10 -4545.36 91.08456 5.59628
## 11 -4545.02 91.34835 5.57435
## 12 -4545.01 91.37438 5.57124
## 13 -4545.01 91.37754 5.57085
## Log-likelihood ratio converged.
##
## Calculating the logLikelihood for the reduced model ...
## (variance component 1 is dropped from the model)
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 69.17535
## logL: -4974.01940
## Running AI-REML algorithm ...
## Iter. logL V(e)
## 1 -4939.88 69.14531
## 2 -4939.88 69.12478
## 3 -4939.88 69.08038
## 4 -4939.88 69.08041
## 5 -4939.88 69.08041
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G) 91.377539 3.879798
## V(e) 5.570853 0.358140
## Vp 96.948392 3.837078
## V(G)/Vp 0.942538 0.004478
##
## Sampling variance/covariance of the estimates of variance components:
## 1.505283e+01 -2.289649e-01
## -2.289649e-01 1.282644e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.417530 8.736724
## X_2 -0.112660 0.022406
## X_3 15.076693 483.651444
## X_4 406.821756 292.025107
## X_5 54.259660 188.838053
## X_6 -65.133325 178.945021
## X_7 3.286590 0.720470
## X_8 -9.516487 0.565417
##
## Summary result of REML analysis has been saved in the file [ozhtPRS2.hsq].
##
## Analysis finished at 00:48:55 CEST on Wed Jun 15 2022
## Overall computational time: 11.14 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:48:55 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs3.qcovar
## --covar ozht.covar
## --grm ozGRM
## --out ozhtPRS3
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs3.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs3.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 2 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 43.7761 34.2641
## logL: -4806.95
## Running AI-REML algorithm ...
## Iter. logL V(G) V(e)
## 1 -4747.50 59.39239 17.64269
## 2 -4646.02 69.67637 11.07406
## 3 -4590.07 75.97320 8.58117
## 4 -4566.24 80.20843 7.42150
## 5 -4555.84 83.24244 6.78120
## 6 -4550.93 85.46101 6.39135
## 7 -4548.50 87.09011 6.13976
## 8 -4547.28 88.28472 5.97131
## 9 -4546.65 89.15818 5.85577
## 10 -4546.33 91.17320 5.60084
## 11 -4545.98 91.44246 5.57847
## 12 -4545.97 91.46926 5.57527
## 13 -4545.97 91.47254 5.57487
## Log-likelihood ratio converged.
##
## Calculating the logLikelihood for the reduced model ...
## (variance component 1 is dropped from the model)
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 68.88252
## logL: -4971.07558
## Running AI-REML algorithm ...
## Iter. logL V(e)
## 1 -4935.95 68.85209
## 2 -4935.95 68.83129
## 3 -4935.94 68.78631
## 4 -4935.94 68.78634
## 5 -4935.94 68.78634
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G) 91.472544 3.884525
## V(e) 5.574867 0.358533
## Vp 97.047411 3.841601
## V(G)/Vp 0.942555 0.004479
##
## Sampling variance/covariance of the estimates of variance components:
## 1.508953e+01 -2.300927e-01
## -2.300927e-01 1.285462e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.639132 8.743245
## X_2 -0.112417 0.022418
## X_3 -16.914727 484.031778
## X_4 417.248000 292.245842
## X_5 47.967643 188.946734
## X_6 -54.884823 179.008572
## X_7 2.898535 0.664494
## X_8 -9.535677 0.565604
##
## Summary result of REML analysis has been saved in the file [ozhtPRS3.hsq].
##
## Analysis finished at 00:49:06 CEST on Wed Jun 15 2022
## Overall computational time: 10.93 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:49:06 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs4.qcovar
## --covar ozht.covar
## --grm ozGRM
## --out ozhtPRS4
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs4.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs4.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 2 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 43.7703 34.2567
## logL: -4806.73
## Running AI-REML algorithm ...
## Iter. logL V(G) V(e)
## 1 -4747.20 59.39007 17.62675
## 2 -4645.59 69.66570 11.06386
## 3 -4589.65 75.95222 8.57561
## 4 -4565.86 80.17987 7.41804
## 5 -4555.48 83.20819 6.77879
## 6 -4550.58 85.42239 6.38956
## 7 -4548.16 87.04814 6.13836
## 8 -4546.94 88.24021 5.97017
## 9 -4546.31 89.11175 5.85481
## 10 -4545.99 91.12224 5.60027
## 11 -4545.65 91.39066 5.57795
## 12 -4545.64 91.41736 5.57477
## 13 -4545.64 91.42062 5.57436
## Log-likelihood ratio converged.
##
## Calculating the logLikelihood for the reduced model ...
## (variance component 1 is dropped from the model)
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 68.84960
## logL: -4970.79105
## Running AI-REML algorithm ...
## Iter. logL V(e)
## 1 -4935.55 68.81912
## 2 -4935.55 68.79829
## 3 -4935.55 68.75324
## 4 -4935.55 68.75327
## 5 -4935.55 68.75327
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G) 91.420622 3.882391
## V(e) 5.574364 0.358485
## Vp 96.994986 3.839481
## V(G)/Vp 0.942529 0.004480
##
## Sampling variance/covariance of the estimates of variance components:
## 1.507296e+01 -2.299303e-01
## -2.299303e-01 1.285118e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.212392 8.737739
## X_2 -0.112576 0.022412
## X_3 2.572887 483.807226
## X_4 412.614296 292.132395
## X_5 48.578056 188.894702
## X_6 -57.183480 178.965315
## X_7 2.801184 0.629326
## X_8 -9.528049 0.565490
##
## Summary result of REML analysis has been saved in the file [ozhtPRS4.hsq].
##
## Analysis finished at 00:49:17 CEST on Wed Jun 15 2022
## Overall computational time: 10.97 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:49:17 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs5.qcovar
## --covar ozht.covar
## --grm ozGRM
## --out ozhtPRS5
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs5.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs5.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 2 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 43.7455 34.2568
## logL: -4806.24
## Running AI-REML algorithm ...
## Iter. logL V(G) V(e)
## 1 -4746.69 59.33900 17.65774
## 2 -4645.51 69.62249 11.07936
## 3 -4589.55 75.91415 8.58302
## 4 -4565.71 80.14304 7.42229
## 5 -4555.30 83.17158 6.78155
## 6 -4550.39 85.38582 6.39149
## 7 -4547.97 87.01165 6.13978
## 8 -4546.74 88.20385 5.97127
## 9 -4546.12 89.07556 5.85569
## 10 -4545.79 91.08660 5.60067
## 11 -4545.45 91.35558 5.57830
## 12 -4545.44 91.38236 5.57510
## 13 -4545.44 91.38565 5.57470
## Log-likelihood ratio converged.
##
## Calculating the logLikelihood for the reduced model ...
## (variance component 1 is dropped from the model)
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 68.71671
## logL: -4969.43466
## Running AI-REML algorithm ...
## Iter. logL V(e)
## 1 -4933.73 68.68605
## 2 -4933.73 68.66510
## 3 -4933.73 68.61979
## 4 -4933.73 68.61982
## 5 -4933.73 68.61982
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G) 91.385646 3.881152
## V(e) 5.574698 0.358524
## Vp 96.960344 3.838207
## V(G)/Vp 0.942505 0.004483
##
## Sampling variance/covariance of the estimates of variance components:
## 1.506334e+01 -2.300251e-01
## -2.300251e-01 1.285397e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.446708 8.737619
## X_2 -0.113225 0.022407
## X_3 -9.622066 483.773546
## X_4 412.972238 292.081569
## X_5 48.675267 188.861598
## X_6 -57.801962 178.935760
## X_7 2.763790 0.613928
## X_8 -9.544225 0.565300
##
## Summary result of REML analysis has been saved in the file [ozhtPRS5.hsq].
##
## Analysis finished at 00:49:28 CEST on Wed Jun 15 2022
## Overall computational time: 10.91 sec.Format GCTA outputs predictivness
# We can read in the gcta results
res <- as.data.frame(read.table("gcta_results.txt"))
names(res) <- c("Model", "beta", "se")
# To calculate R2 we scale beta into a marginal correlation
# and then square it We are calculating this and adding it
# to the res object
res$r2 <- rbind((res[1, 2]/sd(all$ht, na.rm = T) * sd(all$PRS1,
na.rm = T))^2, (res[2, 2]/sd(all$ht, na.rm = T) * sd(all$PRS2,
na.rm = T))^2, (res[3, 2]/sd(all$ht, na.rm = T) * sd(all$PRS3,
na.rm = T))^2, (res[4, 2]/sd(all$ht, na.rm = T) * sd(all$PRS4,
na.rm = T))^2, (res[5, 2]/sd(all$ht, na.rm = T) * sd(all$PRS5,
na.rm = T))^2)
# To test the significance of the betas we can use t-tests
# We are calculating this and adding it to the res object
res$tStat <- res$beta/res$se # Follows a student distribution
res$pval <- (1 - pt(q = res$tStat, df = 1887, lower.tail = T)) *
2
# Now we can look at the results
resMake barplots to summarise associations
r2 <- data.frame(name = c("5e-08", "1e-06", "0.0001", "0.01",
"1"), R2 = c(1.236, 1.913, 2.218, 2.34, 2.561), P = c(0.000195,
5.04e-06, 6.91e-07, 2.67e-07, 6.12e-08))
# Increase bottom margin
par(mar = c(4, 4, 4, 4))
# Barplot jpeg('height.R2.jpg', width = 800, height = 800,
# quality=75, bg = 'white', type = 'cairo')
my_bar <- barplot(r2$R2, border = F, names.arg = r2$name, las = 2,
col = c("skyblue1", "skyblue2", "skyblue3", "skyblue4"),
ylim = c(0, 3), main = "Predicting height from height PRS",
ylab = "% variance explained", xlab = "PRS Thresholds")
# Add P values to the bars
text(my_bar, r2$R2 + 0.04, paste("P=", r2$P, sep = ""), cex = 1)
# dev.off()
# Next we'll plot height by quintile of PRS (often these
# plots are made with deciles but the sample is small so
# we're going to use quintiles here)
library(Hmisc)
library(gplots)
# the following command divides the data into quintiles.
longData$percentiles <- factor(cut2(longData$PRS5, g = 5, levels.mean = TRUE))
# Plot the mean of height by precentile with 95% CIs
# jpeg('Quintile.jpg', width = 800, height = 800,
# quality=75, bg = 'white', type = 'cairo')
plotmeans(ht ~ percentiles, data = longData, xlab = "Quintile of PRS",
ylab = "Height in Cm", legends = c("1", "2", "3", "4", "5"))
# dev.off()Optional analysis - adapt code to fit ACE model in GCTA
Make matrix of variance covariance for the shared environment
# We can make the CRM by modyfing the GRM
hist(GRM$off[which(GRM$off > 0.2)], breaks = 100) # Remember that the sample contains MZ and DZ twins (or siblings)
# Format GRM into a NxN matrix (simpler to write later)
asNxNGRM <- function(BRMbin) {
mat <- matrix(0, nrow = length(BRMbin$diag), ncol = length(BRMbin$diag))
mat[upper.tri(mat)] <- BRMbin$off
mat <- mat + t(mat)
diag(mat) <- BRMbin$diag
colnames(mat) <- BRMbin$id[, 2]
rownames(mat) <- BRMbin$id[, 2]
return(mat)
}
# convert GRM into NxN matrix
GRMmat = asNxNGRM(GRM)
# We can build the CRM as a modified GRM
CRM = GRMmat
CRM[which(CRM > 0.2, arr.ind = T)] = 1 # we assume that individuals with GRM>0.2 have grown up together -> C=1
CRM[which(CRM < 0.2, arr.ind = T)] = 0 # other less related individuals have no shared environment -> C=0
# Check how many 1 and 0 we have
table(CRM[upper.tri(CRM)]) # 978 families##
## 0 1
## 1911012 978table(diag(CRM)) # check diagonal##
## 1
## 1956# We get the info about family ID and IID - needed to
# create the .grm.id file
fam = GRM$id
colnames(fam) = c("fam", "id")
# We get this package which allows writing a GRM
# install.packages('genio')
library(genio) # needed to write_grm
genio::write_grm(kinship = CRM, name = "ozCRM", fam = fam)Fit several random effects in GCTA
# We create this text file, which stores the name of the different GRM/CRM matrices so GCTA knows which to open
echo "ozGRM" > GRM-CRM.txt
echo "ozCRM" >> GRM-CRM.txt
head GRM-CRM.txt
./gcta64 --pheno ozht.phen --qcovar ozht.prs1.qcovar --covar ozht.covar --mgrm GRM-CRM.txt --out ozhtPRS1.ACE --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs2.qcovar --covar ozht.covar --mgrm GRM-CRM.txt --out ozhtPRS2.ACE --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs3.qcovar --covar ozht.covar --mgrm GRM-CRM.txt --out ozhtPRS3.ACE --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs4.qcovar --covar ozht.covar --mgrm GRM-CRM.txt --out ozhtPRS4.ACE --reml --reml-est-fix
./gcta64 --pheno ozht.phen --qcovar ozht.prs5.qcovar --covar ozht.covar --mgrm GRM-CRM.txt --out ozhtPRS5.ACE --reml --reml-est-fix
grep X_7 *ACE.log > gcta_results_ACE.txt## ozGRM
## ozCRM
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:49:34 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs1.qcovar
## --covar ozht.covar
## --mgrm GRM-CRM.txt
## --out ozhtPRS1.ACE
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs1.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs1.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## There are 2 GRM file names specified in the file [GRM-CRM.txt].
## Reading the GRM from the 1th file ...
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading the GRM from the 2th file ...
## Reading IDs of the GRM from [ozCRM.grm.id].
## 1956 IDs are read from [ozCRM.grm.id].
## Reading the GRM from [ozCRM.grm.bin].
## GRM for 1956 individuals are included from [ozCRM.grm.bin].
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 3 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 29.9405 29.1033 24.645
## logL: -4736.31
## Running AI-REML algorithm ...
## Iter. logL V(G1) V(G2) V(e)
## 1 -4700.32 50.36555 18.19224 13.55946
## 2 -4611.85 61.35311 13.83297 9.51921
## 3 -4573.17 67.87119 11.53647 7.88367
## 4 -4557.75 72.23652 10.14153 7.04864
## 5 -4550.93 75.29060 9.25897 6.56151
## 6 -4547.70 77.45579 8.69594 6.25505
## 7 -4546.10 78.99608 8.33749 6.05316
## 8 -4545.30 80.09254 8.11040 5.91617
## 9 -4544.89 82.56309 7.65749 5.61614
## 10 -4544.46 82.77717 7.74327 5.59209
## 11 -4544.46 82.82704 7.72411 5.58854
## 12 -4544.46 82.82629 7.72855 5.58805
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G1) 82.826286 5.203826
## V(G2) 7.728551 3.584645
## V(e) 5.588054 0.359747
## Vp 96.142891 3.815144
## V(G1)/Vp 0.861492 0.037840
## V(G2)/Vp 0.080386 0.037338
##
## Sum of V(G)/Vp 0.941878 0.004543
##
## Sampling variance/covariance of the estimates of variance components:
## 2.707981e+01 -1.251748e+01 -2.574890e-01
## -1.251748e+01 1.284968e+01 2.317447e-02
## -2.574890e-01 2.317447e-02 1.294179e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 169.255966 8.892221
## X_2 -0.112214 0.023057
## X_3 -49.312966 489.526725
## X_4 414.675175 296.962542
## X_5 39.556627 191.180832
## X_6 -31.085455 181.347441
## X_7 3.071223 0.804539
## X_8 -9.529359 0.582336
##
## Summary result of REML analysis has been saved in the file [ozhtPRS1.ACE.hsq].
##
## Analysis finished at 00:49:43 CEST on Wed Jun 15 2022
## Overall computational time: 9.18 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:49:43 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs2.qcovar
## --covar ozht.covar
## --mgrm GRM-CRM.txt
## --out ozhtPRS2.ACE
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs2.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs2.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## There are 2 GRM file names specified in the file [GRM-CRM.txt].
## Reading the GRM from the 1th file ...
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading the GRM from the 2th file ...
## Reading IDs of the GRM from [ozCRM.grm.id].
## 1956 IDs are read from [ozCRM.grm.id].
## Reading the GRM from [ozCRM.grm.bin].
## GRM for 1956 individuals are included from [ozCRM.grm.bin].
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 3 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 29.8995 29.0562 24.6377
## logL: -4733.47
## Running AI-REML algorithm ...
## Iter. logL V(G1) V(G2) V(e)
## 1 -4697.25 50.24276 18.04696 13.61573
## 2 -4609.52 61.21929 13.64557 9.56059
## 3 -4570.92 67.72897 11.33108 7.91190
## 4 -4555.47 72.08346 9.92817 7.06959
## 5 -4548.62 75.12825 9.04166 6.57814
## 6 -4545.34 77.28681 8.47639 6.26889
## 7 -4543.77 78.82292 8.11646 6.06511
## 8 -4542.96 79.91706 7.88831 5.92680
## 9 -4542.55 82.38431 7.43282 5.62374
## 10 -4542.11 82.60353 7.51772 5.59911
## 11 -4542.11 82.65447 7.49848 5.59543
## 12 -4542.11 82.65387 7.50293 5.59492
## 13 -4542.11 82.65550 7.50186 5.59486
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G1) 82.655495 5.190102
## V(G2) 7.501856 3.565254
## V(e) 5.594857 0.360340
## Vp 95.752208 3.801264
## V(G1)/Vp 0.863223 0.037804
## V(G2)/Vp 0.078347 0.037294
##
## Sum of V(G)/Vp 0.941569 0.004571
##
## Sampling variance/covariance of the estimates of variance components:
## 2.693716e+01 -1.242863e+01 -2.588074e-01
## -1.242863e+01 1.271104e+01 2.321948e-02
## -2.588074e-01 2.321948e-02 1.298450e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.624485 8.862613
## X_2 -0.110984 0.022994
## X_3 -17.724628 488.242119
## X_4 409.842478 296.195656
## X_5 43.235878 190.742260
## X_6 -55.042363 180.962099
## X_7 3.174626 0.717059
## X_8 -9.508999 0.580801
##
## Summary result of REML analysis has been saved in the file [ozhtPRS2.ACE.hsq].
##
## Analysis finished at 00:49:53 CEST on Wed Jun 15 2022
## Overall computational time: 10.07 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:49:53 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs3.qcovar
## --covar ozht.covar
## --mgrm GRM-CRM.txt
## --out ozhtPRS3.ACE
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs3.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs3.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## There are 2 GRM file names specified in the file [GRM-CRM.txt].
## Reading the GRM from the 1th file ...
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading the GRM from the 2th file ...
## Reading IDs of the GRM from [ozCRM.grm.id].
## 1956 IDs are read from [ozCRM.grm.id].
## Reading the GRM from [ozCRM.grm.bin].
## GRM for 1956 individuals are included from [ozCRM.grm.bin].
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 3 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 29.8736 29.0475 24.6513
## logL: -4732.9
## Running AI-REML algorithm ...
## Iter. logL V(G1) V(G2) V(e)
## 1 -4696.82 50.02688 18.17504 13.71829
## 2 -4610.66 61.02046 13.79117 9.62093
## 3 -4572.02 67.54188 11.49160 7.94854
## 4 -4556.43 71.90039 10.10100 7.09463
## 5 -4549.51 74.94802 9.22342 6.59669
## 6 -4546.23 77.10985 8.66428 6.28347
## 7 -4544.60 78.64964 8.30847 6.07709
## 8 -4543.78 79.74753 8.08303 5.93699
## 9 -4543.36 82.22597 7.63312 5.62996
## 10 -4542.92 82.45294 7.71726 5.60473
## 11 -4542.91 82.50486 7.69811 5.60091
## 12 -4542.91 82.50445 7.70257 5.60038
## 13 -4542.91 82.50611 7.70149 5.60032
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G1) 82.506111 5.184994
## V(G2) 7.701490 3.567178
## V(e) 5.600317 0.360871
## Vp 95.807918 3.804076
## V(G1)/Vp 0.861162 0.037797
## V(G2)/Vp 0.080385 0.037287
##
## Sum of V(G)/Vp 0.941546 0.004575
##
## Sampling variance/covariance of the estimates of variance components:
## 2.688416e+01 -1.239702e+01 -2.608549e-01
## -1.239702e+01 1.272476e+01 2.379687e-02
## -2.608549e-01 2.379687e-02 1.302279e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.920867 8.871820
## X_2 -0.110746 0.023018
## X_3 -52.818537 488.707817
## X_4 419.092591 296.477111
## X_5 36.334022 190.873399
## X_6 -44.340405 181.048200
## X_7 2.813803 0.661036
## X_8 -9.527184 0.581329
##
## Summary result of REML analysis has been saved in the file [ozhtPRS3.ACE.hsq].
##
## Analysis finished at 00:50:03 CEST on Wed Jun 15 2022
## Overall computational time: 9.54 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:50:03 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs4.qcovar
## --covar ozht.covar
## --mgrm GRM-CRM.txt
## --out ozhtPRS4.ACE
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs4.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs4.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## There are 2 GRM file names specified in the file [GRM-CRM.txt].
## Reading the GRM from the 1th file ...
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading the GRM from the 2th file ...
## Reading IDs of the GRM from [ozCRM.grm.id].
## 1956 IDs are read from [ozCRM.grm.id].
## Reading the GRM from [ozCRM.grm.bin].
## GRM for 1956 individuals are included from [ozCRM.grm.bin].
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 3 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 29.8714 29.0453 24.6475
## logL: -4732.7
## Running AI-REML algorithm ...
## Iter. logL V(G1) V(G2) V(e)
## 1 -4696.57 50.05982 18.14688 13.69875
## 2 -4610.19 61.03608 13.76834 9.61007
## 3 -4571.60 67.54205 11.47352 7.94229
## 4 -4556.05 71.88905 10.08655 7.09055
## 5 -4549.15 74.92784 9.21178 6.59380
## 6 -4545.87 77.08286 8.65476 6.28129
## 7 -4544.25 78.61747 8.30052 6.07537
## 8 -4543.44 79.71148 8.07621 5.93558
## 9 -4543.02 82.18079 7.62885 5.62923
## 10 -4542.58 82.40633 7.71308 5.60410
## 11 -4542.57 82.45812 7.69390 5.60030
## 12 -4542.57 82.45766 7.69839 5.59977
## 13 -4542.57 82.45932 7.69729 5.59971
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G1) 82.459325 5.181032
## V(G2) 7.697294 3.563947
## V(e) 5.599706 0.360812
## Vp 95.756325 3.801953
## V(G1)/Vp 0.861137 0.037783
## V(G2)/Vp 0.080384 0.037273
##
## Sum of V(G)/Vp 0.941521 0.004577
##
## Sampling variance/covariance of the estimates of variance components:
## 2.684309e+01 -1.237322e+01 -2.605583e-01
## -1.237322e+01 1.270172e+01 2.370034e-02
## -2.605583e-01 2.370034e-02 1.301853e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.476664 8.865994
## X_2 -0.110867 0.023012
## X_3 -35.580668 488.477468
## X_4 415.705484 296.370772
## X_5 36.669484 190.821852
## X_6 -47.537680 181.007689
## X_7 2.723517 0.626279
## X_8 -9.520307 0.581206
##
## Summary result of REML analysis has been saved in the file [ozhtPRS4.ACE.hsq].
##
## Analysis finished at 00:50:12 CEST on Wed Jun 15 2022
## Overall computational time: 8.75 sec.
## *******************************************************************
## * Genome-wide Complex Trait Analysis (GCTA)
## * version 1.94.0 beta Mac
## * (C) 2010-present, Jian Yang, The University of Queensland
## * Please report bugs to Jian Yang <jian.yang.qt@gmail.com>
## *******************************************************************
## Analysis started at 00:50:12 CEST on Wed Jun 15 2022.
## Hostname: imb17-013265-lt
##
## Accepted options:
## --pheno ozht.phen
## --qcovar ozht.prs5.qcovar
## --covar ozht.covar
## --mgrm GRM-CRM.txt
## --out ozhtPRS5.ACE
## --reml
## --reml-est-fix
##
## Note: This is a multi-thread program. You could specify the number of threads by the --thread-num option to speed up the computation if there are multiple processors in your machine.
##
## Reading phenotypes from [ozht.phen].
## Non-missing phenotypes of 1895 individuals are included from [ozht.phen].
## Reading quantitative covariate(s) from [ozht.prs5.qcovar].
## 6 quantitative covariate(s) of 1933 individuals are included from [ozht.prs5.qcovar].
## Reading discrete covariate(s) from [ozht.covar].
## 1 discrete covariate(s) of 1933 individuals are included from [ozht.covar].
##
## There are 2 GRM file names specified in the file [GRM-CRM.txt].
## Reading the GRM from the 1th file ...
## Reading IDs of the GRM from [ozGRM.grm.id].
## 1956 IDs are read from [ozGRM.grm.id].
## Reading the GRM from [ozGRM.grm.bin].
## GRM for 1956 individuals are included from [ozGRM.grm.bin].
## Reading the GRM from the 2th file ...
## Reading IDs of the GRM from [ozCRM.grm.id].
## 1956 IDs are read from [ozCRM.grm.id].
## Reading the GRM from [ozCRM.grm.bin].
## GRM for 1956 individuals are included from [ozCRM.grm.bin].
## 6 quantitative variable(s) included as covariate(s).
## 1 discrete variable(s) included as covariate(s).
## 1894 individuals are in common in these files.
##
## Performing REML analysis ... (Note: may take hours depending on sample size).
## 1894 observations, 8 fixed effect(s), and 3 variance component(s)(including residual variance).
## Calculating prior values of variance components by EM-REML ...
## Updated prior values: 29.8578 29.0401 24.6464
## logL: -4732.1
## Running AI-REML algorithm ...
## Iter. logL V(G1) V(G2) V(e)
## 1 -4695.94 49.99573 18.17732 13.71870
## 2 -4609.97 60.96663 13.80366 9.62128
## 3 -4571.39 67.46612 11.51264 7.94896
## 4 -4555.83 71.80684 10.12860 7.09504
## 5 -4548.92 74.84089 9.25582 6.59709
## 6 -4545.63 76.99269 8.70006 6.28385
## 7 -4544.02 78.52530 8.34655 6.07745
## 8 -4543.20 79.61813 8.12264 5.93734
## 9 -4542.78 82.08548 7.67591 5.63026
## 10 -4542.34 82.31263 7.75955 5.60502
## 11 -4542.33 82.36451 7.74045 5.60120
## 12 -4542.33 82.36411 7.74491 5.60067
## 13 -4542.33 82.36578 7.74383 5.60060
## Log-likelihood ratio converged.
##
## Summary result of REML analysis:
## Source Variance SE
## V(G1) 82.365776 5.177148
## V(G2) 7.743826 3.562831
## V(e) 5.600599 0.360904
## Vp 95.710201 3.800263
## V(G1)/Vp 0.860575 0.037789
## V(G2)/Vp 0.080909 0.037278
##
## Sum of V(G)/Vp 0.941484 0.004581
##
## Sampling variance/covariance of the estimates of variance components:
## 2.680287e+01 -1.235535e+01 -2.610650e-01
## -1.235535e+01 1.269376e+01 2.397526e-02
## -2.610650e-01 2.397526e-02 1.302519e-01
## Estimatesof fixed effects:
##
## Source Estimate SE
## mean 168.700024 8.866271
## X_2 -0.111501 0.023010
## X_3 -48.135784 488.464289
## X_4 416.564315 296.337132
## X_5 36.881111 190.792892
## X_6 -48.392397 180.984907
## X_7 2.694306 0.611060
## X_8 -9.536273 0.581096
##
## Summary result of REML analysis has been saved in the file [ozhtPRS5.ACE.hsq].
##
## Analysis finished at 00:50:21 CEST on Wed Jun 15 2022
## Overall computational time: 9.62 sec.Format results
# We can read in the gcta results
res <- as.data.frame(read.table("gcta_results_ACE.txt"))
names(res) <- c("Model", "beta", "se")
# To calculate R2 we scale beta into a marginal correlation
# and then square it We are calculating this and adding it
# to the res object
res$r2 <- rbind((res[1, 2]/sd(all$ht, na.rm = T) * sd(all$PRS1,
na.rm = T))^2, (res[2, 2]/sd(all$ht, na.rm = T) * sd(all$PRS2,
na.rm = T))^2, (res[3, 2]/sd(all$ht, na.rm = T) * sd(all$PRS3,
na.rm = T))^2, (res[4, 2]/sd(all$ht, na.rm = T) * sd(all$PRS4,
na.rm = T))^2, (res[5, 2]/sd(all$ht, na.rm = T) * sd(all$PRS5,
na.rm = T))^2)
# To test the significance of the betas we can use t-tests
# We are calculating this and adding it to the res object
res$tStat <- res$beta/res$se # Follows a student distribution
res$pval <- (1 - pt(q = res$tStat, df = 1887, lower.tail = T)) *
2
# Now we can look at the results
resA work by Sarah Medland & Baptiste Couvy-Duchesne