MigrationEffectOnAdmixture.Rmd

---
title: "Untitled"
author: "Aaron Wolf"
date: "9/27/2018"
output: html_document
---

```{r}
library(data.table)
library(ggplot2)
library(dplyr)
library(tidyr)
library(grid)
library(gridExtra)
```

```{r Coalescent Calls}

mdls <- list.files('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Tenn/', pattern = 'n1_0.05')


dt.Coalescent <- data.table(NULL)

for (mdl in mdls){
    print(mdl)
    
    infile <- paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Tenn/',mdl)
    
    chr_list <- fread(paste0(infile,'/Tenn.chr_list'), col.names = c('chr'), colClasses = c('character'))
    
    ###########
    bedfile.EAS <- fread(input = paste0("gzcat ",infile,"/TreeCalls/*_ASN.gz"),
                   col.names = c('msp_ID_chr','start','end'))
    bedfile.EAS[,len_bp:=end-start]
#    bedfile.EAS <- bedfile.EAS %>% filter(len_bp>50000)
    
    bedfile.EAS %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[3]]) %>% as.data.table()
    bedfile.EAS.chr <- bedfile.EAS %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/1008/1000000) %>% as.data.table()
    
    bedfile.EAS.chr <- right_join(x=bedfile.EAS.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EAS = mean(bedfile.EAS.chr$admix)
    ci.EAS =  paste0(t.test(x = bedfile.EAS.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EAS.chr$admix)$conf.int[2])
    sd.EAS = sd(bedfile.EAS.chr$admix)
    
    
    ###########
    bedfile.EUR <- fread(paste0('gzcat ',infile,'/TreeCalls/*_EUR.gz'),
                   col.names = c('msp_ID_chr','start','end'))
    bedfile.EUR[,len_bp:=end-start]
#    bedfile.EUR <- bedfile.EUR %>% filter(len_bp>50000)
    
    
    bedfile.EUR %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[3]]) %>% as.data.table()
    bedfile.EUR.chr <- bedfile.EUR %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/1006/1000000) %>% as.data.table()
    
    bedfile.EUR.chr <- right_join(x=bedfile.EUR.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EUR = mean(bedfile.EUR.chr$admix)
    ci.EUR = paste0(t.test(x = bedfile.EUR.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EUR.chr$admix)$conf.int[2])
    sd.EUR = sd(bedfile.EUR.chr$admix)
    
    ###########
    
    bedfile.ALL <- rbind(bedfile.EAS, bedfile.EUR)
    bedfile.ALL.chr <- bedfile.ALL %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/2014/1000000) %>% as.data.table()
    
    bedfile.ALL.chr <- right_join(x=bedfile.ALL.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.ALL = mean(bedfile.ALL.chr$admix)
    ci.ALL = paste0(t.test(x = bedfile.ALL.chr$admix)$conf.int[1],' : ',t.test(x = bedfile.ALL.chr$admix)$conf.int[2])
    sd.ALL = sd(bedfile.ALL.chr$admix)
    
    ###########
    row <- data.table(mdl,admix.ALL,ci.ALL,sd.ALL,admix.EAS,ci.EAS,sd.EAS,admix.EUR,ci.EUR,sd.EUR,(admix.EAS/admix.EUR),t.test(bedfile.EAS.chr$admix, bedfile.EUR.chr$admix)$p.value)

    dt.Coalescent <- rbind(dt.Coalescent, row)
}
dt.Coalescent %>% setnames(c('mdl','admix.ALL','ci.ALL','sd.ALL','admix.EAS','ci.EAS','sd.EAS','admix.EUR','ci.EUR','sd.EUR','EAS:EUR','t.test.pval'))

dt.Coalescent[,c('ci.EUR.lwr', 'ci.EUR.upr') := tstrsplit(ci.EUR, ' : ',fixed=TRUE)]
dt.Coalescent[,c('ci.EAS.lwr', 'ci.EAS.upr') := tstrsplit(ci.EAS, ' : ',fixed=TRUE)]
dt.Coalescent[,c('ci.ALL.lwr', 'ci.ALL.upr') := tstrsplit(ci.ALL, ' : ',fixed=TRUE)]

```

```{r Coalescent data graph}

ggplot() + theme_bw() +
    geom_point(data = dt.Coalescent, aes(x=as.numeric(admix.ALL), y=factor(mdl, levels = rev(mdls))), color='blue', size=3 ) +
    geom_errorbarh(data = dt.Coalescent, aes(x=as.numeric(admix.ALL), xmin=as.numeric(ci.ALL.lwr), 
                                             xmax=as.numeric(ci.ALL.upr), 
                                             y=factor(mdl, levels = rev(mdls))), color='blue', size=1 ) +
    geom_point(data = dt.Coalescent, aes(x=as.numeric(admix.EAS), y=factor(mdl, levels = rev(mdls))), color='orange', size=3 ) +
    geom_errorbarh(data = dt.Coalescent, aes(x=as.numeric(admix.EAS), 
                                             xmin=as.numeric(ci.EAS.lwr), xmax=as.numeric(ci.EAS.upr), 
                                             y=factor(mdl, levels = rev(mdls))), color='orange', size=1 ) +
    geom_point(data = dt.Coalescent, aes(x=as.numeric(admix.EUR), y=factor(mdl, levels = rev(mdls))), color='darkgreen', size=3 ) +
    geom_errorbarh(data = dt.Coalescent, aes(x=as.numeric(admix.EUR), 
                                             xmin=as.numeric(ci.EUR.lwr), xmax=as.numeric(ci.EUR.upr), 
                                             y=factor(mdl, levels = rev(mdls))), color='darkgreen', size=1 ) +
    geom_hline(yintercept = 8.5, linetype='dashed', color='darkgrey', size=2) + 
    geom_hline(yintercept = 16.5, linetype='dashed', color='darkgrey', size=2) +
    xlab(label = 'admixture fraction estimate') +
    ylab(label = 'admixture_migration model') +
    ggtitle(label = 'Admixture Estimate from Coalescent Tree Calls') +
    theme()

```



```{r S* calls}
admix <- 'n1_0.05'
mdls <- list.files(paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Sriram/', admix), pattern = 'n1_0.05')


dt.Sstar <- data.table(NULL)

for (mdl in mdls){
    print(mdl)
    
    infile <- paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Sriram/',admix,'/', mdl)
    
    chr_list <- fread(paste0(infile,'/Sriram.chr_list'), col.names = c('chr'), colClasses = c('character'))
    
    ###########
    bedfile.EAS <- fread(input = paste0("gzcat ",infile,"/bedfiles/ALL.merged_ASN.gz"),
                   col.names = c('msp_ID_chr','start','end'))
    bedfile.EAS[,len_bp:=end-start]
    
    bedfile.EAS %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[3]]) %>% as.data.table()
    bedfile.EAS.chr <- bedfile.EAS %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/1008/1000000) %>% as.data.table()
    
    bedfile.EAS.chr <- right_join(x=bedfile.EAS.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EAS = mean(bedfile.EAS.chr$admix)
    ci.EAS =  paste0(t.test(x = bedfile.EAS.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EAS.chr$admix)$conf.int[2])
    sd.EAS = sd(bedfile.EAS.chr$admix)
    
    
    ###########
    bedfile.EUR <- fread(paste0('gzcat ',infile,'/bedfiles/ALL.merged_EUR.gz'),
                   col.names = c('msp_ID_chr','start','end'))
    bedfile.EUR[,len_bp:=end-start]
    
    bedfile.EUR %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[3]]) %>% as.data.table()
    bedfile.EUR.chr <- bedfile.EUR %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/1006/1000000) %>% as.data.table()
    
    bedfile.EUR.chr <- right_join(x=bedfile.EUR.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EUR = mean(bedfile.EUR.chr$admix)
    ci.EUR = paste0(t.test(x = bedfile.EUR.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EUR.chr$admix)$conf.int[2])
    sd.EUR = sd(bedfile.EUR.chr$admix)
    
    ###########
    
    bedfile.ALL <- rbind(bedfile.EAS, bedfile.EUR)
    bedfile.ALL.chr <- bedfile.ALL %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/2014/1000000) %>% as.data.table()
    
    bedfile.ALL.chr <- right_join(x=bedfile.ALL.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.ALL = mean(bedfile.ALL.chr$admix)
    ci.ALL = paste0(t.test(x = bedfile.ALL.chr$admix)$conf.int[1],' : ',t.test(x = bedfile.ALL.chr$admix)$conf.int[2])
    sd.ALL = sd(bedfile.ALL.chr$admix)
    
    ###########
    row <- data.table(mdl,admix.ALL,ci.ALL,sd.ALL,admix.EAS,ci.EAS,sd.EAS,admix.EUR,ci.EUR,sd.EUR,(admix.EAS/admix.EUR),t.test(bedfile.EAS.chr$admix, bedfile.EUR.chr$admix)$p.value)

    dt.Sstar <- rbind(dt.Sstar, row)
}
dt.Sstar %>% setnames(c('mdl','admix.ALL','ci.ALL','sd.ALL','admix.EAS','ci.EAS','sd.EAS','admix.EUR','ci.EUR','sd.EUR','EAS:EUR','t.test.pval'))

dt.Sstar[,c('ci.EUR.lwr', 'ci.EUR.upr') := tstrsplit(ci.EUR, ' : ',fixed=TRUE)]
dt.Sstar[,c('ci.EAS.lwr', 'ci.EAS.upr') := tstrsplit(ci.EAS, ' : ',fixed=TRUE)]
dt.Sstar[,c('ci.ALL.lwr', 'ci.ALL.upr') := tstrsplit(ci.ALL, ' : ',fixed=TRUE)]

```


