Commit 89b22387 authored by Kosmas Hench's avatar Kosmas Hench

incl linkage decay

parent eae8412e
......@@ -134,10 +134,10 @@ trplot <- function(sel){
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
scale_fill_gradientn(name=expression(bar(italic(r)^2)),colours=clr,
values=rescale(c(1,.5,.08,.03,.015,.01,0)),
limits=c(0,1),
guide = 'legend',breaks=c(0,.005,.01,.02,.03,.1,1))+
# scale_fill_gradientn(name=expression(bar(italic(r)^2)),colours=clr,
# values=rescale(c(1,.5,.08,.03,.015,.01,0)),
# limits=c(0,1),
# guide = 'legend',breaks=c(0,.005,.01,.02,.03,.1,1))+
theme_void()+
theme(legend.position = c(.9,.75),
legend.direction = 'vertical')
......@@ -151,14 +151,21 @@ trplot <- function(sel){
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
scale_fill_gradientn(name=expression(bar(r^2)),colours=clr,
values=rescale(c(1,.08,.03,.015,.01,0)),
limits=c(0,1),guide = 'legend',breaks=c(0,.005,.01,.02,.03,.1,1))+
#scale_fill_gradientn(name=expression(bar(r^2)),colours=clr,
# values=rescale(c(1,.08,.03,.015,.01,0)),
# limits=c(0,1),guide = 'legend',breaks=c(0,.005,.01,.02,.03,.1,1))+
theme_void()+
theme(legend.position = c(.7,.75),legend.direction = 'vertical')
}
return(p1)
# cs <- scale_fill_viridis_c(name=expression(bar(r^2)),begin = 0,end = 0.55,direction = -1)
cs <- scale_fill_viridis_c(name=expression(bar(r^2)),
limits = c(0,0.55),values = c(0,.08,1),
option = 'B',direction = -1)
# cs <- scale_fill_gradientn(name=expression(bar(r^2)),colours=clr,
# values=rescale(c(1,.08,.03,.015,.01,0)),
# limits=c(0,0.55))
return(p1+cs)
}
run_one_sample_test <- function(pop,n,data){
......
......@@ -188,9 +188,6 @@ trplot <- function(sel){
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
scale_fill_gradientn(name=expression(bar(italic(r)^2)),colours=clr,
values=rescale(c(1,.08,.03,.015,.01,0)),
limits=c(0,1),guide = 'legend',breaks=c(0,.005,.01,.02,.03,.1,1))+
scale_color_manual(values=viridis_pal(option='A')(6))+theme_void()+
theme(legend.position = c(.7,.75),legend.direction = 'vertical')
} else {
......@@ -203,11 +200,11 @@ trplot <- function(sel){
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
scale_fill_gradientn(name=expression(bar(r^2)),colours=clr,
values=rescale(c(1,.08,.03,.015,.01,0)),
limits=c(0,1),guide = 'legend',breaks=c(0,.005,.01,.02,.03,.1,1))+
scale_color_manual(values=viridis_pal(option='A')(6))+theme_void()+
theme(legend.position = c(.7,.75),legend.direction = 'vertical')
}
return(p1)
cs <- scale_fill_viridis_c(name=expression(bar(r^2)),
limits = c(0,0.55),values = c(0,.08,1),
option = 'B',direction = -1)
return(p1+cs)
}
\ No newline at end of file
......@@ -138,7 +138,7 @@ trplot <- function(sel){
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
scale_fill_viridis_c(name=expression(bar(italic(r)^2)))+
scico::scale_fill_scico(name=expression(bar(italic(r)^2)),palette = 'lapaz',direction = -1)+
theme_void()+
theme(legend.position = c(.9,.75),
legend.direction = 'vertical')
......
#PBS -l elapstim_req=0:20:00
#PBS -l memsz_job=10gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N linkage_decay
#PBS -q clexpress
#PBS -o 2.2.8.1.linkage_decay.stdout
#PBS -e 2.2.8.1.linkage_decay.stderr
cd $WORK/2_output/08_popGen/07_LD
mkdir -p decay
cd decay
RANDOM=27678 #set seed
for k in {1..20};do
LG=$(echo $(( (RANDOM % 24) + 1 )) | awk '$1<10{print "0"$1} $1>=10{print$1}')
RN1=$(( ( RANDOM % 13684534 ) + 1))
RN2=$(($RN2+15000))
echo $RN2
vcftools \
--gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--chr LG$LG \
--from-bp $RN1 \
--to-bp $RN2 \
--hap-r2 \
--out dist.$k.$RN1
done
gzip ./*.hap.ld
......@@ -46,7 +46,7 @@ honGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/hond
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]]+theme(legend.direction = 'horizontal')+ guides(fill = guide_legend(nrow=1)))
leg <- get_legend(plts[[1]]+theme(legend.direction = 'horizontal',legend.key.width = unit(60,'pt')))
sY <- .035;
sX <- -.03;
F4 <- ggdraw(pGr1)+
......@@ -54,7 +54,7 @@ F4 <- ggdraw(pGr1)+
draw_plot(boxGenes+theme(legend.position='none'),.65,.54,.45,.21)+
draw_plot(pGRAD,.307,-.018,.16,.53)+
draw_grob(pGr2,0,0,1,1)+
draw_grob(leg,-.33,.22,1,1)+
draw_grob(leg,-.25,.22,1,1)+
draw_grob(pGr3,0,0,1,1)+
draw_grob(pGr4,0,0,1,1)+
draw_grob(pGr5,0,0,1,1)+
......
#!/usr/bin/env Rscript
library(sp)
library(maptools)
library(scales)
library(PBSmapping)
library(RColorBrewer)
library(tidyverse)
library(ggforce)
library(cowplot)
library(scatterpie)
library(hrbrthemes)
library(cowplot)
library(grid)
library(gridSVG)
library(grImport2)
library(grConvert)
library(ggmap)
library(marmap)
source("../../0_data/0_scripts/F1.functions.R")
dataAll <- read.csv('../../0_data/0_resources/F1.sample.txt',sep='\t') %>% mutate(loc=substrRight(as.character(id),3))
data <- dataAll %>% filter(sample=="sample")
dataSum <- data%>% rowwise()%>% mutate(nn = sum(as.numeric(strsplit(as.character(`Latittude.N`),' ')[[1]])*c(1,1/60,1/3600)),
ww = sum(as.numeric(strsplit(as.character(`Longitude.W`),' ')[[1]])*c(1,1/60,1/3600))) %>%
group_by(loc) %>% summarise(n=mean(nn,na.rm = T),w=-mean(ww,na.rm = T)) %>%
as.data.frame() %>% cbind(.,read.csv('../../0_data/0_resources/F1.pointer.csv'))
crvs <- read.csv('../../0_data/0_resources/F1.curve.csv')
xlimW = c(-100,-53); ylimW = c(7,30.8)
worldmap = map_data("world")
names(worldmap) <- c("X","Y","PID","POS","region","subregion")
worldmap = clipPolys(worldmap, xlim=xlimW,ylim=ylimW, keepExtra=TRUE)
compGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/north-cairo.svg"))))
legGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend-map-cairo.svg"))))
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"))))
legPGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend_pairs-cairo.svg"))))
ypos = c(2,1,0)*1.6;
xpos=c(1,0,1)*4.57
sampSize <- data %>% group_by(spec,loc) %>% summarise(n=length(id)) %>% filter(spec!='gum') %>% arrange(loc) %>%
cbind(y=c(rep(24.4,3)+ypos, # belize
rep(13.6,3)+ypos, # panama
rep(9.2,3)+ypos), # honduras
x=c(rep(-94.7,3)+xpos, # belize
rep(-78,3)+xpos, # panama
rep(-98.75,3)+xpos)) # honduras
uS <- readShapeSpatial('../../0_data/0_resources/F1_shps/F1.npu_distribution.shp');
uS_f <- fortify(uS,region="DN");
names(uS_f) <- c("X","Y","POS","hole","piece","id","PID"); uS_f$PID <- as.numeric(uS_f$PID)
uS_f <- clipPolys(uS_f, xlim=xlimW,ylim=ylimW, keepExtra=TRUE)
hS <- readShapeSpatial('../../0_data/0_resources/F1_shps/F1.hypoplectrus_distribution.shp');
hS_f <- fortify(hS,region="DN");
names(hS_f) <- c("X","Y","POS","hole","piece","id","PID"); hS_f$PID <- as.numeric(hS_f$PID)
hS_f <- clipPolys(hS_f, xlim=xlimW,ylim=ylimW, keepExtra=TRUE)
clrShps <- c("#ffc097","#b1d4f8")
cFILL <- rgb(.7,.7,.7)
ccc<-rgb(0,.4,.8)
p1 <- ggplot()+coord_map(projection = 'mercator',xlim=xlimW,ylim=ylimW)+
labs(y="Latitude",x="Longitude") +
theme_mapK +
geom_polygon(data=worldmap,aes(X,Y,group=PID),fill = cFILL,
col=rgb(0,0,0),lwd=.2)+
geom_polygon(data=hS_f,aes(X,Y,group=PID),fill=clrShps[2],col=rgb(.3,.3,.3),lwd=.2)+
geom_polygon(data=uS_f,aes(X,Y,group=PID),fill=clrShps[1],col=clrShps[1],lwd=.2)+
geom_point(data=sampSize,aes(x=x,y=y),col='white',size=7)+
geom_text(data=sampSize,aes(x=x,y=y,label=n),size=3)+
geom_bezier(data = crvs,aes(x= x, y = y, group = loc),col=ccc,lwd=.4) +
geom_point(data = crvs[c(3,6,9),],aes(x= x, y = y, group = loc),size=5,col=ccc,shape=21,fill='white')+
geom_text(data = crvs[c(3,6,9),],aes(x= x, y = y, label=c('B','P','H')),size=3)+
geom_point(data=dataSum,aes(x=w,y=n),size=1.5,shape=21,fill='white')
fstdata <- read.csv('../../2_output/08_popGen/05_fst/genome_wide_weighted_mean_fst.txt',
sep='\t',skip = 1,head=F,col.names = c('pair','fst')) %>%
filter(!pair %in% c('NU','PN','PU')) %>%
arrange(fst) %>% mutate(x=row_number()) %>% rowwise() %>%
mutate(loc = substr(as.character(pair),2,2),
pop1 = substr(as.character(pair),1,1),
pop2 = substr(as.character(pair),3,3),
pairclass=paste(pop1,pop2,sep='-'))
clrP <- c('#fb8620','#1b519c','#d93327')
p5 <- ggplot(fstdata,aes(x=x,y=fst,fill=pairclass))+
geom_bar(stat='identity')+
annotation_custom(legPGrob,xmin =0.15,xmax =5,ymin =.008,ymax = .035)+
scale_fill_manual(values = clrP)+
scale_y_continuous(name=expression(italic("F"['ST'])),
expand = c(0,0))+
scale_x_continuous(expand = c(.01,0),breaks = 1:9,
labels = c('H','H','B','H','P','B','B','P','P'))+
theme_mapK+
theme(legend.position = 'none',
panel.grid.major.y = element_line(color=rgb(.9,.9,.9)),
axis.title.x = element_blank(),
axis.text.x = element_text(#face = 'bold',
size = 8));
clr<-c('#000000','#d45500','#000000')
fll<-c('#000000','#d45500','#ffffff')
eVclr <- c(rep('darkred',2),rep(cFILL,8))
belPCA <- readRDS('../../2_output/08_popGen/04_pca/bel/belpca.Rds')
honPCA <- readRDS('../../2_output/08_popGen/04_pca/hon/honpca.Rds')
panPCA <- readRDS('../../2_output/08_popGen/04_pca/pan/panpca.Rds')
dataHON <- cbind(dataAll %>% filter(sample=='sample',loc=='hon')%>% select(id,spec),
honPCA$scores); names(dataHON)[3:12]<- paste('PC',1:10,sep='')
exp_varHON <- (honPCA$singular.values[1:10])^2/length(honPCA$maf)
xlabHON <- paste('PC1 (',sprintf("%.1f",exp_varHON[1]*100),'%)')# explained varinace)')
ylabHON <- paste('PC2 (',sprintf("%.1f",exp_varHON[2]*100),'%)')# explained varinace)')
p2 <- ggplot(dataHON,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=2,shape=21)+
scale_color_manual(values=clr,guide=F)+
ggtitle('Honduras')+
scale_fill_manual(values=fll,guide=F)+
scale_x_continuous(name=xlabHON)+
scale_y_continuous(name=ylabHON)
dataBEL <- cbind(dataAll %>% filter(sample=='sample',loc=='bel') %>% select(id,spec),
belPCA$scores); names(dataBEL)[3:12]<- paste('PC',1:10,sep='')
exp_varBEL <- (belPCA$singular.values[1:10])^2/length(belPCA$maf)
xlabBEL <- paste('PC1 (',sprintf("%.1f",exp_varBEL[1]*100),'%)')
ylabBEL <- paste('PC2 (',sprintf("%.1f",exp_varBEL[2]*100),'%)')
p3 <- ggplot(dataBEL,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=2,shape=21)+
scale_color_manual(values=clr,guide=F)+
ggtitle('Belize')+
scale_fill_manual(values=fll,guide=F)+
scale_x_continuous(name=xlabBEL)+
scale_y_continuous(name=ylabBEL)
dataPAN <- cbind(dataAll %>% filter(sample=='sample',loc=='boc',spec!='gum')%>% select(id,spec),
panPCA$scores); names(dataPAN)[3:12]<- paste('PC',1:10,sep='')
exp_varPAN <- (panPCA$singular.values[1:10])^2/length(panPCA$maf)
xlabPAN <- paste('PC1 (',sprintf("%.1f",exp_varPAN[1]*100),'%)')
ylabPAN <- paste('PC2 (',sprintf("%.1f",exp_varPAN[2]*100),'%)')
p4 <- ggplot(dataPAN,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=2,shape=21)+
scale_color_manual(values=clr,guide=F)+
ggtitle('Panama')+
scale_fill_manual(values=fll,guide=F)+
scale_x_continuous(name=xlabPAN)+
scale_y_continuous(name=ylabPAN)
library(hypoimg)
ggflag <- function(name='panama'){(ggplot()+hypo_anno_flag(name)+theme_void()) %>% ggplotGrob()}
tm <- theme_bw()+theme(plot.title = element_text(hjust = 0.5,vjust = 4),
panel.grid = element_blank(),
panel.border = element_blank(),
axis.line = element_line(size=.2),plot.margin = unit(c(30,0,0,0),'pt'))
pca <- cowplot::plot_grid(p2+tm,p3+tm,p4+tm,nrow = 1)
plt <- ggdraw(pca)+
draw_grob(ggflag('honduras'),.03,.87,.1,.1)+
draw_grob(ggflag('belize'),.36,.87,.1,.1)+
draw_grob(ggflag(),.69,.87,.1,.1)
species <- c('nigricans','puella','unicolor')
leg <- hypo_legend_single(species = species, color_map = fll,ncol = 3,
circle_color = 'black', plot_names = TRUE,circle_lwd = .2)
plot_grid(plt,NULL,leg,rel_heights = c(4,.3,.97),ncol = 1) %>%
ggsave(.,filename = '~/Desktop/pca.pdf',width = 10,height = 4.8,device = cairo_pdf)
#!/usr/bin/env Rscript
library(grid)
library(gridSVG)
library(grImport2)
library(grConvert)
library(tidyverse)
library(rtracklayer)
library(hrbrthemes)
library(cowplot)
require(rtracklayer)
source('../../0_data/0_scripts/F3.functions.R');
source('../../0_data/0_scripts/F3.plot_fun.R')
p1 <- create_K_plot(searchLG = "LG09",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(17800000,18000000),
searchgene = c("SOX10_1"),secondary_genes = c("RNASEH2A"),searchsnp = c(17871737,17872597,17873443),
muskID = 'A')
p2 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(20080000,20500000),
searchgene = c("Casz1_2"),secondary_genes = c(),searchsnp = c(20316944,20317120,20323661,20323670,20333895,20347263),
muskID = 'B')
p3 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22100000,22350000),
searchgene = c("hoxc13a"),secondary_genes = c('hoxc10a',"hoxc11a","hoxc12a","calcoco1_1"),
searchsnp = c(),
muskID = 'C')
p4 <- create_K_plot(searchLG = "LG17",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22500000,22670000),
searchgene = c('LWS',"SWS2abeta","SWS2aalpha","SWS2b"),searchsnp = c(22553970,22561903,22566254),
secondary_genes = c("Hcfc1","HCFC1_2","HCFC1_1","GNL3L","TFE3_0","MDFIC2_1","CXXC1_3","CXXC1_1",'Mbd1','CCDC120'),
muskID = 'D')
legGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend-pw-cairo.svg"))))
F3 <- plot_grid(NULL,NULL,NULL,NULL,
p1,NULL,p2,NULL,
NULL,NULL,NULL,NULL,
p3,NULL,p4,NULL,
NULL,NULL,NULL,NULL,
ncol=4,rel_heights = c(.03,1,.03,1,.1),rel_widths = c(1,.025,1,.02),
labels = c('a','','b','',
'','','','',
'c','','d','',
'','','','',
'','','',''),label_size = 10)+
draw_grob(legGrob, 0.1, 0, .8, 0.04)
ggsave(plot = F3,filename = '~/Desktop/ranges.pdf',
width = 16/1.8,height = 9/1.8,device = cairo_pdf)
exp_plt <- function(plt,nm){
plot_grid(plt,legGrob,rel_heights = c(1,.1),ncol = 1) %>%
ggsave(plot = .,filename = stringr::str_c('~/Desktop/ranges.',nm,'.pdf'),
width = (16/1.8),height = (9/1.8),device = cairo_pdf)}
exp_plt(p4,'D')
#ggsave(plot = F3,filename = '../output/F3.png',
# width = 183,height = 125,dpi = 200,units = 'mm')
#!/usr/bin/env Rscript
library(tidyverse)
library(scales)
library(cowplot)
library(grid)
library(gridSVG)
library(grImport2)
source('../../0_data/0_scripts/F4.functions.R')
plts <- trplot(1)
tN <- theme(legend.position = 'none')
pG1 <- ggplotGrob(plts+tN)
pGr1 <- editGrob(pG1, vp=viewport(x=0.5, y=0.91, angle=45,width = .7), name="pG1")
source('../../0_data/0_scripts/F4.genomeWide_box.R')
source('../../0_data/0_scripts/F4.peakArea_box.R')
globGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/caribbean-cairo.svg"))))
leg <- get_legend(plts+theme(legend.direction = 'horizontal'))
bt <-theme(legend.position = 'none',plot.title = element_text(size=8))
bl <- get_legend(pBOX)
sY <- .035;
sX <- -.03;
S16 <- ggdraw(pGr1)+
draw_plot(pBOX+bt,-.012,0,.4,.45)+
draw_plot(boxGenes+bt,.61,0,.4,.45)+
draw_grob(globGrob, 0.01, .88, 0.07, 0.07)+
draw_grob(globGrob, 0.01, .88, 0.07, 0.07)+
draw_plot_label(letters[1:3],
x = c(rep(0.01,2),.7),
y=c(.99,.5,.5),
size = 14)
#plot_grid((ggdraw()+draw_grob(leg,-.45,-.05,1,1)),S16,ncol = 1,rel_heights = c(.1,1)) %>%
ggsave(plot = S16,filename = '~/Desktop/LD.png',width = 16/1.7,height = 9/1.7,dpi = 150)
......@@ -16,5 +16,5 @@ p1 <- ggdraw()+
draw_grob(GA, 0, .84, .25, .2)+
draw_grob(HP, 0.75,.83,.25,.2)
#ggsave(plot = p1,filename = '../output/S01.pdf',width = 183,height = 183,units = 'mm',device = cairo_pdf)
ggsave(plot = p1,filename = '../output/S01.png',width = 183,height = 183,units = 'mm',dpi = 200)
\ No newline at end of file
ggsave(plot = p1,filename = '../output/S01.pdf',width = 183,height = 183,units = 'mm',device = cairo_pdf)
#ggsave(plot = p1,filename = '../output/S01.png',width = 183,height = 183,units = 'mm',dpi = 200)
\ No newline at end of file
......@@ -42,14 +42,14 @@ honGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/hond
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]])
leg <- get_legend(plts[[1]]+theme(legend.key.height = unit(36,'pt')))
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(leg,.24,-.04,1,1)+
draw_grob(pGr3,0,0,1,1)+
draw_grob(pGr4,0,0,1,1)+
draw_grob(pGr5,0,0,1,1)+
......@@ -61,7 +61,7 @@ S09 <- ggdraw(pGr1)+
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_grob(globGrob, 0.01, .78, 0.14, 0.08)+
draw_plot_label(letters[1:3],x = c(rep(0.01,2),.52),
y=c(.99,.57,.57), size = 14)
......
......@@ -32,7 +32,7 @@ cM1 <- hypo_map1.9 %>% group_by(LG) %>%
ungroup() %>%select(Loci,cM,abs_cM)
bp1 <- read.csv('../../1-output/10_recLandscape/MAP1Mb.txt',
sep='\t',head=F,col.names = c("Loci","LG",'bp'))
sep='\t',head=F,col.names = c("Loci","LG",'bp'))
bp1$LG <- factor(as.character(bp1$LG),levels=as.character(seqStart$LG))
bp1 <- bp1 %>%
arrange(LG,bp) %>% group_by(LG) %>%
......@@ -141,7 +141,7 @@ BINWDTH <- 50000
outl_window <- function(LGBIN){
lg = str_sub(LGBIN, 1, 4)
bin = str_sub(LGBIN, 5) %>% as.numeric()
if (lg == "LG09") {
out <- c("Other LGs","Candidate intervals")[between(x = bin,left = 17821000,right = 17930000-BINWDTH)+1]
return(out)
......@@ -169,37 +169,63 @@ an2 <- paste0('r^2:',round(sumMOD$adj.r.squared,3))
an3 <- paste0("italic(F[ST]) == ",modA,modB*10^5," %.% 10^-5"," %.% italic(\u03C1)")
S12c <- ggplot(data_both,aes(x=RHO,y=WEIGHTED_FST))+
geom_point(data=(data_both %>% filter(gr=="Other LGs")),aes(color=gr),size=.7)+
geom_point(data=(data_both %>% filter(gr=="LG08")),aes(color=gr),size=.7)+
geom_point(data=(data_both %>% filter(gr=="Candidate intervals")),aes(color=gr),size=.7)+
geom_point(data=(data_both %>% filter(gr=="Other LGs")),aes(color=gr),size=.5)+
geom_point(data=(data_both %>% filter(gr=="LG08")),aes(color=gr),size=.5)+
geom_point(data=(data_both %>% filter(gr=="Candidate intervals")),aes(color=gr),size=.5)+
geom_smooth(method='lm',
col=modclr,fill=rgb(.9,.9,.9,.001))+
geom_rect(inherit.aes = FALSE,
data=tibble(x1=675,x2=915,y1=.46,y2=.58),
aes(xmin=x1,xmax=x2,ymin=y1,ymax=y2),
col=modclr,fill=NA,lwd=.25)+
annotate(geom = 'text',x=c(795),y=c(.55),label=an3,size=3,parse=TRUE,col=modclr)+
annotate(geom = 'text',x=c(745),y=c(.5),label=an1,size=3)+
annotate(geom = 'text',x=c(865),y=c(.5),label=c(an2),parse=TRUE,size=3)+
annotate(geom = 'text',x=c(795)-150,y=c(.58),label=an3,size=2.5,parse=TRUE,col=modclr)+
annotate(geom = 'text',x=c(745)-215,y=c(.55),label=an1,size=2.5)+
annotate(geom = 'text',x=c(865)-105,y=c(.55),label=c(an2),parse=TRUE,size=2.5)+
scale_x_continuous(name=expression(Population~recombination~rate~(italic("\u03C1"))),
expand = c(.01,.01))+
scale_y_continuous(name=expression(Genetic~differentiation~(weighted~mean~italic(F[ST]))),
expand = c(.005,.005))+
scale_color_manual(name = 'Group',values = c('#d40606',rgb(.5,.5,.5,.8),'black'))+
theme_gw+theme(plot.background = element_blank(),
panel.background = element_blank(),
panel.grid=element_blank(),
panel.border = element_blank(),
plot.margin = unit(c(0,15,0,25),'pt'),
axis.line.x = element_line(color = annoclr),
axis.text.x = element_text(size=6),
axis.title.x = element_text(size=8),
legend.position = c(.15,.82),
legend.background = element_rect(color = secclr,fill=NA),
legend.title = element_blank())
panel.background = element_blank(),
panel.grid=element_blank(),
panel.border = element_blank(),
plot.margin = unit(c(0,15,0,25),'pt'),
axis.line.x = element_line(color = annoclr),
axis.text.x = element_text(size=6),
axis.title.x = element_text(size=8),
legend.position = 'bottom',
legend.title = element_blank())
decay_path <- '../../2_output/08_popGen/07_LD/decay/'
files_d <- dir(path = decay_path,pattern = 'ld.gz')
read_r2 <- function(file){read_delim(file = str_c(decay_path,file),delim = '\t')}
data_d <- files_d %>% purrr::map(.,read_r2) %>% bind_rows()
S12d <- ggplot(data=data_d,aes(x=(POS2-POS1)/1000,y=`R^2`))+
geom_hex(bins = 100,aes(fill=log10(..count..))) +
geom_smooth(col='red',se=FALSE,size=.3,
method = 'gam', formula = y ~ s(x, bs = "cs"))+
scale_fill_gradient(name = expression(Count~(log[10])),
low = rgb(.9,.9,.9), high = "black")+
guides(fill=guide_colorbar(direction = 'horizontal',title.position = 'top'))+
scale_x_continuous(name = "Distance (kb)",expand = c(0,0))+
scale_y_continuous(name = expression(Linkage~(r^2)),
expand = c(0.002,0.002))+
theme_gw+
theme(plot.background = element_blank(),
panel.background = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
plot.margin = unit(c(0,15,0,25),'pt'),
axis.line = element_line(color = annoclr),
axis.text = element_text(size=6),
axis.title.x = element_text(size=8),
axis.title.y = element_text(size=8),
legend.title = element_text(size=8),
legend.text = element_text(size=6),
legend.key.height = unit(4,'pt'),
legend.position = 'bottom')
S12ab <- plot_grid(S12a,S12b,ncol = 1,align = 'v',labels = c('a','b'))
S12 <- plot_grid(S12ab,NULL,S12c,ncol = 1,rel_heights = c(1,.05,.8),labels = c('','','c'))
S12cd <- plot_grid(S12c,S12d,nrow = 1,labels = c('c','d'),align = 'h')
S12 <- plot_grid(S12ab,NULL,S12cd,ncol = 1,rel_heights = c(1,.05,.8))
ggsave(plot = S12,filename = '../output/S12.pdf',width = 183,height = 220,units = 'mm',device = cairo_pdf)
#ggsave(plot = S12,filename = '../output/S12.png',width = 183,height = 220,units = 'mm',dpi = 200)
#ggsave(plot = S12,filename = '../output/S12.png',width = 183,height = 220,units = 'mm',dpi = 200)
......@@ -27,12 +27,13 @@ puGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/uni-p
nuGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/nig-uni-cairo.svg"))))
leg <- get_legend(plts[[1]]+theme(legend.text = element_text(size = 5),
legend.title = element_text(size = 5),
legend.key.size = unit(7,'pt')))
legend.title = element_text(size = 5),
legend.key.size = unit(7,'pt'),
legend.key.height = unit(22,'pt')))
S13 <- ggdraw(pGr8)+
draw_grob(pGr9,0,0,1,1)+
draw_grob(leg,-.85,.04,1,1)+
draw_grob(leg,-.85,.03,1,1)+
draw_grob(pGr10,0,0,1,1)+
draw_grob(npGrob, 0.7, 0.7, 0.28, 0.1)+
draw_grob(nuGrob, 0.7, 0.37, 0.28, 0.1)+
......
......@@ -27,11 +27,12 @@ nuGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/nig-u
leg <- get_legend(plts[[1]]+theme(legend.text = element_text(size = 5),
legend.title = element_text(size = 5),
legend.key.size = unit(7,'pt')))
legend.key.size = unit(7,'pt'),
legend.key.height = unit(22,'pt')))
S14 <- ggdraw(pGr8)+
draw_grob(pGr9,0,0,1,1)+
draw_grob(leg,-.65,-.02,1,1)+
draw_grob(leg,-.65,0.02,1,1)+
draw_grob(pGr10,0,0,1,1)+
draw_grob(npGrob, 0.7, 0.65, 0.28, 0.1)+
draw_grob(nuGrob, 0.7, 0.34, 0.28, 0.1)+
......
......@@ -12,7 +12,6 @@ plts <- trplot(1)
tN <- theme(legend.position = 'none')
pG1 <- ggplotGrob(plts+tN)
pGr1 <- editGrob(pG1, vp=viewport(x=0.5, y=0.91, angle=45,width = .7), name="pG1")
source('../../0_data/0_scripts/S16.genomeWide_box.R')
......@@ -24,7 +23,6 @@ leg <- get_legend(plts+theme(legend.direction = 'horizontal'))
bt <-theme(legend.position = 'none',plot.title = element_text(size=8))
bl <- get_legend(pBOX)
sY <- .035;
sX <- -.03;
S16 <- ggdraw(pGr1)+
......
No preview for this file type
3_figures/output/S16.png

152 KB | W: | H:

3_figures/output/S16.png

150 KB | W: | H:

3_figures/output/S16.png
3_figures/output/S16.png
3_figures/output/S16.png
3_figures/output/S16.png
  • 2-up
  • Swipe
  • Onion skin
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment