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


## Description:
## Reads GWAS results for each sex, and bivariate test.
## Calculates FDR, based on all and independent tests.
## Outputs plots of FDR results, manhattan and QQ-plots.
## Outputs tables of significant SNPs and summary counts.


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.
library(fdrtool) ## FDR calculations.


## General settings for plots.
theme_fdr = function() {
	theme_bw(base_size = 6) +
		theme(
			axis.ticks = element_line(size = .2),
			legend.key.size = unit(.8, 'lines'),
			legend.key.width = unit(1, 'lines'),
			legend.background = element_rect(fill = "white", size = 0.25, linetype = 1, colour = "grey"),
			legend.key = element_rect(colour = "white"),
			panel.border = element_rect(color = "black", size = .2),
			panel.grid.minor = element_blank(),
			panel.grid.major = element_line(size = .1),
			plot.margin = margin(6,6,6,6),
			strip.background = element_rect(color = "grey90", fill = "grey98", size = .2))}



## Load single-sex GWAS results
rawF = fread(verbose = TRUE, "lhm_hc_2017.Female.gwas_res.txt")
rawM = fread(verbose = TRUE, "lhm_hc_2017.Male.gwas_res.txt")


## Load and format bivariate GWAS results
rawB = fread("resultsJoint.wide.txt", verbose = TRUE) %>%
	filter(label == "pval") %>%
	select(chr_pos = rsid, JointModel) %>%
	separate(chr_pos, c("pchr", "BP"), sep = "_") %>%
	mutate(CHR = car::recode(pchr, "1 = 'chr2L'; 2 = 'chr2R'; 3 = 'chr3L'; 4 = 'chr3R'; 5 = 'chrX'; 6 = 'chr4'; 7 = 'chrM'")) %>%
	mutate(BP = as.numeric(BP))

head(rawB)


## Check number of rows are equal.
ifelse(nrow(rawF) == nrow(rawM), "Female and male row number is equal :)", "Unequal female and male row number :(")
ifelse(nrow(rawF) == nrow(rawB), "Female and Bivariate row number is equal :)", "Unequal Female and Bivariate row number :(")


## Load list of LD-independent/tagging SNPs
indepSNPs = fread("data_tagging_snps.txt", header = FALSE) %>% dplyr::rename(SNP = V1) %>% mutate(tagging_snp = "yes")
head(indepSNPs)


## Combine p-values from the different tests into one table.
pdat = full_join(inner_join(
	rawF %>% select(SNP, CHR, BP, Female = P),
	rawM %>% select(SNP, CHR, BP, Male = P)),
	rawB) %>% select(-pchr) %>% rename(Bivariate = JointModel)

head(pdat)


## Make two tables, one for independent snps, and one for all snps.
pall = full_join(pdat, indepSNPs) %>% select(-c(SNP, CHR, BP, tagging_snp)) %>% gather(test, pval)
pindep = inner_join(pdat, indepSNPs) %>% select(-c(SNP, CHR, BP, tagging_snp)) %>% gather(test, pval)

head(pindep)




## Function to calculate GWAS FDR values, using fdrtool package.
calFdr = function(df, phenotype){
	dat = df %>% filter(test == phenotype) %>% select(pval) %>% arrange(pval)
	res = fdrtool(dat$pval, statistic = "pvalue", plot = FALSE)
	tab1 = data.table(pval = res$pval, qval = res$qval, lfdr = res$lfdr) %>%
		mutate(Phenotype = phenotype)
	rm(res)
	tab2 = data.table(dat, tab1) %>% select(-pval)
	rm(dat) ; rm(tab1) ; return(tab2)}


## Calculate FDR for LD-independent SNPs.
fdrIndep = do.call("rbind", list(calFdr(pindep, "Female"),
	calFdr(pindep, "Male"),
	calFdr(pindep, "Bivariate"))) %>%
	mutate(log10pval = -log10(pval))

## FDR using all SNPs.
fdrAll = do.call("rbind", list(calFdr(pall, phenotype = "Female"),
	calFdr(pall, "Male"),
	calFdr(pall, "Bivariate"))) %>%
	mutate(log10pval = -log10(pval))





## Function for p-value distribution plots.
plotPdistr = function(df, myxlab) {
	ggplot(df, aes(pval, colour = test, fill = test)) +
		geom_histogram(size = .1, alpha = .4) +
		scale_colour_manual(values = c("grey20", "red4", "blue4")) +
		scale_fill_manual(values = c("grey20", "red4", "blue4")) +
		scale_y_continuous(expand = c(.03,0), name = "Count") +
		scale_x_continuous(expand = c(0.01,0),
			labels = c("0", "0.25", "0.5", "0.75", "1"), name = myxlab) +
		facet_wrap("test") +
		theme_fdr() +
		theme(
			axis.line.x = element_line(size = .3, colour = "black"),
			axis.line.y = element_line(size = .3, colour = "black"),
			legend.position = "none")}