```{r S* diploid calls}
admix <- 'n1_0.05'
mdls <- list.files(paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Sriram/', admix), pattern = 'n1_0.05')


dt.Sstar.diploid <- data.table(NULL)

for (mdl in mdls){
    print(mdl)
    
    infile <- paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Sriram/',admix,'/', mdl)
    
    chr_list <- fread(paste0(infile,'/Sriram.chr_list'), col.names = c('chr'), colClasses = c('character'))
    
    ###########
    bedfile.EAS <- fread(input = paste0("cat ",infile,"/bedfiles/ALL.merged_ASN.diploid"),
                   col.names = c('msp_ID_chr','start','end'))
    bedfile.EAS[,len_bp:=end-start]
    
    bedfile.EAS %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[2]]) %>% as.data.table()
    bedfile.EAS.chr <- bedfile.EAS %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/504/1000000) %>% as.data.table()
    
    bedfile.EAS.chr <- right_join(x=bedfile.EAS.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EAS = mean(bedfile.EAS.chr$admix)
    ci.EAS =  paste0(t.test(x = bedfile.EAS.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EAS.chr$admix)$conf.int[2])
    sd.EAS = sd(bedfile.EAS.chr$admix)
    
    
    ###########
    bedfile.EUR <- fread(paste0('cat ',infile,'/bedfiles/ALL.merged_EUR.diploid'),
                   col.names = c('msp_ID_chr','start','end'))
    bedfile.EUR[,len_bp:=end-start]
    
    bedfile.EUR %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[2]]) %>% as.data.table()
    bedfile.EUR.chr <- bedfile.EUR %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/503/1000000) %>% as.data.table()
    
    bedfile.EUR.chr <- right_join(x=bedfile.EUR.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EUR = mean(bedfile.EUR.chr$admix)
    ci.EUR = paste0(t.test(x = bedfile.EUR.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EUR.chr$admix)$conf.int[2])
    sd.EUR = sd(bedfile.EUR.chr$admix)
    
    ###########
    
    bedfile.ALL <- rbind(bedfile.EAS, bedfile.EUR)
    bedfile.ALL.chr <- bedfile.ALL %>% rename(chr=msp_ID_chr) %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/1007/1000000) %>% as.data.table()
    
    bedfile.ALL.chr <- right_join(x=bedfile.ALL.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.ALL = mean(bedfile.ALL.chr$admix)
    ci.ALL = paste0(t.test(x = bedfile.ALL.chr$admix)$conf.int[1],' : ',t.test(x = bedfile.ALL.chr$admix)$conf.int[2])
    sd.ALL = sd(bedfile.ALL.chr$admix)
    
    ###########
    row <- data.table(mdl,admix.ALL,ci.ALL,sd.ALL,admix.EAS,ci.EAS,sd.EAS,admix.EUR,ci.EUR,sd.EUR,(admix.EAS/admix.EUR),t.test(bedfile.EAS.chr$admix, bedfile.EUR.chr$admix)$p.value)

    dt.Sstar.diploid <- rbind(dt.Sstar.diploid, row)
}
dt.Sstar.diploid %>% setnames(c('mdl','admix.ALL','ci.ALL','sd.ALL','admix.EAS','ci.EAS','sd.EAS','admix.EUR','ci.EUR','sd.EUR','EAS:EUR','t.test.pval'))

dt.Sstar.diploid[,c('ci.EUR.lwr', 'ci.EUR.upr') := tstrsplit(ci.EUR, ' : ',fixed=TRUE)]
dt.Sstar.diploid[,c('ci.EAS.lwr', 'ci.EAS.upr') := tstrsplit(ci.EAS, ' : ',fixed=TRUE)]
dt.Sstar.diploid[,c('ci.ALL.lwr', 'ci.ALL.upr') := tstrsplit(ci.ALL, ' : ',fixed=TRUE)]

```


```{r S* calls graph}


ggplot() + theme_bw() +
    geom_point(data = dt.Sstar, aes(x=as.numeric(admix.ALL), y=factor(mdl, levels = rev(mdls))), color='blue', size=3 ) +
    geom_errorbarh(data = dt.Sstar, aes(x=as.numeric(admix.ALL), xmin=as.numeric(ci.ALL.lwr), 
                                             xmax=as.numeric(ci.ALL.upr), 
                                             y=factor(mdl, levels = rev(mdls))), color='blue', size=1 ) +
    geom_point(data = dt.Sstar, aes(x=as.numeric(admix.EAS), y=factor(mdl, levels = rev(mdls))), color='orange', size=3 ) +
    geom_errorbarh(data = dt.Sstar, aes(x=as.numeric(admix.EAS), 
                                             xmin=as.numeric(ci.EAS.lwr), xmax=as.numeric(ci.EAS.upr), 
                                             y=factor(mdl, levels = rev(mdls))), color='orange', size=1 ) +
    geom_point(data = dt.Sstar, aes(x=as.numeric(admix.EUR), y=factor(mdl, levels = rev(mdls))), color='darkgreen', size=3 ) +
    geom_errorbarh(data = dt.Sstar, aes(x=as.numeric(admix.EUR), 
                                             xmin=as.numeric(ci.EUR.lwr), xmax=as.numeric(ci.EUR.upr), 
                                             y=factor(mdl, levels = rev(mdls))), color='darkgreen', size=1 ) +
    geom_hline(yintercept = 8.5, linetype='dashed', color='darkgrey', size=2) + 
    geom_hline(yintercept = 16.5, linetype='dashed', color='darkgrey', size=2) +
    xlab(label = 'admixture fraction estimate') +
    ylab(label = 'admixture_migration model') +
    ggtitle(label = 'Admixture Estimate from Sstar Calls') +
    theme()

```



