Commit 9aee6225 authored by Kosmas Hench's avatar Kosmas Hench

update figures and documentation

parent 5b2f663f
......@@ -14,15 +14,15 @@ GeomCustom <- ggproto(
data <- ggproto_parent(Geom, self)$setup_data(data, params)
data
},
draw_group = function(data, panel_scales, coord) {
vp <- grid::viewport(x=data$x, y=data$y)
g <- grid::editGrob(data$grob[[1]], vp=vp)
ggplot2:::ggname("geom_custom", g)
},
required_aes = c("grob","x","y")
)
geom_custom <- function(mapping = NULL,
......@@ -49,64 +49,64 @@ geom_custom <- function(mapping = NULL,
getPofZ <- function(runname,pops){
colN <- c(paste(c(pops[1],pops[1],pops[2],pops[2]),
c('Pure','BC','Pure','BC'),sep = ''))
NHres <-read.csv(paste('../../2_output/08_popGen/11_newHyb/NH_PofZ/newHyb.',runname,'.txt_PofZ.txt',sep=''),sep='\t',
col.names = c('indNR','IndivName',
colN[1],colN[3],
'F1','F2',colN[2],colN[4]))
NHres$IndivName <- read.csv(paste('../../2_output/08_popGen/11_newHyb/NH_PofZ/newHyb.',runname,'_individuals.txt',sep=''),
header = F,col.names = c('IndivName'))[,1] %>%
as.character()
IND_inventory <- read.csv('../../0_data/0_resources/sample_id.csv',sep='\t',header = F,
col.names = c('ID','IndivName','spec','pop','loc'),
stringsAsFactors = F)
data <- NHres %>% merge(.,IND_inventory,by='IndivName',sort=F) %>%
select(IndivName,spec,pop,one_of(colN[1],colN[3],'F1','F2',colN[2],colN[4])) %>%
gather(key = 'BIN',value = "post_prob",4:9) %>%
mutate(RUN=runname)
return(data)
}
plot_fun <- function(data,loc,sel,gl){
pdat <- data %>%
pdat <- data %>%
filter(grepl(loc,RUN)) %>%
arrange(spec,IndivName) %>%
mutate(x2=paste0(spec,IndivName),x = as.numeric(as.factor(x2)))
clr=c('#34b434','#99d499',
'#000000',rgb(.4,.4,.4),
"#d45500",'#cf7c44',
'#ffd736','#ffeeaa')
labsRAW <- substrAllRight(levels(as.factor(pdat$x2)),3)
XL <- list(bel=expression('18151','18153','18155','18156','18157','18158','18159','18162','18165','18171',
'18185','18187','18152','18154',bold('18161'),'18166','18169','18172','18174',
'18175','18176',bold('18178'),'18179',bold('18180'),'18163','18261',bold('18267'),
bold('18274'),bold('18276'),'19881','20092','20120','20126','20128','20135','20149'),
hon=expression('20599','20600','20601','20602','20603','20604','20605','20606','20607','20608',
'20609','20610','20551','20552','20553',bold('20554'),'20555','20556',
bold('20558'),bold('20559'),'20625',bold('20633'),'20635','20638','20560',
'20609','20610','20551','20552','20553','20554','20555','20556',
bold('20558'),'20559','20625',bold('20633'),'20635','20638','20560',
'20561','20562','20563','20564','20565','20566','20567','20568','20571','20572'),
boc=expression("16_21-30","18418","18424","18428","18436","18901","18902","18903",
"18904","18905","18906","18907","18909","18419","18421","18422",
"18426","18427","18429","18430","18432","18434","18912","18915",
"18917",bold("27678"),"16_31-40","18420","18435","18439","18440","18441",
"18442","18445","18446","18447","18448","18450","18454"))
"18904","18905","18906","18907","18909","18419","18421","18422",
"18426","18427","18429","18430","18432","18434","18912","18915",
"18917",bold("27678"),"16_31-40","18420","18435","18439","18440","18441",
"18442","18445","18446","18447","18448","18450","18454"))
groblist <- tibble(RUN=gsub('XX',loc,c('nigXX-pueXX','nigXX-uniXX','pueXX-uniXX')),
grob = gl )
dotdat <- pdat %>%
group_by(IndivName) %>%
dotdat <- pdat %>%
group_by(IndivName) %>%
summarise(spec=spec[1],
x=x[1],
RUN=gsub('XX',loc,c('pueXX-uniXX')))
Grobx=max(pdat$x)+1
p <- ggplot(pdat,
aes(x=x,fill=BIN))+
facet_grid(RUN~.)+
......@@ -130,7 +130,7 @@ plot_fun <- function(data,loc,sel,gl){
legend.position = 'none',
strip.background = element_blank(),
strip.text = element_blank())
return(p)
}
......
......@@ -14,12 +14,12 @@ filter1000snps <- function(name){
mutate(FST_RANK = rank(-WEIR_AND_COCKERHAM_FST,ties.method = "random")) %>%
select(CHROM,POS,WEIR_AND_COCKERHAM_FST,FST_RANK)
filterSet <- raw %>% filter(FST_RANK < 1001)
filterSet <- raw %>% filter(FST_RANK < 801)
h <- hist(filterSet$WEIR_AND_COCKERHAM_FST)
print(min(filterSet$WEIR_AND_COCKERHAM_FST))
write.table(filterSet %>% select(CHROM,POS),paste0(args[7],'2_output/08_popGen/11_newHyb/snpSets/filterSet.1k.',name,'.snps'),
write.table(filterSet %>% select(CHROM,POS),paste0(args[7],'2_output/08_popGen/11_newHyb/snpSets/filterSet.80SNPs.',name,'.snps'),
sep = '\t',row.names = F,quote = F)
return(h)
}
......
......@@ -12,4 +12,4 @@ cd $WORK/2_output/08_popGen/12_fasteprr
mkdir -p step2
mkdir -p step2/XXnameXX
Rscript --vanilla $WORK/0_data/0_scripts/fasteprr_step2_single.R $PWD step1 step2 XXnameXX
Rscript --vanilla $WORK/0_data/0_scripts/fasteprr_step2.R $PWD step1 step2 XXnameXX
......@@ -21,23 +21,23 @@ vcftools \
--gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_$pop1.pop \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_$pop2.pop \
--thin 3000 \
--thin 5000 \
--out $WORK/2_output/08_popGen/11_newHyb/NHrecode/newHyb.$k \
--positions $WORK/2_output/08_popGen/11_newHyb/snpSets/filterSet.1k.$k.snps \
--positions $WORK/2_output/08_popGen/11_newHyb/snpSets/filterSet.80SNPs.$k.snps \
--recode
$WORK/0_data/0_scripts/subsetVcf newHyb.$k.recode.vcf 120
$WORK/0_data/0_scripts/subsetVcf newHyb.$k.recode.vcf 80
mkdir -p NH120
mkdir -p NH80
grep '#CHROM' $WORK/2_output/08_popGen/11_newHyb/NHsubset/newHyb.$k.recode.subset-120.vcf | \
grep '#CHROM' $WORK/2_output/08_popGen/11_newHyb/NHsubset/newHyb.$k.recode.subset-80.vcf | \
cut -f 10- | \
sed 's/\t/\n/g' >$WORK/2_output/08_popGen/11_newHyb/ NH120/newHyb.$k"_individuals.txt"
sed 's/\t/\n/g' >$WORK/2_output/08_popGen/11_newHyb/ NH80/newHyb.$k"_individuals.txt"
java -Xmx1024m -Xms512m -jar $SFTWR/PGDSpider_2.1.1.5/PGDSpider2-cli.jar \
-inputfile $WORK/2_output/08_popGen/11_newHyb/NHsubset/newHyb.$k.recode.subset-120.vcf \
-inputfile $WORK/2_output/08_popGen/11_newHyb/NHsubset/newHyb.$k.recode.subset-80.vcf \
-inputformat VCF \
-outputfile $WORK/2_output/08_popGen/11_newHyb/NH120/newHyb.$k.txt \
-outputfile $WORK/2_output/08_popGen/11_newHyb/NH80/newHyb.$k.txt \
-outputformat NEWHYBRIDS \
-spid $WORK/0_data/0_resources/vcf2nh.spid
......
......@@ -9,10 +9,10 @@
cd $WORK/2_output/08_popGen/11_newHyb
Rscript --vanilla $WORK/0_data/0_scripts/run_parallelNewHybrids.R NH120 $PWD $SFTWR
Rscript --vanilla $WORK/0_data/0_scripts/run_parallelNewHybrids.R NH80 $PWD $SFTWR
mkdir -p NH_PofZ
for k in nigbel-puebel nigbel-unibel puebel-unibel nigboc-pueboc nigboc-uniboc pueboc-uniboc nighon-puehon nighon-unihon puehon-unihon;do
cp N120/NH.Results/newHyb.$k.txt_Results/newHyb.$k.txt_PofZ.txt NH_PofZ/newHyb.$k.txt_PofZ.txt
cp N120/NH.Results/newHyb.$k.txt_Results/newHyb.$k"_individuals.txt" NH_PofZ/newHyb.$k"_individuals.txt"
cp N80/NH.Results/newHyb.80SNPs.$k.txt_Results/newHyb.$k.txt_PofZ.txt NH_PofZ/newHyb.$k.txt_PofZ.txt
cp N80/NH.Results/newHyb.80SNPs.$k.txt_Results/newHyb.$k"_individuals.txt" NH_PofZ/newHyb.$k"_individuals.txt"
done
\ No newline at end of file
......@@ -3,7 +3,7 @@ output: html_document
editor_options:
chunk_output_type: console
---
# Figure 4
# Figure 3
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
......
......@@ -9,19 +9,20 @@ library(ggrepel)
source('../../0_data/0_scripts/S02.functions.R')
run_list <- c("nigbel-puebel","nigbel-unibel","puebel-unibel",
"nigboc-pueboc","nigboc-uniboc","pueboc-uniboc",
"nighon-puehon","nighon-unihon","puehon-unihon")
"nigboc-pueboc","nigboc-uniboc","pueboc-uniboc",
"nighon-puehon","nighon-unihon","puehon-unihon")
pops_list <- data.frame(pop1=c('pue','nig','pue','pue','nig','uni','nig','uni','uni'),
pop2=c('nig','uni','uni','nig','uni','pue','pue','nig','pue'),
pops_list <- data.frame(pop1=c('pue','nig','pue','nig','uni','pue','pue','nig','pue'),
pop2=c('nig','uni','uni','pue','nig','uni','nig','uni','uni'),
stringsAsFactors = F)
pop_order <- c(-1,-1,1,1,-1,-1,1,1,-1)
data <- data.frame(IndivName=c(), spec=c(),pop=c(), BIN=c(),
post_prob=c(),RUN=c(),stringsAsFactors = F)
for (k in 1:9) {
data <- rbind(data,getPofZ(runname = run_list[k],
pops = pops_list[k,]))
pops = pops_list[k,order((1:2)*pop_order[k])]))
}
data$BIN <- factor(as.character(data$BIN),
......@@ -49,7 +50,7 @@ honPCA <- readRDS('../../2_output/08_popGen/04_pca/hon/honpca.Rds')
panPCA <- readRDS('../../2_output/08_popGen/04_pca/pan/panpca.Rds')
introgression_candidates <- c('18161','18178','18180','18267','18274','18276',
'20554','20558','20559','20633','27678')
'20558','20633','27678')
p1 <- plotPCA(pca=belPCA, dataAll=dataAll, locIN='bel',
introgression_candidates=introgression_candidates, clr=clr, fll=fll)
......
#!/usr/bin/env Rscript
#!/usr/bin/env Rscriptlibrary(ggplot2)
library(tidyverse)
library(scales)
library(cowplot)
library(grid)
library(gridSVG)
library(grImport2)
library(grConvert)
library(tidyverse)
library(cowplot)
library(hrbrthemes)
karyo <- read.csv('../../0_data/0_resources/F2.karyo.txt',sep='\t') %>%
mutate(GSTART=lag(cumsum(END),n = 1,default = 0),
GEND=GSTART+END,GROUP=rep(letters[1:2],12)) %>%
select(CHROM,GSTART,GEND,GROUP)
# global -------------
bh <- read.csv('../../2_output/08_popGen/05_fst/bel-hon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BH');
bp <- read.csv('../../2_output/08_popGen/05_fst/bel-hon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BP');
hp <- read.csv('../../2_output/08_popGen/05_fst/hon-boc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HP');
# NIG -------------
bhN <- read.csv('../../2_output/08_popGen/05_fst/nigbel-nighon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BHN');
bpN <- read.csv('../../2_output/08_popGen/05_fst/nigbel-nigboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BPN');
hpN <- read.csv('../../2_output/08_popGen/05_fst/nighon-nigboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HPN');
# PUE -------------
bhP <- read.csv('../../2_output/08_popGen/05_fst/puebel-puehon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BHP');
bpP <- read.csv('../../2_output/08_popGen/05_fst/puebel-pueboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BPP');
hpP <- read.csv('../../2_output/08_popGen/05_fst/puehon-pueboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HPP');
# UNI -------------
bhU <- read.csv('../../2_output/08_popGen/05_fst/unibel-unihon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BHU');
bpU <- read.csv('../../2_output/08_popGen/05_fst/unibel-uniboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BPU');
hpU <- read.csv('../../2_output/08_popGen/05_fst/unihon-uniboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HPU');
# ------------------------------------------
data <- rbind(bh,bp,hp,bhN,bpN,hpN,bhP,bpP,hpP,bhU,bpU,hpU)
data$RUN <- factor(as.character(data$RUN),
levels=c('BH','BP','HP','BHN','BPN','HPN','BHP','BPP','HPP','BHU','BPU','HPU'))
data$COL <- factor(as.character(data$COL),levels=c('BH','BP','HP'))
data2 <- read.csv('../../0_data/0_resources/fst_outlier_windows.txt', sep='\t') %>%
select(-RUN) %>%
group_by(musk) %>%
summarise_all(function(x){x[1]})
musks <- read.csv('../../0_data/0_resources/fst_outlier_ID.txt', sep='\t') %>%
select(-RUN) %>%
group_by(musk) %>%
summarise_all(function(x){x[1]}) %>%
mutate(RUN=factor("BH",levels=c('BH','BP','HP','BHN','BPN','HPN','BHP','BPP','HPP','BHU','BPU','HPU')))
clr <- c('#762a89','#1b7837','#999999')
annoclr <- rgb(.7,.7,.7)
lgclr <- rgb(.9,.9,.9)
yl <- expression(italic("F"[ST]))
yT <- seq(0,.8,length.out = 5)
p1 <- ggplot()+
geom_rect(inherit.aes = F,data=karyo,aes(xmin=GSTART,xmax=GEND,
ymin=-Inf,ymax=Inf,fill=GROUP))+
geom_vline(data=data2,aes(xintercept=xmean),col=annoclr,lwd=.2)+
geom_text(data=musks,aes(x=muskX,y=.7,label=musk))+
geom_point(data=data,aes(x=GPOS,y=WEIGHTED_FST,col=COL),size=.01)+
scale_y_continuous(name = yl,limits=c(-.03,.83),
breaks=yT,labels = c(yT[1],'',yT[3],'',yT[5]))+
scale_x_continuous(expand = c(0,0),breaks = (karyo$GSTART+karyo$GEND)/2,labels = 1:24,position = "top")+
scale_color_manual(name='comparison',values=clr)+
scale_fill_manual(values = c(NA,lgclr),guide=F)+
facet_grid(RUN~.)+
theme_ipsum(grid = F,
axis_title_just = 'c',
axis_title_family = 'bold',
strip_text_size = 9)+
theme(plot.margin = unit(c(2,25,40,0),'pt'),
panel.spacing.y = unit(3,'pt'),
legend.position = 'none',
axis.title.x = element_blank(),
axis.ticks.y = element_line(color = annoclr),
axis.line.y = element_line(color = annoclr),
axis.text.y = element_text(margin = unit(c(0,3,0,0),'pt')),
axis.title.y = element_text(vjust = -1),
strip.background = element_blank(),
strip.text = element_blank())
legGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend-adaptation-cairo.svg"))))
carGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/caribbean-cairo.svg"))))
nigGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/rot-nig-cairo.svg"))))
pueGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/rot-pue-cairo.svg"))))
uniGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/rot-uni-cairo.svg"))))
ys <- .0765
ysc <- .0057
yd <- .2195
boxes = data.frame(x=rep(.9675,4),y=c(ys,ys+yd+ysc,ys+(yd+ysc)*2,ys+(yd+ysc)*3))
lShift <- 0.0225
S09 <- ggdraw(p1)+
geom_rect(data = boxes, aes(xmin = x, xmax = x + .02, ymin = y, ymax = y + yd),
colour = NA, fill = c(rep(lgclr,3),annoclr))+
draw_label('Global', x = .9775, y = boxes$y[4]+.08,
size = 9, angle=-90,col='white')+
draw_label(expression(italic("H. nigricans")), x = .9775, y = boxes$y[3]+.08,
size = 9, angle=-90)+
draw_label(expression(italic("H. puella")), x = .9775, y = boxes$y[2]+.08,
size = 9, angle=-90)+
draw_label(expression(italic("H. unicolor")), x = .9775, y = boxes$y[1]+.08,
size = 9, angle=-90)+
draw_grob(legGrob, 0.01, 0, 1, 0.08)+
draw_grob(carGrob, 0.954, boxes$y[4]+.78*yd, .045, .045)+
draw_grob(nigGrob, 0.954-lShift, boxes$y[3]+.78*yd, .09, .045)+
draw_grob(pueGrob, 0.954-lShift, boxes$y[2]+.78*yd, .09, .045)+
draw_grob(uniGrob, 0.954-lShift, boxes$y[1]+.78*yd, .09, .045)
ggsave(plot = S09,filename = '../output/S09.png',width = 183,height = 183,dpi = 200,units = 'mm')
source('../../0_data/0_scripts/S09.functions.R')
plts <- list()
for(k in 1:7){
plts[[k]] <- trplot(k)
}
tN <- theme(legend.position = 'none')
pG1 <- ggplotGrob(plts[[1]]+tN)
pG2 <- ggplotGrob(plts[[2]]+tN)
pG3 <- ggplotGrob(plts[[3]]+tN)
pG4 <- ggplotGrob(plts[[4]]+tN)
pG5 <- ggplotGrob(plts[[6]]+tN)
pG6 <- ggplotGrob(plts[[5]]+tN)
pG7 <- ggplotGrob(plts[[7]]+tN)
pGr1 <- editGrob(pG1, vp=viewport(x=0.5, y=0.95, angle=45,width = .7), name="pG1")
pGr2 <- editGrob(pG2, vp=viewport(x=0.25, y=0.54, angle=45,width = .28), name="pG2")
pGr3 <- editGrob(pG3, vp=viewport(x=0.25, y=0.357, angle=45,width = .28), name="pG3")
pGr4 <- editGrob(pG4, vp=viewport(x=0.25, y=0.175, angle=45,width = .28), name="pG4")
pGr5 <- editGrob(pG5, vp=viewport(x=0.75, y=0.54, angle=45,width = .28), name="pG5")
pGr6 <- editGrob(pG6, vp=viewport(x=0.75, y=0.357, angle=45,width = .28), name="pG5")
pGr7 <- editGrob(pG7, vp=viewport(x=0.75, y=0.175, angle=45,width = .28), name="pG5")
pGRAD <- ggplot(data = data.frame(x=c(0,1,-1,0),
y=c(0,15,15,0)),
aes(x,y))+geom_polygon(fill='lightgray')+coord_equal()+theme_void()
nigGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/nigricans-cairo.svg"))))
pueGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/puella-cairo.svg"))))
uniGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/unicolor-cairo.svg"))))
belGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/belize-cairo.svg"))))
honGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/honduras-cairo.svg"))))
panGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/panama-cairo.svg"))))
globGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/caribbean-cairo.svg"))))
leg <- get_legend(plts[[1]])
sY <- .035;
sX <- -.03;
S09 <- ggdraw(pGr1)+
draw_plot(pGRAD,.307,-.018,.16,.53)+
draw_grob(pGr2,0,0,1,1)+
draw_grob(leg,.2,-.08,1,1)+
draw_grob(pGr3,0,0,1,1)+
draw_grob(pGr4,0,0,1,1)+
draw_grob(pGr5,0,0,1,1)+
draw_grob(pGr6,0,0,1,1)+
draw_grob(pGr7,0,0,1,1)+
draw_grob(nigGrob, 0.85, 0.38, 0.12, 0.095)+
draw_grob(pueGrob, 0.85, 0.2, 0.12, 0.095)+
draw_grob(uniGrob, 0.85, 0, 0.12, 0.095)+
draw_grob(panGrob, 0.35+sX+.01, 0.37+sY, 0.12, 0.045)+
draw_grob(belGrob, 0.35+sX+.01, 0.19+sY, 0.12, 0.045)+
draw_grob(honGrob, 0.35+sX, 0+sY, 0.14, 0.045)+
draw_grob(globGrob, 0.7, .6, 0.14, 0.08)+
draw_plot_label(letters[1:3],x = c(rep(0.01,2),.52),
y=c(.99,.57,.57), size = 14)
ggsave(plot = S09,filename = '../output/S09.png',width = 183,height = 230,units = 'mm',dpi = 300)
......@@ -6,46 +6,46 @@ library(grConvert)
library(tidyverse)
library(cowplot)
library(hrbrthemes)
karyo <- read.csv('../../0_data/0_resources/F2.karyo.txt',sep='\t') %>%
karyo <- read.csv('../../0_data/0_resources/F2.karyo.txt',sep='\t') %>%
mutate(GSTART=lag(cumsum(END),n = 1,default = 0),
GEND=GSTART+END,GROUP=rep(letters[1:2],12)) %>%
select(CHROM,GSTART,GEND,GROUP)
# global -------------
bh <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/bel-hon-gemma.lm.50k.5k.smooth',sep='\t') %>%
bh <- read.csv('../../2_output/08_popGen/05_fst/bel-hon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BH');
bp <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/bel-hon-gemma.lm.50k.5k.smooth',sep='\t') %>%
bp <- read.csv('../../2_output/08_popGen/05_fst/bel-hon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BP');
hp <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/hon-boc-gemma.lm.50k.5k.smooth',sep='\t') %>%
hp <- read.csv('../../2_output/08_popGen/05_fst/hon-boc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HP');
# NIG -------------
bhN <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/nigbel-nighon-gemma.lm.50k.5k.smooth',sep='\t') %>%
bhN <- read.csv('../../2_output/08_popGen/05_fst/nigbel-nighon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BHN');
bpN <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/nigbel-nigboc-gemma.lm.50k.5k.smooth',sep='\t') %>%
bpN <- read.csv('../../2_output/08_popGen/05_fst/nigbel-nigboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BPN');
hpN <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/nighon-nigboc-gemma.lm.50k.5k.smooth',sep='\t') %>%
hpN <- read.csv('../../2_output/08_popGen/05_fst/nighon-nigboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HPN');
# PUE -------------
bhP <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/puebel-puehon-gemma.lm.50k.5k.smooth',sep='\t') %>%
bhP <- read.csv('../../2_output/08_popGen/05_fst/puebel-puehon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BHP');
bpP <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/puebel-pueboc-gemma.lm.50k.5k.smooth',sep='\t') %>%
bpP <- read.csv('../../2_output/08_popGen/05_fst/puebel-pueboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BPP');
hpP <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/puehon-pueboc-gemma.lm.50k.5k.smooth',sep='\t') %>%
hpP <- read.csv('../../2_output/08_popGen/05_fst/puehon-pueboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HPP');
# UNI -------------
bhU <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/unibel-unihon-gemma.lm.50k.5k.smooth',sep='\t') %>%
bhU <- read.csv('../../2_output/08_popGen/05_fst/unibel-unihon.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BH',RUN='BHU');
bpU <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/unibel-uniboc-gemma.lm.50k.5k.smooth',sep='\t') %>%
bpU <- read.csv('../../2_output/08_popGen/05_fst/unibel-uniboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='BP',RUN='BPU');
hpU <- read.csv('../../2_output/08_popGen/06_GxP/smoothed/50kb/unihon-uniboc-gemma.lm.50k.5k.smooth',sep='\t') %>%
hpU <- read.csv('../../2_output/08_popGen/05_fst/unihon-uniboc.50kb.5kb.windowed.weir.fst',sep='\t') %>%
merge(.,(karyo %>% select(-GEND,-GROUP)),by='CHROM',allx=T) %>% mutate(POS=(BIN_START+BIN_END)/2,GPOS=POS+GSTART,COL='HP',RUN='HPU');
# ------------------------------------------
......@@ -54,13 +54,13 @@ data$RUN <- factor(as.character(data$RUN),
levels=c('BH','BP','HP','BHN','BPN','HPN','BHP','BPP','HPP','BHU','BPU','HPU'))
data$COL <- factor(as.character(data$COL),levels=c('BH','BP','HP'))
data2 <- read.csv('../../0_data/0_resources/fst_outlier_windows.txt', sep='\t') %>%
select(-RUN) %>%
group_by(musk) %>%
data2 <- read.csv('../../0_data/0_resources/fst_outlier_windows.txt', sep='\t') %>%
select(-RUN) %>%
group_by(musk) %>%
summarise_all(function(x){x[1]})
musks <- read.csv('../../0_data/0_resources/fst_outlier_ID.txt', sep='\t') %>%
select(-RUN) %>%
group_by(musk) %>%
musks <- read.csv('../../0_data/0_resources/fst_outlier_ID.txt', sep='\t') %>%
select(-RUN) %>%
group_by(musk) %>%
summarise_all(function(x){x[1]}) %>%
mutate(RUN=factor("BH",levels=c('BH','BP','HP','BHN','BPN','HPN','BHP','BPP','HPP','BHU','BPU','HPU')))
......@@ -68,15 +68,16 @@ clr <- c('#762a89','#1b7837','#999999')
annoclr <- rgb(.7,.7,.7)
lgclr <- rgb(.9,.9,.9)
yl <- expression(paste('-log'[10],' (',italic('p'),'value)',sep=''))
yl <- expression(italic("F"[ST]))
yT <- seq(0,.8,length.out = 5)
p1 <- ggplot()+
geom_rect(inherit.aes = F,data=karyo,aes(xmin=GSTART,xmax=GEND,
ymin=-Inf,ymax=Inf,fill=GROUP))+
geom_vline(data=data2,aes(xintercept=xmean),col=annoclr,lwd=.2)+
geom_text(data=musks,aes(x=muskX,y=8.5,label=musk))+
geom_point(data=data,aes(x=GPOS,y=avgp_wald,col=COL),size=.01)+
scale_y_continuous(name = yl,limits=c(-.3,11),
breaks=0:4*2.5,labels = c(0,'',5,'',10))+
geom_text(data=musks,aes(x=muskX,y=.7,label=musk))+
geom_point(data=data,aes(x=GPOS,y=WEIGHTED_FST,col=COL),size=.01)+
scale_y_continuous(name = yl,limits=c(-.03,.83),
breaks=yT,labels = c(yT[1],'',yT[3],'',yT[5]))+
scale_x_continuous(expand = c(0,0),breaks = (karyo$GSTART+karyo$GEND)/2,labels = 1:24,position = "top")+
scale_color_manual(name='comparison',values=clr)+
scale_fill_manual(values = c(NA,lgclr),guide=F)+
......@@ -94,8 +95,7 @@ p1 <- ggplot()+
axis.text.y = element_text(margin = unit(c(0,3,0,0),'pt')),
axis.title.y = element_text(vjust = -1),
strip.background = element_blank(),
strip.text = element_blank()
)
strip.text = element_blank())
legGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend-adaptation-cairo.svg"))))
carGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/caribbean-cairo.svg"))))
......@@ -107,7 +107,6 @@ ys <- .0765
ysc <- .0057
yd <- .2195
boxes = data.frame(x=rep(.9675,4),y=c(ys,ys+yd+ysc,ys+(yd+ysc)*2,ys+(yd+ysc)*3))
lShift <- 0.0225
S10 <- ggdraw(p1)+
......@@ -120,7 +119,7 @@ S10 <- ggdraw(p1)+
draw_label(expression(italic("H. puella")), x = .9775, y = boxes$y[2]+.08,
size = 9, angle=-90)+
draw_label(expression(italic("H. unicolor")), x = .9775, y = boxes$y[1]+.08,
size = 9, angle=-90)+
size = 9, angle=-90)+
draw_grob(legGrob, 0.01, 0, 1, 0.08)+
draw_grob(carGrob, 0.954, boxes$y[4]+.78*yd, .045, .045)+
draw_grob(nigGrob, 0.954-lShift, boxes$y[3]+.78*yd, .09, .045)+
......@@ -128,4 +127,3 @@ S10 <- ggdraw(p1)+
draw_grob(uniGrob, 0.954-lShift, boxes$y[1]+.78*yd, .09, .045)
ggsave(plot = S10,filename = '../output/S10.png',width = 183,height = 183,dpi = 200,units = 'mm')
This diff is collapsed.
This diff is collapsed.
......@@ -5,21 +5,22 @@ library(cowplot)
library(grid)
library(gridSVG)
library(grImport2)
source('../../0_data/0_scripts/S13.functions.R')
source('../../0_data/0_scripts/F3.functions.R')
plts <- list()
for(k in 1:3){
plts[[k]] <- trplot((8:10)[k])
plts[[k]] <- trplot((8:10)[k])
}
tN <- theme(legend.position = 'none')
pG8 <- ggplotGrob(plts[[1]]+tN)
pG9 <- ggplotGrob(plts[[2]]+tN)
pG10 <- ggplotGrob(plts[[3]]+tN)
pGr8 <- editGrob(pG8, vp=viewport(x=0.5, y=0.95, angle=45,width = .76), name="pG8")
pGr9 <- editGrob(pG9, vp=viewport(x=0.5, y=0.61, angle=45,width = .7), name="pG9")
pGr10 <- editGrob(pG10, vp=viewport(x=0.5, y=0.29, angle=45,width = .7), name="pG10")
pGr8 <- editGrob(pG8, vp=viewport(x=0.5, y=0.97, angle=45,width = .76), name="pG8")
pGr9 <- editGrob(pG9, vp=viewport(x=0.5, y=0.62, angle=45,width = .7), name="pG9")
pGr10 <- editGrob(pG10, vp=viewport(x=0.5, y=0.3, angle=45,width = .7), name="pG10")
npGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/pue-nig-cairo.svg"))))
puGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/uni-pue-cairo.svg"))))
......@@ -31,10 +32,10 @@ leg <- get_legend(plts[[1]]+theme(legend.text = element_text(size = 5),
S13 <- ggdraw(pGr8)+
draw_grob(pGr9,0,0,1,1)+
draw_grob(leg,-.65,-.02,1,1)+
draw_grob(leg,-.85,.04,1,1)+
draw_grob(pGr10,0,0,1,1)+