## Produce prediction R2 by comparing the full model with PRS to the null model without PRS

args = commandArgs(trailingOnly=TRUE)

phenFile = args[1]
covFile = args[2]
prsFile = args[3]

# Read in the phenotype file 
phenotype = read.table(phenFile, header=F)
colnames(phenotype) <- c("FID", "IID", "Trait") 
# Read in the covariates
covariates = read.table(covFile, header=F)
colnames(covariates)[1:2] = c("FID", "IID")
# Now merge the files
pheno <- merge(phenotype, covariates, by=c("FID", "IID"))
# We can then calculate the null model (model with PRS) using a linear regression 
# (as height is quantitative)
null.model <- lm(Trait~., data=pheno)
# And the R2 of the null model is 
null.r2 <- summary(null.model)$r.squared
# Calculate PRS based on the selected P-value threshold obtained from the target sample
prs <- read.table(prsFile, header=T)
# relevant columns
pheno.prs <- merge(pheno, prs[,c("FID","IID", "SCORESUM")], by=c("FID", "IID"))
# Now perform a linear regression on Height with PRS and the covariates
# ignoring the FID and IID from our model
model <- lm(Trait~., data=pheno.prs)
# model R2 is obtained as 
model.r2 <- summary(model)$r.squared
# R2 of PRS is simply calculated as the model R2 minus the null R2
cat("Prediction R2: ", model.r2-null.r2, "\n")