## Make p-value distribution plots.
pDistrAll = plotPdistr(pall, "p-value category (All SNPs)")
pDistrInd = plotPdistr(pindep, "p-value category (Independent SNPs)")





## Function for plotting log10p against q-value.
plot_pVq = function(df, myxlab){
	ggplot(df %>% select(log10pval, qval, Phenotype),
		aes(log10pval, qval, colour = Phenotype)) +
		geom_point(size = .4, alpha = .7) +
		scale_colour_manual(values = c("grey20", "red4", "blue4")) +
		scale_x_continuous(name = myxlab, breaks = seq(2,7,1),
			labels = c(2, 3, 4, 5, 6 ,"")) +
		scale_y_continuous(name = expression(italic(q)[FDR]), expand = c(0,0)) +
		coord_cartesian(ylim = c(0,.35), xlim = c(2, 7)) +
		facet_wrap("Phenotype") +
		theme_fdr() +
		theme(
			axis.line.x = element_line(size = .3, colour = "black"),
			axis.line.y = element_line(size = .3, colour = "black"),
			legend.position = "none")}


## Make log10p verses q-value plots.
pqAll = plot_pVq(fdrAll, "-log10.p (All SNPs)")
pqInd = plot_pVq(fdrIndep, "-log10.p (LD-independent SNPs)")




## Function for generating cumulative counts and plot of SNPs at each FDR q-value.
plotCumq = function(df, myxmax, myymax, myxlab) {
	dat = ddply(df,.(Phenotype), transform, len = length(qval))
	plot = ggplot(dat, aes(x = qval, colour = Phenotype)) +
		geom_step(aes(len = len, y = ..y.. * len), stat = "ecdf") +
		coord_cartesian(xlim = c(0, myxmax), ylim = c(0, myymax)) +
		scale_colour_manual(values = c("grey20", "red4", "blue4")) +
		scale_x_continuous(name = myxlab, breaks = seq(0,.3,.05),
			labels = c("0", "", "0.1", "", "0.2", "", "0.3")) +
		scale_y_continuous(name = "Cumulative SNP count") +
		facet_wrap("Phenotype") +
		theme_fdr() +
		theme(
			axis.line.x = element_line(size = .3, colour = "black"),
			axis.line.y = element_line(size = .3, colour = "black"),
			legend.position = "none") ; return(plot)}


## Make cumulatative q-value plots.
cumqInd = plotCumq(fdrIndep, 0.2, 25, "q-value (LD-independent)")
cumqAll = plotCumq(fdrAll, 0.35, 15000, "q-value (All SNPs)")







## Function to calculate expected p-values (+95% CI).
calQQ = function(df, mytest){
	x = df %>%
		filter(test == mytest) %>%
		arrange(pval) %>%
		mutate(expP = (1:length(pval) / length(pval))) %>%
		mutate(log10p = -log10(pval)) %>%
		mutate(explog10p = -log10(expP))
	c95 = as.numeric(rep(0,length(x$pval)))
	c05 = c95
	for (i in 1:length(x$pval)) {
		c95[i] = as.numeric(-log10(qbeta(0.05,i,length(x$pval) - i + 1)))
		c05[i] = as.numeric(-log10(qbeta(0.95,i,length(x$pval) - i + 1)))}
	final = as.data.table(cbind(x, c05 = as.numeric(c05), c95 = as.numeric(c95)))
	rm(x) ; rm(c05) ; rm(c95) ; return(final) }



## Make QQ-plot data: LD-independent SNPs.
qqInd = do.call("rbind", list(calQQ(pindep, "Female"),
	calQQ(pindep, "Male"), calQQ(pindep, "Bivariate")))


## Make QQ-plot data: All SNPs.
qqAll = do.call("rbind", list(calQQ(pall, "Female"),
	calQQ(pall, "Male"), calQQ(pall, "Bivariate")))



