# Script to create Manhattan and QQ plots (including lambda value).

# Get the qqman library (install if needed before: install.packages("qqman")

rm(list=ls())

library(qqman)

# Indicate working directory 

setwd("XXX") # set your working directory; e.g. ~/session3/casecontrol or e.g. ~/session3/continuous

Data <- read.table('XXX', header=T) # select you file; e.g. plot.adclean.cc.logistic.txt or plot.adclean.cont.linear.txt
head(Data)

# make sure the data is prepared so there is no missing data and the key variables (CHR, BP, P) are numeric.

#Data<- Data[complete.cases(Data),] # make sure there's no missing data
#Data$CHR[which(Data$CHR=="X")]<-23 # if you have data on chrX, it'll have to be changed to a numeric format.
Data$CHR<-as.numeric(Data$CHR)
Data$BP<-as.numeric(Data$BP)
Data$P<-as.numeric(Data$P)


# to generate qq plot (give a name to the file you want to create; e.g. qq_cc.jpg or qq_cont.jpg)

jpeg("qq_XXX.jpg", width = 400, height = 400, quality=75, bg = "white", type = "cairo")
alpha<-median(qchisq(1-Data$P,1))/qchisq(0.5,1)
qq(Data$P)
text(0.7,5, paste("lambda","=",  signif(alpha, digits = 3)) )
dev.off()


# to generate manhattan plot (give a name to the file you want to create; e.g. manhattan_cc.jpg or manhattan_cont.jpg)

jpeg("manhattan_XXX.jpg", width = 800, height = 400, quality=75, bg = "white", type = "cairo")
manhattan(Data, chr = "CHR", bp = "BP", p = "P", snp = "SNP")
dev.off()