Last active
June 29, 2021 21:12
-
-
Save nievergeltlab/fd79d545973f159b53cb7d25c16d60f8 to your computer and use it in GitHub Desktop.
User R to analyze a single SNP from PLINK. The SNP is analyzed using an ordered logit model, with the SNP as an outcome (usually it's done as OLS or logistic with the SNP as a predictor).
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#Test the association between a SNP and PTSD, in R, assuming ordered logit of the SNP | |
#Do the following in PLINK | |
plink --bfile --snp rs4680 --recodeA --out rs4680 | |
#Open the corresponding .raw file and note the column name for the SNP. you'll need this for R | |
#Then load R | |
R | |
library(MASS) | |
setwd('C:/Users/adam/Desktop/PGC-PTSD/phenotype_harmonization/shareefa_check/') | |
#Merge genotype, phenotype, and covariates | |
#Assuming that each of the sheets has columns FID and IID which identify subjects | |
gene <- read.table('rs4680.raw', header=T) | |
pheno <- read.table('C:/Users/adam/Desktop/PGC-PTSD/phenotype_harmonization/shareefa_check/Copy of C7Y6N1S2_6_2_15_COVs_PCAs_withSA_SD_harm_Jan20162.csv', sep=";", header=T,na.strings=c("NA","#N/A","-9")) | |
dat <- merge(gene,pheno,all=TRUE,by=c("FID","IID")) | |
m <- polr(as.factor(rs4680_G) ~ as.numeric(PTSDCurrentContinuous), data = dat, Hess=TRUE) | |
(ctable <- coef(summary(m))) | |
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2 | |
#View beta, se, test stat, P value | |
(ctable <- cbind(ctable, "p value" = p)) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment