#########################################################################################
# Plots histograms of length of IBD segments on chromosome X, 
# densities of lengths of IBD segments on chromosome X that match the Neandertal or 
# Denisovan genome and are shared by a particular population or groups of populations
#########################################################################################

library(ggplot2)
library(doBy)

wd <- paste0("chrX/")
setwd(wd)

source("../1000G_pop_functions.R")

chr <- "chrX"

load(file = paste("DenisovaMatNew.Rda", sep = ""))
load(file = paste("NeanderMatNew.Rda", sep = ""))
load(file = paste("AncDenisovaMatNew.Rda", sep = ""))
load(file = paste("AncNeanderMatNew.Rda", sep = ""))
load(file = paste("PopulationMat.Rda", sep = ""))
load("SubPopMat.Rda")

real1 <- which(NeanderMat[,"IBDlengthCorrected"] > 10)
real2 <- which(NeanderMat[,"#SNVs"] > 7)
#real2 <- which(NeanderMat[,"#SNVs"] > 20)
real <- intersect(real1, real2)

DenisovaMat <- DenisovaMat[real,]
NeanderMat <- NeanderMat[real,]
AncDenisovaMat <- AncDenisovaMat[real,]
AncNeanderMat <- AncNeanderMat[real,]
PopulationMat <- PopulationMat[real,]
SubPopMat <- SubPopMat[real,]

de1 <- which(DenisovaMat[,"AltR"] > 0.15)
de2 <- which(DenisovaMat[,"ancientRatio"] > 0.3)
de <- intersect(de1, de2)

ne1 <- which(NeanderMat[,"AltR"] > 0.15)
ne2 <- which(NeanderMat[,"ancientRatio"] > 0.3)
ne <- intersect(ne1, ne2)
 
load(file = paste("DenisovaMat_wChrNew.Rda",sep = ""))
DenisovaMat <- DenisovaMat[real,]

##### Ancestor correction ####
nanca <- which(AncNeanderMat[,"ancestorRatio"] < 0.3)
danca <- which(AncDenisovaMat[,"ancestorRatio"] < 0.3)

de <- intersect(de, danca)
ne <- intersect(ne, nanca)

###### Populations ######

# population is to 1/5 observed
AFR <- which(PopulationMat[,"AFR"] > 0.2) 
EAS <- which(PopulationMat[,"EAS"] > 0.2) 
SAS <- which(PopulationMat[,"SAS"] > 0.2) 
EUR <- which(PopulationMat[,"EUR"] > 0.2) 

# only this population: exclusively
getoPops(PopulationMat)

# not this population: not observed
getNPops(PopulationMat)

# at least once this population: observed
getJPops(PopulationMat)

######### Subpopulations ############

# population is to 1/5 observed
ACB <- which(SubPopMat[,"ACBr"] > 0.2)
ASW <- which(SubPopMat[,"ASWr"] > 0.2)
BEB <- which(SubPopMat[,"BEBr"] > 0.2)
CDX <- which(SubPopMat[,"CDXr"] > 0.2)
CEU <- which(SubPopMat[,"CEUr"] > 0.2)
CHB <- which(SubPopMat[,"CHBr"] > 0.2)
CHS <- which(SubPopMat[,"CHSr"] > 0.2)
CLM <- which(SubPopMat[,"CLMr"] > 0.2)
ESN <- which(SubPopMat[,"ESNr"] > 0.2)
FIN <- which(SubPopMat[,"FINr"] > 0.2)
GBR <- which(SubPopMat[,"GBRr"] > 0.2)
GIH <- which(SubPopMat[,"GIHr"] > 0.2)
GWD <- which(SubPopMat[,"GWDr"] > 0.2)
IBS <- which(SubPopMat[,"IBSr"] > 0.2)
ITU <- which(SubPopMat[,"ITUr"] > 0.2)
JPT <- which(SubPopMat[,"JPTr"] > 0.2)
KHV <- which(SubPopMat[,"KHVr"] > 0.2)
LWK <- which(SubPopMat[,"LWKr"] > 0.2)
MSL <- which(SubPopMat[,"MSLr"] > 0.2)
MXL <- which(SubPopMat[,"MXLr"] > 0.2)
PEL <- which(SubPopMat[,"PELr"] > 0.2)
PJL <- which(SubPopMat[,"PJLr"] > 0.2)
PUR <- which(SubPopMat[,"PURr"] > 0.2)
STU <- which(SubPopMat[,"STUr"] > 0.2)
TSI <- which(SubPopMat[,"TSIr"] > 0.2)
YRI <- which(SubPopMat[,"YRIr"] > 0.2)

# only this population: exclusively
getoSubPops(SubPopMat)

# not this population: not observed
getNSubPops(SubPopMat)

# neither observed in EAS, SAS nor EUR
NEASEUR <- intersect(NEAS, NEUR)
NSASEASEUR <- intersect(NSAS, NEASEUR)

# at least once this subpopulation: observed
getJSubPops(SubPopMat)

# at least once in AFR without ASW and ACB
jafr <- sort(unique(c(JESN, JGWD, JLWK, JMSL, JYRI)))
oafr <- intersect(jafr, NSASEASEUR)
# only AFR without considering AMR, ASW and ACB
oAFR <- oafr

##### get lengths ####
neL <- NeanderMat[ne,"ancientLength"]
deL <- DenisovaMat[de,"ancientLength"]

neLZ <- which(neL > 10)
neL <- neL[neLZ]
ne <- ne[neLZ]

deLZ <- which(deL > 10)
deL <- deL[deLZ]
de <- de[deLZ]

not_neL <- NeanderMat[-ne,"IBDlengthCorrected"]
not_neLZ <- which(not_neL > 10)
not_neL <- not_neL[not_neLZ]

not_deL <- DenisovaMat[-de,"IBDlengthCorrected"]
not_deLZ <- which(not_deL > 10)
not_deL <- not_deL[not_deLZ]


LLL <- DenisovaMat[,"IBDlengthCorrected"]

################### Denisova ###############################
DdeL <- density(deL, bw = 2000, from = 0, to = 60000)

deEAS <- intersect(EAS, de)
deEASL <- DenisovaMat[deEAS,"ancientLength"]
DdeEASL <- density(deEASL, bw = 2000, from = 0, to = 60000)

deSAS <- intersect(SAS, de)
deSASL <- DenisovaMat[deSAS,"ancientLength"]
DdeSASL <- density(deSASL, bw = 2000, from = 0, to = 60000)

deAFR <- intersect(AFR, de)
deAFRL <- DenisovaMat[deAFR,"ancientLength"]
DdeAFRL <- density(deAFRL, bw = 2000, from = 0, to = 60000)

deEUR <- intersect(EUR, de)
deEURL <- DenisovaMat[deEUR,"ancientLength"]
DdeEURL <- density(deEURL, bw = 2000, from = 0, to = 60000)

################## Neandertal ################################
DneL <- density(neL, bw = 2000, from = 0, to = 60000)

neEAS <- intersect(EAS, ne)
neEASL <- NeanderMat[neEAS,"ancientLength"]
DneEASL <- density(neEASL, bw = 2000, from = 0, to = 60000)

neSAS <- intersect(SAS, ne)
neSASL <- NeanderMat[neSAS,"ancientLength"]
DneSASL <- density(neSASL, bw = 2000, from = 0, to = 60000)

neAFR <- intersect(AFR, ne)
neAFRL <- NeanderMat[neAFR,"ancientLength"]
DneAFRL <- density(neAFRL, bw = 2000, from = 0, to = 60000)

neEUR <- intersect(EUR, ne)
neEURL <- NeanderMat[neEUR,"ancientLength"]
DneEURL <- density(neEURL, bw = 2000, from = 0, to = 60000)

########## Africans-NonAfricans ############
deoAFR <- intersect(oAFR, de)
deoAFRL <- DenisovaMat[deoAFR,"ancientLength"]
DdeoAFRL <- density(deoAFRL, bw = 2000, from = 0, to = 60000)

neoAFR <- intersect(oAFR, ne)
neoAFRL <- NeanderMat[neoAFR,"ancientLength"]
DneoAFRL <- density(neoAFRL, bw = 2000, from = 0, to = 60000)