## Function for QQ plots.
plotQQ = function(df, myxlab){
	ggplot(df %>% select(test, log10p, explog10p, c05, c95),
		aes(explog10p, log10p, colour = test)) +
		geom_ribbon(aes(ymin = c05, ymax = c95), fill = "grey70", alpha = .6, linetype = 0) +
		geom_abline(intercept = 0, slope = 1, size = .1, colour = "grey20") +
		geom_point(size = .6, fill = "white", stroke = .2) +
		scale_colour_manual(values = c("grey20", "red4", "blue4")) +
		# coord_fixed(ratio = 1, expand = c(0,0), xlim = c(0,7.1), ylim = c(0,7.1)) +
		scale_x_continuous(name = myxlab, 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("test") +
		theme_fdr() +
		theme(
			axis.line.x = element_line(size = .2, colour = "black"),
			axis.line.y = element_line(size = .2, colour = "black"),
			legend.position = "none",
			panel.grid = element_blank())}


## Make QQ plots.
qqIndep_plot = plotQQ(qqInd, "Expected -log10.p (LD-independent)")
qqAll_plot = plotQQ(qqAll, "Expected -log10.p (All SNPs)")





## Save FDR and QQ plots to file.
save_plot("plot_fdr_gwas_lhm.png",
	plot_grid(
		pDistrAll, pDistrInd,
		pqAll, pqInd,
		cumqAll, cumqInd,
		qqAll_plot, qqIndep_plot,
		align = "hv", labels = "AUTO",
		label_size = 6, nrow = 4), base_width = 7, base_height = 5)





## Function to get highest p-value for a given FDR threshold.
getPt = function(df, phenotype, qthreshold){
	x = df %>% filter(Phenotype == phenotype) %>%
		filter(qval < qthreshold)
	y = max(x$pval) ; return(y)}


## Get p-values corresponding to FDR < 0.1 for each phenotype.
thresholds = data.frame(test = c("Bivariate", "Female", "Male"),
	pval = as.numeric(c(
		getPt(fdrIndep, "Bivariate", 0.1),
		getPt(fdrIndep, "Female", 0.1),
		getPt(fdrIndep, "Male", 0.1)))) %>%
	mutate(log10p = -log10(pval))

head(thresholds)



## Make Manhattan plots for each phenotype with lines at FDR<0.1.
mh = ggplot(pdat %>% gather(test, pval, -c(CHR, BP, SNP)),
	aes(BP/1e+06, -log10(pval), colour = test)) +
	geom_point(size = .1, alpha = .8, shape = 20) +
	geom_hline(data = thresholds, aes(yintercept = log10p),
		linetype = 1, size = .5, alpha = .3, colour = "green4") +
	scale_colour_manual(values = c("grey20", "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 = "Significance of association (-log10.p)", expand = c(0,0.02),
		limits = c(0, 7.1), breaks = seq(0,7,1), labels = c(0,"",2,"",4,"",6,"")) +
	facet_grid(test~CHR, space = "free_x", scales = "free_x", switch = "both") +
	theme_fdr() +
	theme(
		axis.line.x = element_line(size = .2, colour = "black"),
		axis.line.y = element_line(size = .2, colour = "black"),
		legend.position = "none",
		panel.border = element_blank(),
		panel.spacing.x = unit(c(0,.3,0,.3), "lines"),
		plot.margin = unit(c(.3,.5,0,.15),"lines"),
		strip.background = element_blank(),
		strip.placement = "outside")


## Save Manhattan plots to file.
save_plot("plot_manhattans_gwas_lhm.png", mh, base_width = 5, base_height = 3)





## Table for plotting and printing to file, and graphing. Classify SNPs by significance.
gwas_res = full_join(inner_join(
	rawF %>% select(SNP, CHR, BP, betaF = BETA, pF = P, gene, snp_funct, rs_id, MAF),
	rawM %>% select(SNP, CHR, BP, betaM = BETA, pM = P)),
	rawB) %>%
	select(CHR, BP, SNP, rs_id, betaF, pF, betaM, pM, pB = JointModel, MAF, gene, snp_funct) %>%
	mutate(sigF = ifelse(pF < thresholds[2,2], "yes", "no")) %>%
	mutate(sigM = ifelse(pM < thresholds[3,2], "yes", "no")) %>%
	mutate(sigB = ifelse(pB < thresholds[1,2], "yes", "no")) %>%
	mutate(Class = paste(sigF, sigM, sigB, sep = "_")) %>%
	mutate(Class = car::recode(Class,
		"'no_no_no' = 'None'; 'yes_no_no' = 'Female'; 'no_no_yes' = 'Bivariate'; 'no_yes_no' = 'Male'; 'yes_yes_no' = 'Female+Male'; 'yes_no_yes' = 'Female+Bivariate'; 'no_yes_yes' = 'Male+Bivariate'; 'yes_yes_yes' = 'All'"))

head(gwas_res)
names(gwas_res)



## Scatter plot of Female and Male betas, coloured by bivariate significance.
beta_plot =	ggplot() +
	geom_point(data = gwas_res %>% filter(Class == "None"),
		aes(betaF, betaM, colour = Class, shape = Class), size = .2) +
	geom_point(data = gwas_res %>% filter(Class != "None"),
		aes(betaF, betaM, colour = Class, shape = Class), stroke = .3, size = .8) +
	scale_colour_manual(
		name = "Significance class:",
		guide = guide_legend(nrow = 1, keywidth = .5, keyheight = .5, label.hjust = .8),
		values = c("green4", "red4", "red1", "blue4", "blue1", "grey85"),
		labels = c("Bivariate     ", "Female     ", "Female+Bivariate     ", "Male     ",
			"Male+Bivariate     ", "None     ")) +
	scale_shape_manual(
		name = "Significance class:",
		guide = guide_legend(nrow = 1, override.aes = list(size = 1, stroke = .4, label.hjust = .8), keywidth = .5, keyheight = .5),
		values = c(2, 1, 1, 0, 0, 16),
	labels = c("Bivariate     ", "Female     ", "Female+Bivariate     ", "Male     ",
			"Male+Bivariate     ", "None     ")) +
	geom_hline(yintercept = 0, size = .05, colour = "black") +
	geom_vline(xintercept = 0, size = .05, colour = "black") +
	coord_fixed(ratio = 1) +
	scale_x_continuous(name = expression("Effect-size in Females ("*beta[F]*")") , expand = c(.01,.01), limits = c(-.4, .4)) +
	scale_y_continuous(name = expression("  Effect-size in Males ("*beta[M]*")") , expand = c(.01,.01), limits = c(-.4, .4)) +
	facet_grid(~CHR, space = "fixed", scales = "fixed") +
	theme_fdr() +
	theme(
		axis.line.x = element_line(size = .2, colour = "black"),
		axis.line.y = element_line(size = .2, colour = "black"),
		legend.background = element_rect(size = .05),
		legend.margin = margin(rep(1.5, 4)),
		legend.box.margin = margin(-7,-7,-7,-7),
		legend.position = "bottom",
		legend.title = element_text(size = 5),
		panel.grid = element_blank(),
		panel.spacing.x = unit(.5, "lines"))

## Save effect-size scatter plot to file.
save_plot("plot_betascatter.png", beta_plot, base_width = 6, base_height = 2)




## Write result tables to file.
fwrite(gwas_res, "data_Alltests_AllSNPs.txt", sep = " ")


## Get tables of significant SNPs for each group only.
sigBivariate = gwas_res %>% filter (sigB == "yes")
sigFemale = gwas_res %>% filter (sigF == "yes")
sigMale = gwas_res %>% filter (sigM == "yes")


## Write significant SNPs to file.
fwrite(sigBivariate, "data_Bivariate_SigSNPs.txt", sep = " ")
fwrite(sigFemale, "data_Female_SigSNPs.txt", sep = " ")
fwrite(sigMale, "data_Male_SigSNPs.txt", sep = " ")


## Create summary table. Top SNP p-value, counts of significant SNPs, and genes.
countsSummary = full_join(full_join(
	data.frame(do.call("rbind", list(
		c("Female", min(gwas_res$pF)),
		c("Male", min(gwas_res$pM)),
		c("Bivariate", min(gwas_res$pB))))) %>% rename(Test = X1, Min_P = X2),
	data.frame(do.call("rbind", list(
		c("Female", nrow(sigFemale)),
		c("Male", nrow(sigMale)),
		c("Bivariate", nrow(sigBivariate))))) %>% rename(Test = X1, `Significant_SNPs` = X2)),
	data.frame(do.call("rbind", list(
		c("Female", nrow(sigFemale %>% group_by(gene) %>% slice(which.min(pF)))),
		c("Male", nrow(sigMale %>% group_by(gene) %>% slice(which.min(pM)))),
		c("Bivariate", nrow(sigBivariate %>% group_by(gene) %>% slice(which.min(pB))))))) %>%
		rename(Test = X1, `Significant_Genes` = X2))
countsSummary


## Write counts summary to file.
fwrite(countsSummary, "data_results_summary.txt", sep = " ")




## Morrow lab, Sussex Uni http://www.sussex.ac.uk/lifesci/morrowlab/
## Zenodo community page https://zenodo.org/communities/sussex_drosophila_gwas/
## FDR with 'fdrtool' https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-9-303
## Much of the above code could be shortened with use of a custom package, loops, and lapply.
## Completed to the sound of Mozart's Don Giovanni.

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