GMRP.R
In GMRP: GWAS-based Mendelian Randomization and Path Analyses

### R code from vignette source 'GMRP.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex(use.unsrturl=FALSE)


###################################################
### code chunk number 2: GMRP.Rnw:47-50
###################################################
#library(knitr)
set.seed(102)
options(width = 90)


###################################################
### code chunk number 3: GMRP.Rnw:53-54
###################################################
library(GMRP)


###################################################
### code chunk number 4: fmerge (eval = FALSE)
###################################################
## data1 <- matrix(NA, 20, 4)
## data2 <- matrix(NA, 30, 7)
## SNPID1 <- paste("rs", seq(1:20), sep="")
## SNPID2 <- paste("rs", seq(1:30), sep="")
## data1[,1:4] <- c(round(runif(20), 4), round(runif(20), 4), round(runif(20), 4), round(runif(20), 4))
## data2[,1:4] <- c(round(runif(30), 4),round(runif(30), 4), round(runif(30), 4), round(runif(30), 4))
## data2[,5:7] <- c(round(seq(30)*runif(30), 4), round(seq(30)*runif(30), 4), seq(30))
## data1 <- cbind(SNPID1, as.data.frame(data1))
## data2 <- cbind(SNPID2, as.data.frame(data2))
## dim(data1)
## dim(data2)
## colnames(data1) <- c("SNP", "var1", "var2", "var3", "var4")
## colnames(data2) <- c("SNP", "var1", "var2", "var3", "var4", "V1", "V2", "V3")
## data1<-DataFrame(data1)
## data2<-DataFrame(data2)
## data12 <- fmerge(fl1=data1, fl2=data2, ID1="SNP", ID2="SNP", A=".dat1", B=".dat2", method="No")


###################################################
### code chunk number 5: GMRP.Rnw:99-103
###################################################
data(cad.data)
#cad <- DataFrame(cad.data)
cad<-cad.data
head(cad)


###################################################
### code chunk number 6: GMRP.Rnw:106-110
###################################################
data(lpd.data)
#lpd <- DataFrame(lpd.data)
lpd<-lpd.data
head(lpd)


###################################################
### code chunk number 7: Step1 (eval = FALSE)
###################################################
## pvalue.LDL <- lpd$P.value.LDL
## pvalue.HDL <-lpd$P.value.HDL
## pvalue.TG <- lpd$P.value.TG
## pvalue.TC <- lpd$P.value.TC
## pv <- cbind(pvalue.LDL, pvalue.HDL, pvalue.TG, pvalue.TC)
## pvj <- apply(pv, 1, min)


###################################################
### code chunk number 8: Step2 (eval = FALSE)
###################################################
## beta.LDL <- lpd$beta.LDL
## beta.HDL <- lpd$beta.HDL
## beta.TG <- lpd$beta.TG
## beta.TC <- lpd$beta.TC
## beta <- cbind(beta.LDL, beta.HDL, beta.TG, beta.TC)


###################################################
### code chunk number 9: Step3 (eval = FALSE)
###################################################
## a1.LDL <- lpd$A1.LDL
## a1.HDL <- lpd$A1.HDL
## a1.TG <- lpd$A1.TG
## a1.TC <- lpd$A1.TC
## alle1 <- cbind(a1.LDL, a1.HDL, a1.TG, a1.TC)


###################################################
### code chunk number 10: Step4 (eval = FALSE)
###################################################
## N.LDL <- lpd$N.LDL
## N.HDL <- lpd$N.HDL
## N.TG <- lpd$N.TG
## N.TC <- lpd$N.TC
## ss <- cbind(N.LDL, N.HDL, N.TG, N.TC)
## sm <- apply(ss,1,sum)
## pcj <- round(sm/max(sm), 6)


###################################################
### code chunk number 11: Step5 (eval = FALSE)
###################################################
## freq.LDL<-lpd$Freq.A1.1000G.EUR.LDL
## freq.HDL<-lpd$Freq.A1.1000G.EUR.HDL
## freq.TG<-lpd$Freq.A1.1000G.EUR.TG
## freq.TC<-lpd$Freq.A1.1000G.EUR.TC
## freq<-cbind(freq.LDL,freq.HDL,freq.TG,freq.TC)