deNAFR <- intersect(NAFR, de)
deNAFRL <- DenisovaMat[deNAFR,"ancientLength"]
DdeNAFRL <- density(deNAFRL, bw = 2000, from = 0, to = 60000)

neNAFR <- intersect(NAFR, ne)
neNAFRL <- NeanderMat[neNAFR,"ancientLength"]
DneNAFRL <- density(neNAFRL, bw = 2000, from = 0, to = 60000)

#### only AFR vs. obs. AFR #####
deJAFR <- intersect(JAFR, de)
deJAFRL <- DenisovaMat[deJAFR,"ancientLength"]
DdeJAFRL <- density(deJAFRL, bw = 2000, from = 0, to = 60000)

neJAFR <- intersect(JAFR, ne)
neJAFRL <- NeanderMat[neJAFR,"ancientLength"]
DneJAFRL <- density(neJAFRL, bw = 2000, from = 0, to = 60000)

JwooAFR <- setdiff(JAFR,  oAFR)
deJwooAFR <- intersect(JwooAFR, de)
deJwooAFRL <- DenisovaMat[deJwooAFR,"ancientLength"]
DdeJwooAFRL <- density(deJwooAFRL, bw = 2000, from = 0, to = 60000)

neJwooAFR <- intersect(JwooAFR, ne)
neJwooAFRL <- NeanderMat[neJwooAFR,"ancientLength"]
DneJwooAFRL <- density(neJwooAFRL, bw = 2000, from = 0, to = 60000)

#### only AFR vs. OoA #####
OoA <- unique(c(JEAS, JSAS, JEUR))
deOoA <- intersect(OoA, de)
deOoAL <- DenisovaMat[deOoA,"ancientLength"]
DdeOoAL <- density(deOoAL, bw = 2000, from = 0, to = 60000)

neOoA <- intersect(OoA, ne)
neOoAL <- NeanderMat[neOoA,"ancientLength"]
DneOoAL <- density(neOoAL, bw = 2000, from = 0, to = 60000)

################################################
##### Histogramm length ####
################################################

name = paste0("IBDlengthHistogram_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

df <- data.frame(counts = LLL)

ggplot(df, aes(x = counts)) +
    geom_histogram(aes(fill = ..count.. < 10000), colour = "darkgreen", 
                   fill = "coral", binwidth = 1000) +
    scale_y_continuous(breaks = seq(0, 2000, 500), expand = c(0, 0)) +
    scale_x_continuous(limits = c(0, 60000), expand = c(0, 0)) +
    geom_vline(xintercept = c(18000),  linetype = "dashed", size = 0.5) +
    labs(title = paste0("Histogram IBD Lengths Chr. X")) +
    theme(legend.position = "none") +
    coord_cartesian(ylim = c(0, 2500),xlim = c(0, 60000)) +
    ylab("Counts") +
    xlab("IBD Length")

dev.off()

################################################
#### DENISOVA LENGTH per Population ####
################################################

name = paste0("IBDlengthDensityDenisovaPop_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)


df <- data.frame(Length = DdeL$x, Denisova1_African = DdeAFRL$y, Denisova4_Asian = DdeEASL$y, 
                 Denisova3_European = DdeEURL$y, Denisova5_South_Asian = DdeSASL$y)

ggplot(df, aes(x = Length)) +
    geom_line(aes(y = Denisova1_African,colour = "Denisova 1African"), size = 1.2) +
    geom_line(aes(y = Denisova4_Asian,colour = "Denisova 4Asian"), size = 1.2) +
    geom_line(aes(y = Denisova3_European,colour = "Denisova 3European"), size = 1.2) +
    geom_line(aes(y = Denisova5_South_Asian,colour = "Denisova 5South_Asian"), size = 1.2) +
    geom_vline(data = df, aes(xintercept = c(7500), colour = c("Denisova 1African")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(13000), colour = c("Denisova 4Asian")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(14500), colour = c("Denisova 5South_Asian")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(22500), colour = c("Denisova 5South_Asian")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(8000), colour = c("Denisova 3European")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(23000), colour = c("Denisova 3European")), 
               linetype = "dashed", size = 0.7) +
    scale_y_continuous(limits = c(0,1.02*max(df[,2:5])), expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    labs(title = paste0("Denisova-matching IBD Length Distributions per Population")) +
    theme(legend.position = c(.77, .85)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, 
                                           linetype = "solid", colour = "black")) +
    coord_cartesian(ylim  = c(0, 1.03*max(df[,2:5])), xlim = c(0, 61000)) +  
    ylab("Density") +
    xlab("IBD Length")  +
    scale_color_hue(name = "Populations", breaks = c("Denisova 1African", "Denisova 2American", "Denisova 3European", 
                                                  "Denisova 4Asian", "Denisova 5South_Asian"),
                    labels = c("Denisova African", "Denisova American", "Denisova European", "Denisova East Asian",
                               "Denisova South Asian")) +
    theme(legend.text = element_text(size = 16), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.4), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0)))

