nievergeltlab
2/23/2017 - 10:38 PM

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 a

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

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