#This file contains R code to produce funnel plots for some genotyping datasets
#A more generic + better documented set of code is in progress!
#See snpfunnel.pdf for interpretation and exposition of funnel plots

#First example

hclean <- read.table("hakonarsenclean.txt", header=T)

pdf("snpfunnelex1.pdf",w=8,h=6)
plot(y=log(hclean$or), x=1/( (log(hclean$highci)-log(hclean$lowci)) / 2 / qnorm(0.975) )^2, xlab="1/std error^2", ylab="log OR est", ylim=2.5*c(-1,1)*max(log(hclean$or)) , main="")
lines(x=1:150, y=qnorm(0.975)/sqrt(1:150), lty=1)
lines(x=1:150, y=-qnorm(0.975)/sqrt(1:150), lty=1)
legend(x=20, y=0.7, c("p=0.05"), lty=1, bty="n")
abline(h=0, lty=3)
dev.off()

#Second example

elyall <- read.table("elybasic.txt")

glmfunc=function(x){
	summary(glm(matrix(as.numeric(x[c(3,5,7,9,11,13)]),3,2)~I(1:3), family=binomial))$coefficients[2,1:2]
	}

elylist1 <- t(apply(elyall, 1,FUN=function(x){summary(glm(matrix(as.numeric(x[c(3,5,7,9,11,13)]),3,2)~I(1:3), family=binomial))$coefficients[2,1:2]
}))

pdf("snpfunnelex2.pdf",w=8,h=6)
plot(x=1/elylist1[,2]^2, y=elylist1[,1], ylim=1.0*c(-1,1)*max(abs(elylist1[,1])), xlab="1/std error^2", ylab="log OR est", main="")
lines(x=1:130, y=qnorm(0.975)/sqrt(1:130), lty=1)
lines(x=1:130, y=-qnorm(0.975)/sqrt(1:130), lty=1)
lines(x=1:130, y=qnorm(0.99)/sqrt(1:130), lty=2)
lines(x=1:130, y=-qnorm(0.99)/sqrt(1:130), lty=2)
lines(x=1:130, y=qnorm(1-0.05/2/71)/sqrt(1:130), lty=3)
lines(x=1:130, y=-qnorm(1-0.05/2/71)/sqrt(1:130), lty=3)
abline(h=0, lty=4)
legend(x=60, y=1.3, lty=1:3, legend=c("p=0.05","p=0.01","p=0.05 (Bonferroni)"), bty="n")
text(x = 1/elylist1[c(39,59,63),2]^2, y = elylist1[c(39,59,63),1], labels = elyall[c(39,59,63),1], cex=0.7, pos=c(3,1,1))
dev.off()

#Third example

rheumall <- read.table("rheumsnp.txt")

#trend
rheumlist1 <- t(apply(rheumall, 1,FUN=function(x){
	summary(glm(matrix(as.numeric(x),3,2)~I(1:3), family=binomial))$coefficients[2,1:2]
	}
))

#dominant
rheumlist2 <- t(apply(rheumall, 1,FUN=function(x){
	summary(glm(matrix(as.numeric(x),3,2)~c(0,1,1), family=binomial))$coefficients[2,1:2]
	}
))

#recessive
rheumlist3 <- t(apply(rheumall, 1,FUN=function(x){
	summary(glm(matrix(as.numeric(x),3,2)~c(0,0,1), family=binomial))$coefficients[2,1:2]
	}
))

pdf("snpfunnelex3.pdf",w=8,h=6)
plot(x=1/rheumlist1[,2]^2, y=rheumlist1[,1], xlim=c(1,200), ylim=c(-.8,.8), xlab="1/std error^2", ylab="log OR est", main="", pch=22)
points(x=1/rheumlist2[,2]^2, y=rheumlist2[,1], pch=24)
points(x=1/rheumlist3[,2]^2, y=rheumlist3[,1], pch=19)
lines(x=1:200, y=qnorm(0.975)/sqrt(1:200), lty=1)
lines(x=1:200, y=-qnorm(0.975)/sqrt(1:200), lty=1)
lines(x=1:200, y=qnorm(0.995)/sqrt(1:200), lty=2)
lines(x=1:200, y=-qnorm(0.995)/sqrt(1:200), lty=2)
lines(x=1:200, y=qnorm(1-0.05/18/3/2)/sqrt(1:200), lty=3)
lines(x=1:200, y=-qnorm(1-0.05/18/3/2)/sqrt(1:200), lty=3)
abline(h=0,lty=4)
text(x=120, y=-0.6-0.07,"slc2F1", cex=0.7)
lines(x=c(120, 1/rheumlist1[7,2]^2), y=c(-0.6,rheumlist1[7,1]) )
lines(x=c(120, 1/rheumlist2[7,2]^2), y=c(-0.6,rheumlist2[7,1]) )
lines(x=c(120, 1/rheumlist3[7,2]^2), y=c(-0.6,rheumlist3[7,1]) )
text(x=70, y= 0.6+0.07,"slc5R1", cex=0.7)
lines(x=c(70, 1/rheumlist1[18,2]^2), y=c(0.6,rheumlist1[18,1]) )
lines(x=c(70, 1/rheumlist2[18,2]^2), y=c(0.6,rheumlist2[18,1]) )
lines(x=c(70, 1/rheumlist3[18,2]^2), y=c(0.6,rheumlist3[18,1]) )
legend(x=140, y=0.8, lty=c(1:3,0,0,0), legend=c("p=0.05","p=0.01","p=0.05 (Bonferroni)","dominant","recessive","trend"), bty="n", pch=c(32,32,32,24,19,22))
dev.off()