Skip to content

Instantly share code, notes, and snippets.

@nievergeltlab
Last active June 29, 2021 21:12
Show Gist options
  • Save nievergeltlab/fd79d545973f159b53cb7d25c16d60f8 to your computer and use it in GitHub Desktop.
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).
#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