```{r IBDmix Calls}

admix <- 'n1_0.05'
mdls <- list.files(paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Sriram/',admix), pattern = 'n1_0.05')


dt.IBDmix <- data.table(NULL)

for (mdl in mdls){
    print(mdl)
    
    infile <- paste0('~/DATALab/SimulatedDemographic/Sstar/chr1_variable_ref/simulations/Sriram/',admix,'/',mdl)
    
    chr_list <- fread(paste0(infile,'/Sriram.chr_list'), col.names = c('chr'), colClasses = c('integer')) %>% separate(col = 'chr', into = c('chr','pool')) %>% rowwise() %>% mutate(chr=as.integer(chr)) %>% as.data.table()
    
    ###########
    bedfile.EAS <- fread(input = paste0("cat ",infile,"/IBDmixCalls/IntroSeg.ASN.D4E0.002.ALL.txt"),
                   col.names = c('Ind','chr','start','end','LOD1','LOD2'))
    bedfile.EAS[,len_bp:=end-start]
    bedfile.EAS <- bedfile.EAS %>% filter(len_bp>=30000)
    
    #bedfile.EAS %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[3]]) %>% as.data.table()
    bedfile.EAS.chr <- bedfile.EAS %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/504/1000000) %>% as.data.table()
    
    bedfile.EAS.chr <- right_join(x=bedfile.EAS.chr, y = chr_list, by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EAS = mean(bedfile.EAS.chr$admix)
    ci.EAS =  paste0(t.test(x = bedfile.EAS.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EAS.chr$admix)$conf.int[2])
    sd.EAS = sd(bedfile.EAS.chr$admix)
    
    
    ###########
    bedfile.EUR <- fread(paste0('cat ',infile,'/IBDmixCalls/IntroSeg.EUR.D4E0.002.ALL.txt'),
                   col.names = c('Ind','chr','start','end','LOD1','LOD2'))
    bedfile.EUR[,len_bp:=end-start]
    bedfile.EUR <- bedfile.EUR %>% filter(len_bp>=30000)
    
    #bedfile.EUR %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[3]]) %>% as.data.table()
    bedfile.EUR.chr <- bedfile.EUR %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/503/1000000) %>% as.data.table()
    
    bedfile.EUR.chr <- right_join(x=bedfile.EUR.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.EUR = mean(bedfile.EUR.chr$admix)
    ci.EUR = paste0(t.test(x = bedfile.EUR.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.EUR.chr$admix)$conf.int[2])
    sd.EUR = sd(bedfile.EUR.chr$admix)
    
    
    ###########
    bedfile.AFR <- fread(paste0('cat ',infile,'/IBDmixCalls/IntroSeg.AFR.D4E0.002.ALL.txt'),
                   col.names = c('Ind','chr','start','end','LOD1','LOD2'))
    bedfile.AFR[,len_bp:=end-start]
    bedfile.AFR <- bedfile.AFR %>% filter(len_bp>=30000)
    
    #bedfile.EUR %<>% rowwise() %>% mutate(chr=strsplit(msp_ID_chr,split = '_')[[1]][[3]]) %>% as.data.table()
    bedfile.AFR.chr <- bedfile.AFR %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/108/1000000) %>% as.data.table()
    
    bedfile.AFR.chr <- right_join(x=bedfile.AFR.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.AFR = mean(bedfile.AFR.chr$admix)
    ci.AFR = paste0(t.test(x = bedfile.AFR.chr$admix)$conf.int[1], ' : ',t.test(x = bedfile.AFR.chr$admix)$conf.int[2])
    sd.AFR = sd(bedfile.AFR.chr$admix)
    
    ###########
    
    bedfile.ALL <- rbind(bedfile.EAS, bedfile.EUR, bedfile.AFR)
    bedfile.ALL.chr <- bedfile.ALL %>% group_by(chr) %>% summarise(admix = sum(as.numeric(len_bp))/1115/1000000) %>% as.data.table()
    
    bedfile.ALL.chr <- right_join(x=bedfile.ALL.chr, y = chr_list,by = c('chr')) %>% replace_na(list(admix=0)) %>% as.data.table()
    
    admix.ALL = mean(bedfile.ALL.chr$admix)
    ci.ALL = paste0(t.test(x = bedfile.ALL.chr$admix)$conf.int[1],' : ',t.test(x = bedfile.ALL.chr$admix)$conf.int[2])
    sd.ALL = sd(bedfile.ALL.chr$admix)
    
    ###########
    row <- data.table(mdl,admix.ALL,ci.ALL,sd.ALL,admix.EAS,ci.EAS,sd.EAS,admix.EUR,ci.EUR,sd.EUR,admix.AFR,ci.AFR,sd.AFR, (admix.EAS/admix.EUR),t.test(bedfile.EAS.chr$admix, bedfile.EUR.chr$admix)$p.value)

    dt.IBDmix <- rbind(dt.IBDmix, row)
}
dt.IBDmix %>% setnames(c('mdl','admix.ALL','ci.ALL','sd.ALL','admix.EAS','ci.EAS','sd.EAS','admix.EUR','ci.EUR','sd.EUR','admix.AFR','ci.AFR','sd.AFR','EAS:EUR','t.test.pval'))

dt.IBDmix[,c('ci.EUR.lwr', 'ci.EUR.upr') := tstrsplit(ci.EUR, ' : ',fixed=TRUE)]
dt.IBDmix[,c('ci.EAS.lwr', 'ci.EAS.upr') := tstrsplit(ci.EAS, ' : ',fixed=TRUE)]
dt.IBDmix[,c('ci.AFR.lwr', 'ci.AFR.upr') := tstrsplit(ci.AFR, ' : ',fixed=TRUE)]
dt.IBDmix[,c('ci.ALL.lwr', 'ci.ALL.upr') := tstrsplit(ci.ALL, ' : ',fixed=TRUE)]


```



