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


## DESCRIPTION
## Brief exploratory plots of data in the Sebida gene expression database.
# install.packages("corrplot")

## Load packages.
lapply(c(
  "plyr", 
  "tidyverse", ## Easy coding and fancy plots.
	"cowplot", ## Combine fancy plots into one figure.
  "dtplyr", ## dplyr-data.table compatibility
  "data.table",  ## Load big data.
  "car", ## Recode character variables.
  "corrplot" ## Correlation matrix plot.
  ), library, character.only = TRUE)


setwd("~/Dropbox/gwas_paper/accessory_data/")

## Consistency between studies is clumpy, from pca of m/f ratios.
## Gibson similar to McIntyre, both OreR.
## Brain - head studies similar.
## https://www.ncbi.nlm.nih.gov/pubmed/11726925
## https://www.ncbi.nlm.nih.gov/pubmed/16934145
## Innocenti LHM, Ayroles DGGRP, Stolc 2004 strain unknown, Wyman 2010 2xUSA, broadly similar from pca.
## https://www.ncbi.nlm.nih.gov/pubmed/15499012
## Also maybe Ayroles, Innocenti, check Ingleby.



## Load data.
expr_raw = fread("sebida_melanogaster_3.2.clean.txt", verbose = TRUE) %>%
	arrange(Full_Name)

head(expr_raw)
length(expr_raw)
length(expr_raw$Accession)


## Select columns containing male-female ratio results.
cow = dplyr::select(expr_raw, matches("MF"))
head(cow)


## Format data for plotting.
mcow = cow %>% gather(study_name, ratio) %>%
		mutate(sn =
					gsub(" et al\\. ", " ",
				gsub("MF", "",
			gsub("\\_"," ",
						study_name)))) %>%
	mutate(log_ratio = log(ratio)) %>%
	mutate(log_rat = gsub("-Inf", 0, log_ratio)) %>%
	mutate(lr = as.numeric(log_rat))

head(mcow)
unique(mcow$sn)
summary(mcow$lr)
## Problems with using tidyr spread :(.
# x = mcow %>% dplyr::select(sn, lr) %>% spread(key = sn, value = lr)
# head(x)



## Plot distributions of male-female ratio for each study (log-transformed)
expr_distr =	ggplot(mcow %>% select(sn, lr) %>% filter(lr != "NA")
		, aes(lr, fill = cut(lr, 30))) +
		geom_histogram(show.legend = FALSE) +
		scale_x_continuous(
			name = expression(italic(log)~"(M:F) gene expression level")) +
		scale_y_continuous(name = "Count (of genes)") +
		scale_fill_discrete(h = c(360, 200), l = 30) +
		facet_wrap("sn", scales = "free") +
		theme_bw(base_size = 4) +
		coord_cartesian(expand = c(.01, 0)) +
		theme(
			panel.grid = element_blank(),
			strip.background = element_rect(fill = "grey90", colour = "white"),
			legend.position = "none")




## Principle components analysis to see how study results group together.
## extracting the first two eigen vectors.
cow.cor = cor(cow, use = 'pair')
cow.cor[1:4,1:4]
cow.eig = eigen(cow.cor)$vectors[,1:2]
e1 = cow.eig[,1]
e2 = cow.eig[,2]



## Tedious modification of study names.
study_names1 <- c(
	"Meta-analysis", "Innocenti 2010", "Wyman 2010",
	"Ayroles 2009", "Goldman 2007 (head)",
	"McIntyre 2006 (OreR)", "McIntyre 2006 (2b)", "McIntyre 2006 (OReR+2b)",
	"Gibson 2004 (OreR)", "Gibson 2004 (2b)",
	"Gibson 2004 (OreR+2b)", "Stolc 2004", "Parisi 2003 (testes:ovaries)", "Parisi 2004 (whole)", "Parisi 2004 (testes:ovaries)",
	"Parisi gonadectomized", "Ranz 2003", "Ranz 2003 (D.simulans)",
	"(brain, Catalan)", "Huylmans (tubule)")
colnames(cow.cor) <- study_names1
rownames(cow.cor) <- study_names1


## Plotting correlation matrix. How to print with ggplots??
png("plots_sebida_corr_matrix.png", width = 900, height = 900, units = "px")
corr_plot = corrplot.mixed(cow.cor, tl.pos = "lt", tl.col = "black",
                           lower = "number", upper = "circle", order = "hclust")
dev.off()




## Format the study names using dplyr+gsub.
edat = data.table(
	study_name = row.names(cow.cor), e1, e2) %>%
	arrange(study_name) %>%
		separate(study_name, into = c("lead_author", "tmp2"), sep = " ", remove = FALSE) ## tedious.

head(edat)
tail(edat)
str(edat)



## Plot PCA results labelled by study.
pca_plot = ggplot(edat, aes(e1, e2, label = study_name, colour = lead_author)) +
		geom_point() +
		geom_text(hjust = 0, nudge_x = 0.01, nudge_y = 0.005, size = 1.5) +
		scale_x_continuous(name = "PC 1") + scale_y_continuous(name = "PC 2") +
		coord_cartesian(xlim = c(-.3, .1)) +
	theme_bw(base_size = 8) +
		theme(legend.position = "none", panel.grid = element_blank())



## Save ggplots to file.
 save_plot("plots_sebida.png",
          plot_grid(expr_distr, pca_plot, labels = "AUTO", 
                    ncol = 2, label_size = 8), base_height = 3, base_width = 6)


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