###################################################
### code chunk number 12: Step6 (eval = FALSE)
###################################################
## sd.LDL <- rep(37.42, length(pvj))
## sd.HDL <- rep(14.87, length(pvj))
## sd.TG <-rep(92.73, length(pvj))
## sd.TC <- rep(42.74, length(pvj))
## sd <- cbind(sd.LDL, sd.HDL, sd.TG, sd.TC)


###################################################
### code chunk number 13: Step7 (eval = FALSE)
###################################################
## hg19 <- lpd$SNP_hg19.HDL
## rsid <- lpd$rsid.HDL


###################################################
### code chunk number 14: Step8 (eval = FALSE)
###################################################
## chr<-chrp(hg=hg19)


###################################################
### code chunk number 15: Step9 (eval = FALSE)
###################################################
## newdata<-cbind(freq,beta,sd,pvj,ss,pcj)
## newdata<-cbind(chr,rsid,alle1,as.data.frame(newdata))
## dim(newdata)


###################################################
### code chunk number 16: Step10 (eval = FALSE)
###################################################
## hg18.d <- cad$chr_pos_b36
## SNP.d <- cad$SNP #SNPID
## a1.d<- tolower(cad$reference_allele)
## freq.d <- cad$ref_allele_frequency
## pvalue.d <- cad$pvalue
## beta.d <- cad$log_odds
## N.case <- cad$N_case
## N.ctr <- cad$N_control
## N.d <- N.case+N.ctr
## freq.case <- N.case/N.d


###################################################
### code chunk number 17: Step11 (eval = FALSE)
###################################################
## newcad <- cbind(freq.d, beta.d, N.case, N.ctr, freq.case)
## newcad <- cbind(hg18.d, SNP.d, a1.d, as.data.frame(newcad))
## dim(newcad)


###################################################
### code chunk number 18: Step12 (eval = FALSE)
###################################################
## varname <-c("CAD", "LDL", "HDL", "TG", "TC")


###################################################
### code chunk number 19: Step13 (eval = FALSE)
###################################################
## mybeta <- mktable(cdata=newdata, ddata=newcad, rt="beta", varname=varname, LG=1, Pv=0.00000005, Pc=0.979, Pd=0.979)
## dim(mybeta)
## beta <- mybeta[,4:8]   #  standard beta table for path analysis
## snp <- mybeta[,1:3]   #  snp data for annotation analysis
## beta<-DataFrame(beta)
## head(beta)


###################################################
### code chunk number 20: GMRP.Rnw:268-275
###################################################
data(beta.data)
beta.data<-DataFrame(beta.data)
CAD <- beta.data$cad
LDL <- beta.data$ldl
HDL <- beta.data$hdl
TG <- beta.data$tg
TC <- beta.data$tc


###################################################
### code chunk number 21: TwoScatterPlot
###################################################
par(mfrow=c(2, 2), mar=c(5.1, 4.1, 4.1, 2.1), oma=c(0, 0, 0, 0))
plot(LDL,CAD, pch=19, col="blue", xlab="beta of SNPs on LDL", ylab="beta of SNP on CAD", cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~LDL), col="red", lwd=2)
plot(HDL, CAD, pch=19,col="blue", xlab="beta of SNPs on HDL", ylab="beta of SNP on CAD", cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~HDL), col="red", lwd=2)
plot(TG, CAD, pch=19, col="blue", xlab="beta of SNPs on TG", ylab="beta of SNP on CAD",cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~TG), col="red", lwd=2)
plot(TC,CAD, pch=19, col="blue", xlab="beta of SNPs on TC", ylab="beta of SNP on CAD", cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~TC), col="red", lwd=2)


###################################################
### code chunk number 22: figTwoScatterPlot
###################################################
par(mfrow=c(2, 2), mar=c(5.1, 4.1, 4.1, 2.1), oma=c(0, 0, 0, 0))
plot(LDL,CAD, pch=19, col="blue", xlab="beta of SNPs on LDL", ylab="beta of SNP on CAD", cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~LDL), col="red", lwd=2)
plot(HDL, CAD, pch=19,col="blue", xlab="beta of SNPs on HDL", ylab="beta of SNP on CAD", cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~HDL), col="red", lwd=2)
plot(TG, CAD, pch=19, col="blue", xlab="beta of SNPs on TG", ylab="beta of SNP on CAD",cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~TG), col="red", lwd=2)
plot(TC,CAD, pch=19, col="blue", xlab="beta of SNPs on TC", ylab="beta of SNP on CAD", cex.lab=1.5, cex.axis=1.5, cex.main=2)
abline(lm(CAD~TC), col="red", lwd=2)