dev.off()

################################################
#### Neander LENGTH per Population ####
################################################

name = paste0("IBDlengthDensityNeanderPop_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

df <- data.frame(Length = DneL$x, Neandertal1_African = DneAFRL$y, Neandertal4_Asian = DneEASL$y, 
                 Neandertal3_European = DneEURL$y, Neandertal5_South_Asian = DneSASL$y)

ggplot(df, aes(x = Length)) +
    geom_line(aes(y = Neandertal1_African, colour = "Neandertal 1African"), size = 1.2) +
    geom_line(aes(y = Neandertal4_Asian, colour = "Neandertal 4Asian"), size = 1.2) +
    geom_line(aes(y = Neandertal3_European, colour = "Neandertal 3European"), size = 1.2) +
    geom_line(aes(y = Neandertal5_South_Asian, colour = "Neandertal 5South_Asian"), size = 1.2) +
    geom_vline(data = df, aes(xintercept = c(6500), colour = c("Neandertal 1African")), 
               linetype = "dashed",  size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(18000), colour = c("Neandertal 4Asian")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(13800), colour = c("Neandertal 5South_Asian")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(21000), colour = c("Neandertal 5South_Asian")), 
               linetype = "dashed", size = 0.7) +
    geom_vline(data = df, aes(xintercept = c(17000), colour = c("Neandertal 3European")), 
               linetype = "dashed", size = 0.7) +
    scale_y_continuous(limits = c(0, 1.02*max(df[,2:5])), expand = c(0, 0)) +
    scale_x_continuous(expand = c(0, 0)) +
    labs(title = paste0("Neandertal-matching IBD Length Distributions per Population")) +
    theme(legend.position = c(.77, .85)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, linetype = "solid", colour = "black")) +
    coord_cartesian(ylim  = c(0, 1.03*max(df[,2:5])), xlim = c(0, 61000)) +  
    ylab("Density") +
    xlab("IBD Length")  +
    scale_color_hue(name = "Populations", breaks = c("Neandertal 1African", "Neandertal 2American", "Neandertal 3European", 
                                                  "Neandertal 4Asian", "Neandertal 5South_Asian"),
                    labels = c("Neandertal African", "Neandertal American", "Neandertal European", 
                               "Neandertal East Asian", "Neandertal South Asian")) +
     theme(legend.text = element_text(size = 16), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.4), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0)))

dev.off()

################################################
#### Denisova LENGTH Africans-NonAfrica ########
################################################

name = paste0("IBDlengthDensityDenisova_Africans-NonAfricans_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

l1 <- length(deoAFRL)
l2 <- length(deNAFRL)

df <- data.frame(Population = factor( c(rep("Africans", l1), rep("Non-Africans", l2)) ), 
   proportion = c(deoAFRL, deNAFRL))

cdf <- summaryBy(data = df, proportion ~ Population, FUN = mean, na.rm = TRUE)

ggplot(df, aes(x = proportion, fill = Population)) + geom_density(alpha = .3, adjust = 0.7) +
    geom_vline(data = cdf, aes(xintercept = c(5800, 13200), colour = Population),
                linetype = "dashed", size = 1.2) +
    geom_vline(data = cdf, aes(xintercept = c(10500, 23500), colour = Population),
               linetype = "dashed", size = 1.2) +
    labs(title = paste0("Denisova IBD Length Distributions: Africans vs. Non-Africans")) +
    scale_y_continuous(breaks = c(0, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 7e-05)) +
    scale_x_continuous(limits = c(0, 60000)) +
    theme(legend.position = c(.8, .7)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, linetype = "solid", colour = "black")) +
    coord_cartesian(ylim = c(0, 7.5e-05), xlim = c(0, 61000)) +
    theme(legend.text = element_text(size = 18), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.5), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0))) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()

