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