# 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.~/gwas2/continuous/ or e.g. ~/gwas2/casecontrol/, or use "Session" in the drop-down menu above to set your directory

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

colnames(Data)<-c("CHROM", "POS", "ID", "P") # give it the right column names
head(Data)


# make sure the data is prepared so there is no missing data and the key variables (CHR, BP, P) are numeric.
# How do you check the structure of the data?
XXX(Data) # What is the command you can use

#Data<- Data[complete.cases(Data),] # make sure there's no missing data
#Data$CHROM[which(Data$CHROM=="X")]<-23 # if you have data on chrX, it'll have to be changed to a numeric format.
Data$CHROM<-as.numeric(Data$CHROM)
Data$POS<-as.numeric(Data$POS)
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) ### This is getting the observed lambda in your data
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 = "CHROM", bp = "POS", p = "P", snp = "ID")
dev.off()