################################################
#### Neandertal LENGTH Africans-NonAfrica ######
################################################

name = paste0("IBDlengthDensityNeandertal_Africans-NonAfricans_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

l1 <- length(neoAFRL)
l2 <- length(neNAFRL)

df <- data.frame(Population = factor( c(rep("Africans", l1), rep("Non-Africans", l2)) ), 
   proportion = c(neoAFRL, neNAFRL))

cdf <- summaryBy(data = df, proportion ~ Population, FUN = mean, na.rm = TRUE)

ggplot(df, aes(x = proportion, fill = Population)) + geom_density(alpha = .3, adjust = 0.7) +
    geom_vline(data = cdf, aes(xintercept = c(7000, 15500), colour = Population),
    linetype = "dashed", size = 1.2) +
    labs(title = paste0("Neandertal IBD Length Distributions: Africans vs. Non-Africans")) +
    scale_y_continuous(breaks = c(0, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05)) +
    scale_x_continuous(limits = c(0, 60000)) +
    theme(legend.position = c(.8, .7)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, linetype = "solid", colour = "black")) +
    coord_cartesian(ylim = c(0, 7.5e-05), xlim = c(0, 61000)) +
    theme(legend.text = element_text(size = 18), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.5), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0))) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()

################################################
#### Denisova LENGTH only AFR vs. obs. AFR #####
################################################

name = paste0("IBDlengthDensityDenisova_exclAfricans-obsAfricanswoexcl_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

l1 <- length(deoAFRL)
l2 <- length(deJwooAFRL)

df <- data.frame(Population = factor( c(rep("1only Africans", l1), rep("2Africans", l2))),
                 proportion = c(deoAFRL, deJwooAFRL))

cdf <- summaryBy(data = df, proportion ~ Population, FUN = mean, na.rm = TRUE)

ggplot(df, aes(x = proportion, fill = Population)) + geom_density(alpha = .3, adjust = 0.7) +
    scale_fill_discrete(name = "Population", breaks = c("1only Africans", "2Africans"),
                        labels = c("only Africans", "Africans")) +
    geom_vline(data = cdf, aes(xintercept = c(5800, 9000), colour = Population),
               linetype = "dashed", size = 1.2) +
    geom_vline(data = cdf, aes(xintercept = c(10500,21000), colour = Population),
               linetype = "dashed", size = 1.2) +
    scale_colour_discrete(name = "Population", breaks = c("1only Africans", "2Africans"),
                          labels = c("only Africans", "Africans")) +
    labs(title = paste0("Denisova IBD Length Distributions: only in Africans vs. observed in Africans")) +
    scale_y_continuous(breaks = c(0, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05), 
                       expand = c(0, 0)) +
    scale_x_continuous(limits = c(0, 61000), expand = c(0, 0)) +
    theme(legend.position = c(.8, .7)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, linetype = "solid", colour = "black")) +
    coord_cartesian(ylim = c(0, 7.7e-05), xlim = c(0, 61000)) +
    theme(legend.text = element_text(size = 18), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.4), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0))) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()

#################################################
#### Neandertal LENGTH only AFR vs. obs. AFR ####
#################################################

