1
0
Entering edit mode
8 weeks ago

Hi everyone,

Some questions: this is the code on the PLINK tutorial to calculate PRS. https://choishingwan.github.io/PRS-Tutorial/plink/ (go down). This is done in R-studio. The phenotype data I'm working with is not continuous but binary. So 1 stands for control and 2 stands for cases. However here, if I'm understanding correctly, they are working with continuous phenotype. In this script they are doing the linear regression, but since my phenotype is control/case it's binary right. And that means I should be working with logistic regression. However, can someone help me to start?

  library(data.table)
library(magrittr)
p.threshold <- c(0.001,0.05,0.1,0.2,0.3,0.4,0.5)
setnames(., colnames(.), c("FID", "IID", paste0("PC",1:6)) )
pheno <- merge(phenotype, covariate) %>%
merge(., pcs)
null.r2 <- summary(lm(Height~., data=pheno[,-c("FID", "IID")]))$r.squared prs.result <- NULL for(i in p.threshold){ pheno.prs <- paste0("EUR.", i, ".profile") %>% fread(.) %>% .[,c("FID", "IID", "SCORE")] %>% merge(., pheno, by=c("FID", "IID")) model <- lm(Height~., data=pheno.prs[,-c("FID","IID")]) %>% summary model.r2 <- model$r.squared
prs.r2 <- model.r2-null.r2
prs.coef <- model$coeff["SCORE",] prs.result %<>% rbind(., data.frame(Threshold=i, R2=prs.r2, P=as.numeric(prs.coef[4]), BETA=as.numeric(prs.coef[1]), SE=as.numeric(prs.coef[2]))) } print(prs.result[which.max(prs.result$R2),])
q() # exit R

regression plink PRS logistic • 264 views
0
Entering edit mode
7 weeks ago
Sam ★ 4.6k

replace lm with glm(formula, data, family=binomial) (replace formula with your regression formula, and data with your data frame)

Make sure your case control is coded as 0 (control) and 1 (case)

to get the R2, use fmsb::NagelkerkeR2(glm.model)

where glm.model is the glm object.

If these are too complicated, I'd recommend you to use dedicated software such as PRSice-2 to get your PRS results.