library("qqman")

#https://cran.r-project.org/web/packages/qqman/vignettes/qqman.html

#SELECT FOLDER TO PRACTICUM FILES!

#File containing info for creating MH plot
gwas_example_bim <- read.delim("gwas-example-bim.txt", header = T)

#we need the following columns
gwas_example_mh <- gwas_example_bim[c(1,2,4)]

results1 <- read.delim("results1.txt", header = T)
head(results1)
tail(results1)
colnames(results1)[3] <- "BP"

results2 <- read.delim("results2.txt", header = T, sep = " ")
head(results2)
tail(results2)
colnames(results2)[3] <- "BP"

#read in meta-analysis files
meta_SE <- read.delim("meta_res_SE1.txt", header = T, sep = "\t")
colnames(meta_SE)[1] <- "SNP"
head(meta_SE)


meta_Z <- read.delim("meta_res_Z1.txt", header= T, sep = "\t")
colnames(meta_Z)[1] <- "SNP"
head(meta_Z)

#merge file to include CHR and BP
meta_Z_MHplot <- merge(gwas_example_mh, meta_Z, by = "SNP")
colnames(meta_Z_MHplot)[8] <- "P"
head(meta_Z_MHplot)

#make MH plots

png(file = "MHplots.png", width = 8000, height = 5000, units='px', res = 300, pointsize=17)
layout(matrix(c(1,2,3), 3, 1, byrow = TRUE))
manhattan(results1, main="results1")
manhattan(results2, main="results2")
manhattan(meta_Z_MHplot, main="meta_Z_MHplot")
dev.off()