library(PPMR)
|
|
|
|
| #load the simulated scaled genenotype matrix in eQTL data (e.g. the 465*50 cis-SNP genotype matrix of BACE1 gene from GEUVADIS data)
|
| zx<-read.table("zx.txt")
|
| zx<-as.matrix(zx)
|
|
|
| #load the simulated scaled genenotype matrix in GWAS data (e.g. the same cis-SNPs from GERA data)
|
| zy<-read.table("zy.txt")
|
| zy<-as.matrix(zy)
|
|
|
| #load the simulated exposure or gene expression vector
|
| x<-read.table("x.txt")
|
| x<-as.vector(x[[1]])
|
|
|
| #load the simulated phenotype vector
|
| y<-read.table("y.txt")
|
| y<-as.vector(y[[1]])
|
|
|
| #run the PPMR model using PMR_individual function
|
| fmH1 = PMR_individual(x, y, zx, zy,gammain=0,alphain = 0,max_iterin =1000, epsin=1e-5)
|
|
|
| #get the estimate of the causal effect
|
| alpha<-fmH1$alpha
|
|
|
| #get the estimate of the pleiotropy effect
|
| gamma<-fmH1$gamma
|
|
|
| fmH0gamma = PMR_individual(x, y, zx, zy,gammain=1,alphain = 0,max_iterin =1000, epsin=1e-5)
|
| fmH0alpha = PMR_individual(x, y, zx, zy,gammain=0,alphain = 1,max_iterin =1000, epsin=1e-5)
|
| loglikH1=max(fmH1$loglik,na.rm=T)
|
| loglikH0gamma=max(fmH0gamma$loglik,na.rm=T)
|
| loglikH0alpha=max(fmH0alpha$loglik,na.rm=T)
|
| stat_alpha = 2 * (loglikH1 - loglikH0alpha)
|
|
|
| #get the pvalue for the causal effect
|
| pvalue_alpha = pchisq(stat_alpha,1,lower.tail=F)
|
| stat_gamma = 2 * (loglikH1 - loglikH0gamma)
|
|
|
| #get the pvalue for the pleiotropy effect
|
| pvalue_gamma = pchisq(stat_gamma,1,lower.tail=F)
|
