#
##
###
## Start log
logname = "log_r_basic_gwas.txt"
mylog = file(logname, open = "wt")
sink(mylog, append = TRUE, type = "message")
Sys.time() ; sessionInfo()


## DESCRIPTION
## Reads GWAS results for each sex.
## Generates manhattan, p-value distribution and QQ-plots, for each sex.


## Load packages
library(plyr) ## For cumulative calculation of FDR values.
library(tidyverse) ## Data formatting, and fancy plots.
library(cowplot) ## Put fancy plots in one figure.
library(dtplyr) ## dplyr-data.table compatibility
library(data.table) ## Loading big data.
library(car) ## Recoding character variables.



## Load Female and Male GWAS results, combine, and format the phenotype identifiers.
df_raw = 	do.call("rbind", lapply(
	list.files(pattern = "*.gwas_res.txt", recursive = TRUE),
	function(x) cbind(Phenotype = gsub(".gwas_res.txt", "", x),
		fread(x, verbose = TRUE, sep = " ")))) %>%
	mutate(Phenotype = gsub("lhm_hc_2017.", "", Phenotype))



## General plot settings.
theme_w = function() {
	theme_bw(base_size = 10) +
		theme(
			axis.ticks = element_line(size = .3),
			panel.grid = element_blank(),
			strip.background = element_rect(color = "grey90", fill = "grey93", size = .3))}



## Plot p-value distributions.
pdist = ggplot(df_raw, aes(P, colour = Phenotype)) +
	geom_histogram(fill = "white", bins = 60) +
	scale_colour_manual(values = c("red4", "blue4")) +
	scale_x_continuous(name = "p-value category", expand = c(.02,0),
		limits = c(0, 1),
		breaks = seq(0,1, .1),
		labels = c(0,"",.2,"",.4,"",.6,"", .8, "", 1)) +
	scale_y_continuous(name = "Count of tests", expand = c(.01,.01))+
	facet_wrap("Phenotype") +
	theme_w() +
	theme(
		axis.line.x = element_line(size = .3, colour = "black"),
		axis.line.y = element_line(size = .3, colour = "black"),
		legend.position = "none",
		plot.margin = unit(c(2,.5,2,.5),"lines"))



## Manhattan plot
mh = ggplot(df_raw %>% select(
	Phenotype, CHR, BP, P) , aes(BP/1e+06, -log10(P), colour = Phenotype)) +
	geom_point(shape = 21, stroke = .2, size = .3) +
	scale_colour_manual(values = c("red4", "blue4")) +
	scale_x_continuous(name = "Genomic position (Mb)", expand = c(0.01,0),
		breaks = seq(0,40,5), labels = c(0, "", 10, "", 20, "", 30, "", 40)) +
	scale_y_continuous(name = "-log10.p (basic GWAS test)", expand = c(0,0.02), limits = c(0, 7.1),
		breaks = seq(0,7,1), labels = c(0,"",2,"",4,"",6,"")) +
	facet_grid(Phenotype~CHR, space = "free_x", scale = "free_x") +
	theme_w() +
	theme(
		axis.line.x = element_line(size = .2, colour = "black"),
		axis.line.y = element_line(size = .2, colour = "black"),
		legend.position = "none",
		strip.placement = "outside")




## FUNCTION : Calculate expected p-values (+95% CI) for QQ plot.
fun_cal_qq = function(df, mytest){
	x = df %>%
		sample_frac(1) %>%
		filter(Phenotype == mytest) %>%
		select(P) %>%
		arrange(P) %>%
		mutate(expP = (1:length(P) / length(P))) %>%
		mutate(log10p = -log10(P)) %>%
		mutate(explog10p = -log10(expP))
	c95 = as.numeric(rep(0,length(x$P)))
	c05 = c95
	for (i in 1:length(x$P)) {
		c95[i] = as.numeric(-log10(qbeta(0.05,i,length(x$P) - i + 1)))
		c05[i] = as.numeric(-log10(qbeta(0.95,i,length(x$P) - i + 1)))}
	final = as.data.table(cbind(x, c05 = as.numeric(c05), c95 = as.numeric(c95), Phenotype = mytest))
	rm(x) ; rm(c05) ; rm(c95) ; return(final) }



Sys.time() ; print("Making data for QQ plots...")

df_qq = do.call("rbind", list(
	fun_cal_qq(df_raw, "Female"),
	fun_cal_qq(df_raw, "Male")))

Sys.time() ; print("Completed QQ plot data.")
head(df_qq)


## QQ plot.
qqplot =	ggplot(df_qq, aes(explog10p, log10p)) +
	geom_ribbon(aes(ymin = as.numeric(c05), ymax = as.numeric(c95)),
		fill = "grey70", alpha = .8, linetype = 0) +
	geom_abline(intercept = 0, slope = 1, size = .1, colour = "grey40") +
	geom_point(aes(colour = Phenotype), size = .6, stroke = .2) +
	scale_colour_manual(values = c("red4", "blue4")) +
	coord_fixed(ratio = 1, expand = c(0,0), xlim = c(0,7.1), ylim = c(0,7.1)) +
	scale_x_continuous(name = "Expected -log10.p", breaks = seq(0,7,1), labels = c(0,"",2,"",4,"",6,"")) +
	scale_y_continuous(name = "Observed -log10.p", breaks = seq(0,7,1), labels = c(0,"",2,"",4,"",6,"")) +
	facet_wrap("Phenotype") +
	theme_w() +
	theme(legend.position = "none")




Sys.time() ; print("Starting plots....")

save_plot("plot_basic_gwas_lhm.png",
	ggdraw() +
		draw_plot(mh, 0, .5, 1, .5) +
		draw_plot(pdist, 0, 0, .5, .5) +
		draw_plot(qqplot, .5, 0, .5, .5) +
		draw_plot_label(c("A", "B", "C"), c(0, 0, .55), c(.98, .48, .48), size = 10)
	, base_height = 6, base_width = 8)

Sys.time() ; print("Plots saved to file.")



## Links
## Morrow lab, Sussex Uni http://www.sussex.ac.uk/lifesci/morrowlab/
## Project https://zenodo.org/communities/sussex_drosophila_gwas/
## Plink https://www.cog-genomics.org/plink2


ls(all.names = TRUE)
rm(list = ls())
Sys.time()
print("William P. Gilks, wpgilks@gmail.com, University of Sussex 2017. End of script")
sink() ; unlink(mylog)
###
##
#