```{r Combined graph}

ggplot() + theme_bw() +
  
    geom_point(data = dt.Sstar, aes(x=as.numeric(admix.ALL), y=factor(mdl, levels = rev(mdls)), color='ALL'), size=3 ) +
    geom_errorbarh(data = dt.Sstar, aes(x=as.numeric(admix.ALL), xmin=as.numeric(ci.ALL.lwr), 
                                             xmax=as.numeric(ci.ALL.upr), 
                                             y=factor(mdl, levels = rev(mdls)), color='ALL'), size=1 ) +
  
    geom_point(data = dt.Sstar, aes(x=as.numeric(admix.EAS), y=factor(mdl, levels = rev(mdls)), color='EAS'), size=3 ) +
    geom_errorbarh(data = dt.Sstar, aes(x=as.numeric(admix.EAS), 
                                             xmin=as.numeric(ci.EAS.lwr), xmax=as.numeric(ci.EAS.upr), 
                                             y=factor(mdl, levels = rev(mdls)), color='EAS'), size=1 ) +
  
    geom_point(data = dt.Sstar, aes(x=as.numeric(admix.EUR), y=factor(mdl, levels = rev(mdls)), color='EUR'), size=3 ) +
    geom_errorbarh(data = dt.Sstar, aes(x=as.numeric(admix.EUR), 
                                             xmin=as.numeric(ci.EUR.lwr), xmax=as.numeric(ci.EUR.upr), 
                                             y=factor(mdl, levels = rev(mdls)), color='EUR'), size=1 ) +
  ###############################  
    geom_point(data = dt.Coalescent, aes(x=as.numeric(admix.ALL), y=factor(mdl, levels = rev(mdls)), color='ALL'), size=3 ) +
    geom_errorbarh(data = dt.Coalescent, aes(x=as.numeric(admix.ALL), xmin=as.numeric(ci.ALL.lwr), 
                                             xmax=as.numeric(ci.ALL.upr), 
                                             y=factor(mdl, levels = rev(mdls)), color='ALL'), size=1 ) +
  
    geom_point(data = dt.Coalescent, aes(x=as.numeric(admix.EAS), y=factor(mdl, levels = rev(mdls)), color='EAS'), size=3 ) +
    geom_errorbarh(data = dt.Coalescent, aes(x=as.numeric(admix.EAS), 
                                             xmin=as.numeric(ci.EAS.lwr), xmax=as.numeric(ci.EAS.upr), 
                                             y=factor(mdl, levels = rev(mdls)), color='EAS'), size=1 ) +
  
    geom_point(data = dt.Coalescent, aes(x=as.numeric(admix.EUR), y=factor(mdl, levels = rev(mdls)), color='EUR'), size=3 ) +
    geom_errorbarh(data = dt.Coalescent, aes(x=as.numeric(admix.EUR), 
                                             xmin=as.numeric(ci.EUR.lwr), xmax=as.numeric(ci.EUR.upr), 
                                             y=factor(mdl, levels = rev(mdls)), color='EUR'), size=1 ) +
  
  ################################
    geom_hline(yintercept = 8.5, linetype='dashed', color='darkgrey', size=2) + 
    geom_hline(yintercept = 16.5, linetype='dashed', color='darkgrey', size=2) +
    xlab(label = 'admixture fraction estimate') +
    ylab(label = 'admixture_migration model') +
    ggtitle(label = 'Admixture Estimate from Sstar and Coalescent Tree Calls') +
    scale_color_manual("",values = c('blue','darkgreen','orange'), breaks=c('ALL','EAS','EUR')) +
    theme()

```