###################################################
### code chunk number 23: path (eval = FALSE)
###################################################
## data(beta.data)
## mybeta <- DataFrame(beta.data)
## mod <- CAD~LDL+HDL+TG+TC
## pathvalue <- path(betav=mybeta, model=mod, outcome="CAD")


###################################################
### code chunk number 24: GMRP.Rnw:312-319
###################################################
mypath <- matrix(NA,3,4)
mypath[1,] <- c(1.000000, -0.066678, 0.420036, 0.764638)
mypath[2,] <- c(-0.066678, 1.000000, -0.559718, 0.496831)
mypath[3,] <- c(0.420036, -0.559718, 1.000000, 0.414346)
colnames(mypath) <- c("LDL", "HDL", "TG", "path")
mypath<-as.data.frame(mypath)
mypath


###################################################
### code chunk number 25: GMRP.Rnw:332-333
###################################################
library(diagram)


###################################################
### code chunk number 26: PathDiagram
###################################################
pathdiagram(pathdata=mypath, disease="CAD", R2=0.988243, range=c(1:3))


###################################################
### code chunk number 27: figPathDiagram
###################################################
pathdiagram(pathdata=mypath, disease="CAD", R2=0.988243, range=c(1:3))


###################################################
### code chunk number 28: GMRP.Rnw:365-373
###################################################
pathD<-matrix(NA,4,5)
pathD[1,] <- c(1,	-0.070161, 0.399038, 0.907127, 1.210474)
pathD[2,] <- c(-0.070161,	1, -0.552106, 0.212201, 0.147933)
pathD[3,] <- c(0.399038, -0.552106, 1, 0.44100, 0.64229)
pathD[4,] <- c(0.907127, 0.212201, 0.441007, 1, -1.035677)
colnames(pathD) <- c("LDL", "HDL", "TG", "TC", "path")
pathD<-as.data.frame(pathD)
pathD


###################################################
### code chunk number 29: PathDiagram2
###################################################
pathdiagram2(pathD=pathD,pathO=mypath,rangeD=c(1:4),rangeO=c(1:3),disease="CAD", R2D=0.536535,R2O=0.988243)


###################################################
### code chunk number 30: figPathDiagram2
###################################################
pathdiagram2(pathD=pathD,pathO=mypath,rangeD=c(1:4),rangeO=c(1:3),disease="CAD", R2D=0.536535,R2O=0.988243)


###################################################
### code chunk number 31: GMRP.Rnw:397-400
###################################################
data(SNP358.data)
SNP358 <- as.data.frame(SNP358.data)
head(SNP358)


###################################################
### code chunk number 32: GMRP.Rnw:404-405
###################################################
library(graphics)


###################################################
### code chunk number 33: ChromHistogram
###################################################
snpPositAnnot(SNPdata=SNP358,SNP_hg19="chr",main="A")


###################################################
### code chunk number 34: figChromHistogram
###################################################
snpPositAnnot(SNPdata=SNP358,SNP_hg19="chr",main="A")


###################################################
### code chunk number 35: GMRP.Rnw:435-438
###################################################
data(SNP368annot.data)
SNP368<-as.data.frame(SNP368annot.data)
SNP368[1:10, ]


###################################################
### code chunk number 36: GMRP.Rnw:442-443
###################################################
library(plotrix)


###################################################
### code chunk number 37: PIE3D1
###################################################
ucscannot(UCSCannot=SNP368,SNPn=368)


###################################################
### code chunk number 38: figPIE3D1
###################################################
ucscannot(UCSCannot=SNP368,SNPn=368)


###################################################
### code chunk number 39: PIE3D2
###################################################
ucscannot(UCSCannot=SNP368,SNPn=368,A=3,B=2,C=1.3,method=2)


###################################################
### code chunk number 40: figPIE3D2
###################################################
ucscannot(UCSCannot=SNP368,SNPn=368,A=3,B=2,C=1.3,method=2)


###################################################
### code chunk number 41: GMRP.Rnw:495-496
###################################################
sessionInfo()