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


## Description
## Reads results file from Plink IBD-PCA analysis.
## See https://www.cog-genomics.org/plink/1.9/ibd
## Assigns population sub-groups to individuals for later Fst analysis.
## Based-on manually-assigned PC thresholds (see ifelse statements)

## Packages.
library(tidyverse) ## fancy data formatting and plots.
library(cowplot) ## combining fancy plots into one figure.


# setwd("~/Dropbox/gwas_paper/popgen/") ## NOT REQUIRED FROM SHELL.



## Population strucutre. IBD results.

pca_dat = read.table("lhm_hc_2017.mds", header = TRUE)
head(pca_dat)

popstruc_plot = ggplot(pca_dat, aes(C1, C2, colour = IID)) +
	geom_point(shape = 1) +
	theme_bw(base_size = 8) + labs(x = "PC1", y = "PC2") +
	theme(
	legend.position = "none",
		axis.line.x = element_line(size = .3, colour = "black"),
	axis.line.y = element_line(size = .3, colour = "black"),
		plot.margin = unit(c(1,1,1,1.5), "lines"))


## Generate a 'within' file which defines population sub-groups for Plink Fst analysis.
## Non-automated way to define sub-populations.
subpops =	cbind(pca_dat %>% 
                  mutate(group = ifelse(C1 < -0.1, "A",
                                        ifelse(C1 > 0.05, "D", 
                                               ifelse(C2 < -.15, "C", "B")))))
head(subpops)
subpops

write.table(subpops %>% select(FID, IID, group), "lhm_subpops_plink.txt", quote = FALSE, row.names = FALSE, sep = "\t")



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