```{r Recent migration graph}



modern.mdls <- mdls[c(1:9)]
ancient.mdls <- mdls[c(1,10:17)]
modern.mig <- data.table(mig=c(0.0, 0.00000001, 0.0000001, 0.000001, 0.00001, 0.000025, 0.00005, 0.0001, 0.0005))
ancient.mig <- data.table(mig=c(0.0, 0.00000001, 0.0000001, 0.000001, 0.00001, 0.00005, 0.0001, 0.00015, 0.0005))

# modern.mdls <- mdls[c(1:8)]
# ancient.mdls <- mdls[c(1,9:15)]
# modern.mig <- data.table(mig=c(0.0, 0.0000001, 0.000001, 0.00001, 0.000025, 0.00005, 0.0001, 0.0005))
# ancient.mig <- data.table(mig=c(0.0, 0.0000001, 0.000001, 0.00001, 0.00005, 0.0001, 0.00015, 0.0005))

dt.Sstar.modern <- cbind(dt.Sstar[mdl %in% modern.mdls], modern.mig)
#dt.Coalescent.modern <- cbind(dt.Coalescent[mdl %in% modern.mdls], modern.mig)
dt.IBDmix.modern <- cbind(dt.IBDmix[mdl %in% modern.mdls], modern.mig)
dt.Sstar.diploid.modern <- cbind(dt.Sstar.diploid[mdl %in% modern.mdls], modern.mig)


dt.Sstar.ancient <- cbind(dt.Sstar[mdl %in% ancient.mdls], ancient.mig)
#dt.Coalescent.ancient <- cbind(dt.Coalescent[mdl %in% ancient.mdls], ancient.mig)
dt.IBDmix.ancient <- cbind(dt.IBDmix[mdl %in% ancient.mdls], modern.mig)
dt.Sstar.diploid.ancient <- cbind(dt.Sstar.diploid[mdl %in% ancient.mdls], ancient.mig)

################################
################################
plot.modern_migration <- 
  ggplot() + theme_bw() +
  
    # geom_point(data = dt.Sstar.modern, aes(y=as.numeric(admix.ALL), x=mdl, color='ALL'), size=3 ) +
    # geom_errorbar(data = dt.Sstar.modern, aes(y=as.numeric(admix.ALL), 
    #                                          ymin=as.numeric(ci.ALL.lwr), ymax=as.numeric(ci.ALL.upr), 
    #                                          x=mdl, color='ALL'), size=1 ) +
    ##############################
    geom_point(data = dt.IBDmix.modern, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.IBDmix'), size=6 ) +
    geom_errorbar(data = dt.IBDmix.modern, aes(y=as.numeric(admix.EAS),
                                             ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr),
                                             x=mdl, color='EAS.IBDmix'), size=1.5 ) +
    geom_line(data = dt.IBDmix.modern, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.IBDmix', group=1), size=1 ) +
  
    geom_point(data = dt.IBDmix.modern, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.IBDmix'), size=6 ) +
    geom_errorbar(data = dt.IBDmix.modern, aes(y=as.numeric(admix.EUR), 
                                             ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
                                             x=mdl, color='EUR.IBDmix'), size=1.5 ) +
    geom_line(data = dt.IBDmix.modern, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.IBDmix', group=1), size=1 ) +
  
    geom_point(data = dt.IBDmix.modern, aes(y=as.numeric(admix.AFR), x=mdl, color='AFR.IBDmix'), size=6 ) +
    geom_errorbar(data = dt.IBDmix.modern, aes(y=as.numeric(admix.AFR), 
                                             ymin=as.numeric(ci.AFR.lwr), ymax=as.numeric(ci.AFR.upr), 
                                             x=mdl, color='AFR.IBDmix'), size=1.5 ) +
    geom_line(data = dt.IBDmix.modern, aes(y=as.numeric(admix.AFR), x=mdl, color='AFR.IBDmix', group=1), size=1 ) +
  

   ###############################  
    # geom_point(data = dt.Coalescent.modern, aes(y=as.numeric(admix.ALL), x=mdl, color='ALL'), size=3 ) +
    # geom_errorbar(data = dt.Coalescent.modern, aes(y=as.numeric(admix.ALL),
    #                                                ymin=as.numeric(ci.ALL.lwr), ymax=as.numeric(ci.ALL.upr), 
    #                                                x=mdl, color='ALL'), size=1 ) +
  
    # geom_point(data = dt.Coalescent.modern, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS'), size=5 ) +
    # geom_errorbar(data = dt.Coalescent.modern, aes(y=as.numeric(admix.EAS), 
    #                                          ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr), 
    #                                          x=mdl, color='EAS'), size=3 ) +
    # 
    # geom_point(data = dt.Coalescent.modern, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR'), size=5 ) +
    # geom_errorbar(data = dt.Coalescent.modern, aes(y=as.numeric(admix.EUR), 
    #                                          ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
    #                                          x=mdl, color='EUR'), size=3 ) +
  
    geom_point(data = dt.Sstar.diploid.modern, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar.diploid'), size=6 ) +
    geom_errorbar(data = dt.Sstar.diploid.modern, aes(y=as.numeric(admix.EAS),
                                             ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr),
                                             x=mdl, color='EAS.Sstar.diploid'), size=1.5 ) +
    geom_line(data = dt.Sstar.diploid.modern, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar.diploid', group=1), size=1 ) +
  
    geom_point(data = dt.Sstar.diploid.modern, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar.diploid'), size=6 ) +
    geom_errorbar(data = dt.Sstar.diploid.modern, aes(y=as.numeric(admix.EUR), 
                                             ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
                                             x=mdl, color='EUR.Sstar.diploid'), size=1.5 ) +
    geom_line(data = dt.Sstar.diploid.modern, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar.diploid', group=1), size=1 ) +

    ##############################
    ##############################
    geom_point(data = dt.Sstar.modern, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar'), size=6 ) +
    geom_errorbar(data = dt.Sstar.modern, aes(y=as.numeric(admix.EAS),
                                             ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr),
                                             x=mdl, color='EAS.Sstar'), size=1.5 ) +
    geom_line(data = dt.Sstar.modern, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar', group=1), size=1 ) +
  
    geom_point(data = dt.Sstar.modern, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar'), size=6 ) +
    geom_errorbar(data = dt.Sstar.modern, aes(y=as.numeric(admix.EUR), 
                                             ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
                                             x=mdl, color='EUR.Sstar'), size=1.5 ) +
    geom_line(data = dt.Sstar.modern, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar', group=1), size=1 ) +
  
  
    
    geom_hline(yintercept = 0.05, linetype='dashed',color='black') +
  
  ################################
    ylab(label = 'admixture fraction estimate') +
    xlab(label = 'migration rate per generation') +
    ggtitle(label = 'Admixture Estimate: Sstar and IBDmix Calls; Modern Migration') +
#    scale_color_manual(values = c('EUR.IBDmix' = 'turquoise4','EUR.Sstar'='steelblue','EAS.IBDmix'='orangered4','EAS.Sstar'='red1', 'AFR.IBDmix'='orange3')) +
    #scale_color_manual(values = c('EUR.IBDmix' = 'turquoise3','EUR.Sstar'='turquoise4','EAS.IBDmix'='red3','EAS.Sstar'='red1', 'AFR.IBDmix'='orange3')) +
    #scale_color_hue(l=50, c=100) +
    #scale_x_discrete(breaks=dt.Sstar.modern$mdl, labels=c(0.0, 0.0000001, 0.000001, 0.00001, 0.000025, 0.00005, 0.0001, 0.0005)) +
    scale_x_discrete(breaks=dt.Sstar.modern$mdl, labels=c(0.0, 0.00000001, 0.0000001, 0.000001, 0.00001, 0.000025, 0.00005, 0.0001, 0.0005)) + 
    scale_y_continuous(breaks = seq(0,0.16,by = 0.01), labels=seq(0,0.16,by=0.01)) + 
#    scale_y_continuous(trans = "log10", limits = c(NA,-2)) +
#    coord_cartesian(ylim = c(0,0.16), expand = c(0)) +
    coord_cartesian(ylim = c(0,0.08), expand = c(0)) +
    theme(
      plot.title=element_text(face="bold", size = 14),
      axis.text=element_text(face="bold", size = 14),
      axis.title=element_text(face="bold", size = 14),
      legend.title=element_text(face="bold", size = 14),
      legend.text=element_text(face="bold", size = 14)
    )
##################################
##################################


##################################
##################################

plot.ancient_migration <- 
  ggplot() + theme_bw() +
  
    # geom_point(data = dt.Sstar.ancient, aes(y=as.numeric(admix.ALL), x=mdl, color='ALL'), size=3 ) +
    # geom_errorbar(data = dt.Sstar.ancient, aes(y=as.numeric(admix.ALL), 
    #                                          ymin=as.numeric(ci.ALL.lwr), ymax=as.numeric(ci.ALL.upr), 
    #                                          x=mdl, color='ALL'), size=1 ) +

    geom_point(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.IBDmix'), size=6 ) +
    geom_errorbar(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.EAS),
                                             ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr),
                                             x=mdl, color='EAS.IBDmix'), size=1.5 ) +
    geom_line(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.IBDmix', group=1), size=1 ) +
  
    geom_point(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.IBDmix'), size=6 ) +
    geom_errorbar(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.EUR), 
                                             ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
                                             x=mdl, color='EUR.IBDmix'), size=1.5 ) +
    geom_line(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.IBDmix', group=1), size=1 ) +
  
    geom_point(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.AFR), x=mdl, color='AFR.IBDmix'), size=6 ) +
    geom_errorbar(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.AFR), 
                                             ymin=as.numeric(ci.AFR.lwr), ymax=as.numeric(ci.AFR.upr), 
                                             x=mdl, color='AFR.IBDmix'), size=1.5 ) +
    geom_line(data = dt.IBDmix.ancient, aes(y=as.numeric(admix.AFR), x=mdl, color='AFR.IBDmix', group=1), size=1 ) +
  
    geom_hline(yintercept = 0.05, linetype='dashed',color='black') +
  
   ###############################  
    # geom_point(data = dt.Coalescent.ancient, aes(y=as.numeric(admix.ALL), x=mdl, color='ALL'), size=3 ) +
    # geom_errorbar(data = dt.Coalescent.ancient, aes(y=as.numeric(admix.ALL),
    #                                                ymin=as.numeric(ci.ALL.lwr), ymax=as.numeric(ci.ALL.upr), 
    #                                                x=mdl, color='ALL'), size=1 ) +
  
    # geom_point(data = dt.Coalescent.ancient, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS'), size=5 ) +
    # geom_errorbar(data = dt.Coalescent.ancient, aes(y=as.numeric(admix.EAS), 
    #                                          ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr), 
    #                                          x=mdl, color='EAS'), size=3 ) +
    # 
    # geom_point(data = dt.Coalescent.ancient, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR'), size=5 ) +
    # geom_errorbar(data = dt.Coalescent.ancient, aes(y=as.numeric(admix.EUR), 
    #                                          ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
    #                                          x=mdl, color='EUR'), size=3 ) +

    geom_point(data = dt.Sstar.diploid.ancient, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar.diploid'), size=6 ) +
    geom_errorbar(data = dt.Sstar.diploid.ancient, aes(y=as.numeric(admix.EAS),
                                             ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr),
                                             x=mdl, color='EAS.Sstar.diploid'), size=1.5 ) +
    geom_line(data = dt.Sstar.diploid.ancient, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar.diploid', group=1), size=1 ) +
  
    geom_point(data = dt.Sstar.diploid.ancient, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar.diploid'), size=6 ) +
    geom_errorbar(data = dt.Sstar.diploid.ancient, aes(y=as.numeric(admix.EUR), 
                                             ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
                                             x=mdl, color='EUR.Sstar.diploid'), size=1.5 ) +
    geom_line(data = dt.Sstar.diploid.ancient, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar.diploid', group=1), size=1 ) +

    ##############################

    geom_point(data = dt.Sstar.ancient, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar'), size=6 ) +
    geom_errorbar(data = dt.Sstar.ancient, aes(y=as.numeric(admix.EAS),
                                             ymin=as.numeric(ci.EAS.lwr), ymax=as.numeric(ci.EAS.upr),
                                             x=mdl, color='EAS.Sstar'), size=1.5 ) +
    geom_line(data = dt.Sstar.ancient, aes(y=as.numeric(admix.EAS), x=mdl, color='EAS.Sstar', group=1), size=1 ) +
  
    geom_point(data = dt.Sstar.ancient, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar'), size=6 ) +
    geom_errorbar(data = dt.Sstar.ancient, aes(y=as.numeric(admix.EUR), 
                                             ymin=as.numeric(ci.EUR.lwr), ymax=as.numeric(ci.EUR.upr), 
                                             x=mdl, color='EUR.Sstar'), size=1.5 ) +
    geom_line(data = dt.Sstar.ancient, aes(y=as.numeric(admix.EUR), x=mdl, color='EUR.Sstar', group=1), size=1 ) +
  
  ################################
    ylab(label = 'admixture fraction estimate') +
    xlab(label = 'migration rate per generation') +
    ggtitle(label = 'Admixture Estimate: Sstar and IBDmix Calls; Ancient Migration') +
#    scale_color_manual(values = c('EUR.IBDmix' = 'turquoise4','EUR.Sstar'='steelblue','EAS.IBDmix'='orangered4','EAS.Sstar'='red1', 'AFR.IBDmix'='orange3')) +
    #scale_color_manual(values = c('EUR.IBDmix' = 'turquoise3','EUR.Sstar'='turquoise4','EAS.IBDmix'='red3','EAS.Sstar'='red1', 'AFR.IBDmix'='orange3')) +
    #scale_color_hue(l=50, c=100) +
    #scale_x_discrete(breaks=dt.Sstar.ancient$mdl, labels=c(0.0, 0.0000001, 0.000001, 0.00001, 0.00005, 0.0001, 0.00015, 0.0005)) +
    scale_x_discrete(breaks=dt.Sstar.ancient$mdl, labels=c(0.0, 0.00000001, 0.0000001, 0.000001, 0.00001, 0.00005, 0.0001, 0.00015, 0.0005)) + 
    scale_y_continuous(breaks = seq(0,0.16,by = 0.01), labels=seq(0,0.16,by=0.01)) +
#    scale_y_continuous(trans = "log10", limits = c(NA,-2), breaks = seq(0,0.12,by = 0.01), labels=seq(0,0.12,by=0.01)) +
#    coord_cartesian(ylim = c(0,0.16), expand = c(0)) +
    coord_cartesian(ylim = c(0,0.08), expand = c(0)) +
    theme(      
      plot.title=element_text(face="bold", size = 14),
      axis.text=element_text(face="bold", size = 14),
      axis.title=element_text(face="bold", size = 14),
      legend.title=element_text(face="bold", size = 14),
      legend.text=element_text(face="bold", size = 14)
      )
##################################
##################################



grid.arrange(plot.ancient_migration, plot.modern_migration, ncol=2)


```