name = paste0("IBDlengthDensityNeandertal_exclAfricans-obsAfricanswoexcl_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

l1 <- length(neoAFRL)
l2 <- length(neJwooAFRL)

df <- data.frame(Population = factor( c(rep("1only Africans", l1), rep("2Africans", l2))), 
                 proportion = c(neoAFRL, neJwooAFRL))

cdf <- summaryBy(data = df, proportion ~ Population, FUN = mean, na.rm = TRUE)

ggplot(df, aes(x = proportion, fill = Population)) + geom_density(alpha = .3, adjust = 0.7) +
    scale_fill_discrete(name = "Population", breaks = c("1only Africans", "2Africans"),
                        labels = c("only Africans", "Africans")) +
    geom_vline(data = cdf, aes(xintercept = c(7000, 14000), colour = Population), 
               linetype = "dashed", size = 1.2) +
    scale_colour_discrete(name = "Population", breaks = c("1only Africans", "2Africans"), 
                          labels = c("only Africans", "Africans")) +
    labs(title = paste0("Neandertal IBD Length Distributions: only in Africans vs. observed in Africans")) +
    scale_y_continuous(breaks = c(0, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05), 
                       expand = c(0, 0)) +
    scale_x_continuous(limits = c(0, 61000), expand = c(0, 0)) +
    theme(legend.position = c(.8, .7)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, linetype = "solid", colour = "black")) +
    coord_cartesian(ylim = c(0, 7.5e-05), xlim = c(0, 61000)) +
    theme(legend.text = element_text(size = 18), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.4), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0))) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()

################################################
##### Denisova LENGTH only AFR vs. OoA #########
################################################

name = paste0("IBDlengthDensityDenisova_exclAfricans-OoA_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

l1 <- length(deoAFRL)
l2 <- length(deOoAL)

df <- data.frame(Population = factor( c(rep("only Africans", l1), rep("outside of Africa", l2))), 
                 proportion = c(deoAFRL, deOoAL))

cdf <- summaryBy(data = df, proportion ~ Population, FUN = mean, na.rm = TRUE)

ggplot(df, aes(x = proportion, fill = Population)) + geom_density(alpha = .3, adjust = 0.7) +
    geom_vline(data = cdf, aes(xintercept = c(5800, 11000), colour = Population), 
               linetype = "dashed", size = 1.2) +
    geom_vline(data = cdf, aes(xintercept = c(10500, 22000), colour = Population), 
               linetype = "dashed", size = 1.2) +
    labs(title = paste0("Denisova IBD Length Distributions: only in Africans vs. outside of Africa")) +
    scale_y_continuous(breaks = c(0, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05), 
                       expand = c(0, 0)) +
    scale_x_continuous(limits = c(0, 61000), expand = c(0, 0)) +
    theme(legend.position = c(.8, .7)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, linetype = "solid", colour = "black")) +
    coord_cartesian(ylim = c(0, 7.7e-05), xlim = c(0, 61000)) +
    theme(legend.text = element_text(size = 18), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.4), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0))) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()

################################################
##### Neandertal LENGTH only AFR vs. OoA #######
################################################

name = paste0("IBDlengthDensityNeandertal_exclAfricans-OoA_", chr)
#pdf(paste(name, ".pdf", sep = ""), width = 8)
#xfig(paste(name, ".fig", sep = ""), width = 8)

l1 <- length(neoAFRL)
l2 <- length(neOoAL)

df <- data.frame(Population = factor( c(rep("only Africans", l1), rep("outside of Africa", l2))),
                 proportion = c(neoAFRL, neOoAL))

cdf <- summaryBy(data = df, proportion ~ Population, FUN = mean, na.rm = TRUE)

ggplot(df, aes(x = proportion, fill = Population)) + geom_density(alpha = .3, adjust = 0.7) +
    geom_vline(data = cdf, aes(xintercept = c(7000, 16000), colour = Population), 
               linetype = "dashed", size = 1.2) +
    labs(title = paste0("Neandertal IBD Length Distributions: only in Africans vs. outside of Africa")) +
    scale_y_continuous(breaks = c(0, 1e-05, 2e-05, 3e-05, 4e-05, 5e-05, 6e-05, 7e-05, 8e-05), 
                       expand = c(0, 0)) +
    scale_x_continuous(limits = c(0, 61000), expand = c(0, 0)) +
    theme(legend.position = c(.8, .7)) +
    theme(legend.background = element_rect(fill = "gray90", size = 0.3, linetype = "solid", colour = "black")) +
    coord_cartesian(ylim = c(0, 7.5e-05), xlim = c(0, 61000)) +
    theme(legend.text = element_text(size = 18), legend.title = element_text(size = 18)) +
    theme(plot.title = element_text(size = rel(1.4), vjust = 1), axis.title = element_text(size = rel(1.2))) +
    theme(axis.text = element_text(size = rel(1.0))) +
    ylab("Density") +
    xlab("IBD Length")

dev.off()