From b0a2f65c746740bb112345841e4ad1d01c2f5920 Mon Sep 17 00:00:00 2001
From: "M.E. Ketelaar" <m.e.ketelaar@rug.nl>
Date: Mon, 17 Feb 2020 19:37:56 +0000
Subject: [PATCH]

---
 eQTL_Scripts.md | 53 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 53 insertions(+)
 create mode 100644 eQTL_Scripts.md

diff --git a/eQTL_Scripts.md b/eQTL_Scripts.md
new file mode 100644
index 0000000..bbc293a
--- /dev/null
+++ b/eQTL_Scripts.md
@@ -0,0 +1,53 @@
+#This script belongs to the paper of Ketelaar, Portelli, Dijk et al. 
+#Entitled 'Phenotypic and functional translation of IL33 genetics in asthma', submitted to JACI
+
+#Below example scripts performing eQTL analyses using R statistics
+#Core Team R. R: A language and environment for statistical computing. 2013; Available at: http://www.R-project.org/.
+
+# Due to ethical restrictions we can only share code and summary statistics, no raw input data files
+
+## set working directory
+setwd("H:/x/x/x")   
+
+## Load snpdata         
+source('./snpdata.load.R') 
+  
+
+## Load and combine covariate and phenotype data:
+data<- read.table("covariatedata.txt",sep="\t", header=TRUE)
+pheno<- read.table("pheno.txt",sep="\t",quote="", header=TRUE)
+##save(data,pheno,file="data_pheno.Rdata")
+##load ("data_pheno.Rdata") ## fast loading data
+eqtl.index<- which(pheno$IncludeForGWAS==TRUE)     ## select samples
+data<-data[,eqtl.index]; pheno<-pheno[eqtl.index,]
+
+
+## Get individual covariates
+smoking<- pheno$Smoking_status
+smoking[which(is.na(smoking)==TRUE)]<-"3"  ## put missing value as "3"
+gender <- pheno$Gender
+age <- as.numeric(pheno$Age)
+diseases <- pheno$Diseases
+tid<- pheno$TID
+cov<- as.data.frame(cbind(smoking,gender,age,diseases,tid)) 
+temp<- as.character(pheno$Sample_ID)
+centers<- apply(t(temp),2, function (x) strsplit(x,"_")[[1]][1])
+rm(temp)
+
+## Perform regression analyses:
+
+lm.regression<- function(data,snp,cov,pc) {
+  
+  gender <- cov$gender
+  smoking <- cov$smoking
+  age   <- as.numeric(cov$age)
+  diseases <- cov$diseases
+  snp<- unlist(snp)
+  lm.reg <- lm(data~snp+gender+smoking+age+diseases+pc)
+  out<- summary(lm.reg)
+  beta<- out$coefficients[2,"Estimate"]
+  sd <-  out$coefficients[2,"Std. Error"]
+  pvalue <- out$coefficients[2, "Pr(>|t|)"]
+  output<- c(beta,sd,pvalue)
+  return(output)
+}
\ No newline at end of file