#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))