parent
a5d28ace7b
commit
b0a2f65c74
53
eQTL_Scripts.md
Normal file
53
eQTL_Scripts.md
Normal file
@ -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)
|
||||
}
|
Loading…
x
Reference in New Issue
Block a user