Commit 83590835 authored by Kosmas Hench's avatar Kosmas Hench

include spotlight LD

include fst.local_allSpec_genome_Wide
parent 0ccac08b
LG09 17850000 17900000 SOX10_1
LG12 20150000 20350000 Casz1_2
LG12 22200000 22275000 hoxc13a
LG17 22525000 22600000 LWS
#PBS -l elapstim_req=12:00:00
#PBS -l memsz_job=70gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N rawLDXXnameXX
#PBS -q clmedium
#PBS -o 2.2.8.1.raw.LD.XXnameXX.stdout
#PBS -e 2.2.8.1.raw.LD.XXnameXX.stderr
cd $WORK/2_output/08_popGen/07_LD
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_XXpop1XX.pop \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_XXpop2XX.pop \
--bed $WORK/0_data/0_resources/genes.bed \
--interchrom-hap-r2 \
--out spotlight.XXnameXX
gzip XXnameXX.interchrom.hap.ld
#PBS -l elapstim_req=12:00:00
#PBS -l memsz_job=70gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N rawLDXXnameXX
#PBS -q clmedium
#PBS -o 2.2.8.1.raw.LD.XXnameXX.stdout
#PBS -e 2.2.8.1.raw.LD.XXnameXX.stderr
cd $WORK/2_output/08_popGen/07_LD
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_XXnameXX.pop \
--bed $WORK/0_data/0_resources/genes.bed \
--interchrom-hap-r2 \
--out spotlight.XXnameXX
gzip XXnameXX.interchrom.hap.ld
#!/bin/bash
# starter script to rewrite and submit the 2.2.8.1.rawLD_pair_temp.sh script
# once per species-pair
cd $WORK/2_reSeq-scripts/2-2_popGen/06_LD_spotlight
for k in nig-pue nig-uni pue-uni;do
POP1=$(echo $k | sed "s/-/\t/" | cut -f 1);
POP2=$(echo $k | sed "s/-/\t/" | cut -f 2);
echo -e "--- run $k ---\npop1: $POP1\npop2: $POP2"
mkdir -p $k
sed "s/XXnameXX/$k/g; s/XXpop1XX/$POP1/g; s/XXpop2XX/$POP2/g;" 00_templates/2.2.8.1.rawLD_pair_temp.sh > ./$k/2.2.8.1.rawLD.$k.sh
cd ./$k
qsub 2.2.8.1.rawLD.$k.sh
cd ..
done
#!/bin/bash
# starter script to rewrite and submit the 2.2.8.1.rawLD_single_temp.sh script
# once per species and once per location
cd $WORK/2_reSeq-scripts/2-2_popGen/06_LD_spotlight
for k in bel hon boc nig pue uni;do
mkdir -p $k
sed "s/XXnameXX/$k/g" 00_templates/2.2.8.1.rawLD_single_temp.sh > ./$k/2.2.8.1.rawLD.$k.sh
cd ./$k
qsub 2.2.8.1.rawLD.$k.sh
cd ..
done
#PBS -l elapstim_req=12:00:00
#PBS -l memsz_job=70gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N rawLDglobal
#PBS -q clmedium
#PBS -o 2.2.8.1.raw.LD.global.stdout
#PBS -e 2.2.8.1.raw.LD.global.stderr
cd $WORK/2_output/08_popGen/07_LD
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--bed $WORK/0_data/0_resources/genes.bed \
--interchrom-hap-r2 \
--out spotlight.global
gzip global.interchrom.hap.ld
#PBS -l elapstim_req=3:00:00
#PBS -l memsz_job=70gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N fstLAsGw
#PBS -q clmedium
#PBS -o local_global_global.stdout
#PBS -e local_global_global.stderr
cd $WORK/2_output/08_popGen/05_fst
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_nigbel.pop \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_puebel.pop \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_unibel.pop \
--out local_allSpec_genome_wide.bel
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_nighon.pop \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_puehon.pop \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_unihon.pop \
--out local_allSpec_genome_wide.hon
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_nigboc.pop \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_pueboc.pop \
--weir-fst-pop $WORK/0_data/0_resources/vcfpops/vcftools_uniboc.pop \
--out local_allSpec_genome_wide.boc
\ No newline at end of file
......@@ -18,3 +18,4 @@ for k in {01..24};
head -n 1 HP.annotation.named.gff > subsets/HP.annotation.named.LG$k.gff;
grep "LG"$k HP.annotation.named.gff >> subsets/HP.annotation.named.LG$k.gff;
done
\ No newline at end of file
#!/usr/bin/env Rscriptlibrary(ggplot2)
library(dplyr)
library(scales)
library(reshape2)
library(cowplot)
library(grid)
library(gridSVG)
library(grImport2)
#library(grConvert)
#install.packages("grConvert")
#data <- read.csv('sox-hox-10000-bins.txt',sep='\t')
sel <- 1
trplot <- function(sel){
nm <- c('no_maf','_maf01','_maf01.boc','_maf01.bel','_maf01.hon','_maf01.pue','_maf01.nig','_maf01.uni')
#nm <- c('no_maf','_maf01_incl_perfect','_maf01.boc','_maf01.bel','_maf01.hon')
if(sel == 0){
check = 'no maf';
data <- read.csv('sox-hox-10000-bins.txt',sep='\t')
datadens <- read.csv(gzfile('sox-hox-ops.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 1){
check = 'global';
data <- read.csv('altData/sox-hox-10000-bins.txt',sep='\t')
#datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 2){
check = 'bocas';
data <- read.csv('altData/sox-hox_boc-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.boc.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 3){
check = 'belize';
data <- read.csv('altData/sox-hox_bel-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.bel.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 4){
check = 'honduras';
data <- read.csv('altData/sox-hox_hon-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 5){
check = 'puella';
data <- read.csv('altData/sox-hox_pue-10000-bins.txt',sep='\t')
#datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 6){
check = 'nigricans';
data <- read.csv('altData/sox-hox_nig-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 7){
check = 'unicolor';
data <- read.csv('altData/sox-hox_uni-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 8){
check = 'nig-pue';
data <- read.csv('altData/sox-hox_nig-pue-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 9){
check = 'nig-uni';
data <- read.csv('altData/sox-hox_nig-uni-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 10){
check = 'pue-uni';
data <- read.csv('altData/sox-hox_pue-uni-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}
print(check)
plim <- c(-20,250)
s1=17620000;s2=19910000;s3=21960000;s4=22320000;e2=20660000;e3=22460000;stp=20000;l1=500000;l2=750000;l3=500000;l4=500000
scling <- c(-s1,-s2+(l1+stp),-s3+(l1+l2+2*stp),-s4+(l1+l2+l3+(stp*3)))
dt <- data %>% mutate(miX = floor(Mx/10000),miY=floor(My/10000))
genes <- data.frame(start=c(17871610,20186151,22225149,22553187,22556763,22561894,22573388),
end=c(17873915,20347811,22228342,22555052,22559742,22566321,22575503),
sclr=c(1,2,3,4,4,4,4),
LG=c("LG09","LG12-1","LG12-2","LG17","LG17","LG17","LG17"),
name=c("sox10",'casz1',"hoxc13a","sws2a\u03B1","sws2a\u03B2","sws2b","lws"))
genes <- genes %>% mutate(Nx1 = (start+scling[sclr])/10000,Nx2 = (end+scling[sclr])/10000,
labPOS = (Nx1+Nx2)/2) #%>%
genes$labPOS[genes$name %in% c("sws2a\u03B1","sws2a\u03B2","sws2b","lws")] <- genes$labPOS[genes$name %in% c("sws2a\u03B1","sws2a\u03B2","sws2b","lws")]+c(-16,-1,12,20)
#melt(id.vars=c("start","end","sclr","LG","name")) %>% arrange(LG,name)
GZrS <- 40000;GZrS2 <- 20000; #gene zoom offset
BS <- c(-5,5,-4,3,-8,8,-20,-10,-6,5,7,16,17,22)*10000 # Backshifter for gene zoom
GZdf <- data.frame(x=c(genes$start[1],genes$end[1],genes$end[1]+GZrS+BS[1],genes$start[1]+GZrS+BS[2],genes$start[1],
genes$start[2],genes$end[2],genes$end[2]+GZrS+BS[3],genes$start[2]+GZrS+BS[4],genes$start[2],
genes$start[3],genes$end[3],genes$end[3]+GZrS+BS[5],genes$start[3]+GZrS+BS[6],genes$start[3],
genes$start[4],genes$end[4],genes$end[4]+GZrS+BS[7],genes$start[4]+GZrS+BS[8],genes$start[4],
genes$start[5],genes$end[5],genes$end[5]+GZrS+BS[9],genes$start[5]+GZrS+BS[10],genes$start[5],
genes$start[6],genes$end[6],genes$end[6]+GZrS+BS[11],genes$start[6]+GZrS+BS[12],genes$start[6],
genes$start[7],genes$end[7],genes$end[7]+GZrS+BS[13],genes$start[7]+GZrS+BS[14],genes$start[7])+GZrS2,
y=c(genes$start[1],genes$end[1],genes$end[1]-GZrS+BS[1],genes$start[1]-GZrS+BS[2],genes$start[1],
genes$start[2],genes$end[2],genes$end[2]-GZrS+BS[3],genes$start[2]-GZrS+BS[4],genes$start[2],
genes$start[3],genes$end[3],genes$end[3]-GZrS+BS[5],genes$start[3]-GZrS+BS[6],genes$start[3],
genes$start[4],genes$end[4],genes$end[4]-GZrS+BS[7],genes$start[4]-GZrS+BS[8],genes$start[4],
genes$start[5],genes$end[5],genes$end[5]-GZrS+BS[9],genes$start[5]-GZrS+BS[10],genes$start[5],
genes$start[6],genes$end[6],genes$end[6]-GZrS+BS[11],genes$start[6]-GZrS+BS[12],genes$start[6],
genes$start[7],genes$end[7],genes$end[7]-GZrS+BS[13],genes$start[7]-GZrS+BS[14],genes$start[7])-GZrS2,
grp=rep(letters[1:7],each=5)) %>%
mutate(sclr=rep(c(1,2,3,4,4,4,4),each=5),
Nx1 = (x+scling[sclr])/10000,
Nx2 = (y+scling[sclr])/10000)
clr = c(viridis::inferno(5)[c(1,1:5)])
#clr = c('black',colorRampPalette(c('#542788',"#fe9929"))(5))
#clr = c('black',rev(RColorBrewer:: brewer.pal(5, "PuOr"))[c(1:5)])
#clr = c('black','black',rev(rev(viridis::viridis(4))))
Gcol <- '#3bb33b'#viridis_pal(option = 'D')(5)[1]
Zcol = rgb(.94,.94,.94)
LGoffset <- 15
GLABoffset <- 8
if(sel %in% c(0,1,8)){
rS <- 100000 # width of grey annotation band
zmRange <- data.frame(x=c(s1,s1+l1,s1+l1+rS,s1+rS,s1,
s2,s2+l2,s2+l2+rS,s2+rS,s2,
s3,s3+l3,s3+l3+rS,s3+rS,s3,
s4,s4+l4,s4+l4+rS,s4+rS,s4),
y=c(s1,s1+l1,s1+l1-rS,s1-rS,s1,
s2,s2+l2,s2+l2-rS,s2-rS,s2,
s3,s3+l3,s3+l3-rS,s3-rS,s3,
s4,s4+l4,s4+l4-rS,s4-rS,s4),
grp=rep(letters[1:4],each=5)) %>%
mutate(sclr=rep(1:4,each=5),
Nx1 = (x+scling[sclr])/10000-1,
Nx2 = (y+scling[sclr])/10000)
rS2 <- .75*rS
zmRange2 <- data.frame(x=c(s1,s1+l1,s1+l1+rS2,s1+rS2,s1,
s2,s2+l2,s2+l2+rS2,s2+rS2,s2,
s3,s3+l3,s3+l3+rS2,s3+rS2,s3,
s4,s4+l4,s4+l4+rS2,s4+rS2,s4),
y=c(s1,s1+l1,s1+l1-rS2,s1-rS2,s1,
s2,s2+l2,s2+l2-rS2,s2-rS2,s2,
s3,s3+l3,s3+l3-rS2,s3-rS2,s3,
s4,s4+l4,s4+l4-rS2,s4-rS2,s4),
grp=rep(letters[1:4],each=5)) %>%
mutate(sclr=rep(1:4,each=5),
Nx1 = (x+scling[sclr])/10000-1,
Nx2 = (y+scling[sclr])/10000)
rS3 <- .07*rS
zmRange3 <- data.frame(x=c(s1,s1+l1,s1+l1+rS3,s1+rS3,s1,
s2,s2+l2,s2+l2+rS3,s2+rS3,s2,
s3,s3+l3,s3+l3+rS3,s3+rS3,s3,
s4,s4+l4,s4+l4+rS3,s4+rS3,s4),
y=c(s1,s1+l1,s1+l1-rS3,s1-rS3,s1,
s2,s2+l2,s2+l2-rS3,s2-rS3,s2,
s3,s3+l3,s3+l3-rS3,s3-rS3,s3,
s4,s4+l4,s4+l4-rS3,s4-rS3,s4),
grp=rep(letters[1:4],each=5)) %>%
mutate(sclr=rep(1:4,each=5),
Nx1 = (x+scling[sclr])/10000-1,
Nx2 = (y+scling[sclr])/10000)
zmLab <- data.frame(x=c(s1+.5*l1,s2+.5*l2,s3+.5*l3,s4+.5*l4),
label=c('LG09 (A)','LG12 (B)','LG12 (C)','LG17 (D)')) %>%
mutate(sclr=1:4, Nx= (x+scling[sclr])/10000)
zmEND <- data.frame(x=c(s1,s1+l1,
s2,s2+l2,
s3,s3+l3,
s4,s4+l4)) %>%
mutate(sclr=rep(1:4,each=2),
Nx = ((x+scling[sclr])/10000)+rep(c(2.5,-4),4),
label=round((x/1000000),2))
DG <- rgb(.4,.4,.4)
if(sel %in% c(0,1)){
p1 <- ggplot(dt %>% filter(!is.na(Mval)),aes(fill=Mval))+#geom_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),alpha=1)+
coord_equal()+
geom_polygon(inherit.aes = F,data=zmRange,
aes(x=Nx1+1,y=Nx2-1,group=grp),fill='lightgray')+
geom_polygon(inherit.aes = F,data=zmRange2,
aes(x=Nx1+11,y=Nx2-11,group=grp),fill=DG)+
geom_polygon(inherit.aes = F,data=zmRange3,
aes(x=Nx1+1.1,y=Nx2-1.1,group=grp),fill=DG)+
geom_polygon(inherit.aes = F,data=GZdf,
aes(x=Nx1,y=Nx2,group=grp),fill=Zcol)+
geom_segment(inherit.aes = F,data=genes,
aes(x=Nx1+1,y=Nx1-1,xend=Nx2+1,yend=Nx2-1),
col=Gcol,size=1.5)+
geom_tile(aes(x=miX,y=miY))+
geom_text(inherit.aes = F,data=zmEND,
aes(x=Nx+LGoffset,y=Nx-LGoffset,label=paste(label,'\n(Mb)')),
angle=45,size=1.8)+
geom_text(inherit.aes = F,data=genes,
aes(x=labPOS+GLABoffset,y=labPOS-GLABoffset,label=name),
angle=-45,fontface='italic',size=2)+
geom_text(inherit.aes = F,data=zmLab,
aes(x=Nx+LGoffset-.8,y=Nx-LGoffset+.8,label=label),
angle=-45,fontface='bold',size=3.5,col='white')+
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
# stat_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),aes(alpha=..level..), geom="polygon")+
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_brewer(name="",type="qual",palette = 6)+
scale_color_manual(values=viridis_pal(option='A')(6))+theme_void()+
theme(legend.position = c(.7,.75),legend.direction = 'vertical')
} else if(sel %in% c(8)){
p1 <- ggplot(dt %>% filter(!is.na(Mval)),aes(fill=Mval))+#geom_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),alpha=1)+
coord_equal()+
geom_polygon(inherit.aes = F,data=zmRange,
aes(x=Nx1+1,y=Nx2-1,group=grp),fill='lightgray')+
geom_polygon(inherit.aes = F,data=zmRange2,
aes(x=Nx1+11,y=Nx2-11,group=grp),fill=DG)+
geom_polygon(inherit.aes = F,data=zmRange3,
aes(x=Nx1+1.1,y=Nx2-1.1,group=grp),fill=DG)+
geom_polygon(inherit.aes = F,data=GZdf,
aes(x=Nx1,y=Nx2,group=grp),fill=Zcol)+
geom_segment(inherit.aes = F,data=genes,
aes(x=Nx1+1,y=Nx1-1,xend=Nx2+1,yend=Nx2-1),
col=Gcol,size=1.5)+
geom_tile(aes(x=miX,y=miY))+
geom_text(inherit.aes = F,data=zmEND,
aes(x=Nx+LGoffset,y=Nx-LGoffset,label=paste(label,'\n(Mb)')),
angle=45,size=.7)+
geom_text(inherit.aes = F,data=genes,
aes(x=labPOS+GLABoffset,y=labPOS-GLABoffset,label=name),
angle=-45,fontface='italic',size=1)+
geom_text(inherit.aes = F,data=zmLab,
aes(x=Nx+LGoffset-.8,y=Nx-LGoffset+.8,label=label),
angle=-45,fontface='bold',size=1.75,col='white')+
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
# stat_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),aes(alpha=..level..), geom="polygon")+
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_brewer(name="",type="qual",palette = 6)+
scale_color_manual(values=viridis_pal(option='A')(6))+theme_void()+
theme(legend.position = c(.7,.75),legend.direction = 'vertical')
}
#dfsub <- dt %>% filter(!is.na(Mval) )
#p2 <- ggplot(dfsub)+#geom_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),alpha=1)+ +
# geom_density(aes(x=Mval))+ggtitle(paste(check,round(mean(dfsub$Mval),5),round(median(dfsub$Mval),5)))
} else {
p1 <- ggplot(dt %>% filter(!is.na(Mval)),aes(fill=Mval))+#geom_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),alpha=1)+
coord_equal()+
geom_segment(inherit.aes = F,data=genes,
aes(x=Nx1+1.5,y=Nx1-1.5,xend=Nx2+1.5,yend=Nx2-1.5),
col=Gcol,size=1)+
geom_tile(aes(x=miX,y=miY))+
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0),
trans = 'reverse')+
# stat_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),aes(alpha=..level..), geom="polygon")+
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_brewer(name="",type="qual",palette = 6)+
scale_color_manual(values=viridis_pal(option='A')(6))+theme_void()+
theme(legend.position = c(.7,.75),legend.direction = 'vertical')
# dfsub <- dt %>% filter(!is.na(Mval) )
# p2 <- ggplot(dfsub)+#geom_density2d(inherit.aes = F,data = d2t,aes(x=x,y=y),alpha=1)+ +
# geom_density(aes(x=Mval))+ggtitle(paste(check,round(mean(dfsub$Mval),5),round(median(dfsub$Mval),5)))
}
return(p1)
#print(p1)
#ggsave(plot = ,filename = paste('ld-triangle',nm[sel+1],'.pdf',sep=''),width = 15,height = 15)
}
plts <- list()
for(k in 0:10){
plts[[k+1]] <- trplot(k)
}
# plot_grid(plts[[1]][2][[1]]+scale_x_continuous(limits = c(0,.1)),
# plts[[2]][2][[1]]+scale_x_continuous(limits = c(0,.1)),
# plts[[3]][2][[1]]+scale_x_continuous(limits = c(0,.1)),
# plts[[4]][2][[1]]+scale_x_continuous(limits = c(0,.1)),
# plts[[5]][2][[1]]+scale_x_continuous(limits = c(0,.1)),
# plts[[6]][2][[1]]+scale_x_continuous(limits = c(0,.1)),
# plts[[7]][2][[1]]+scale_x_continuous(limits = c(0,.1)),
# plts[[8]][2][[1]]+scale_x_continuous(limits = c(0,.1)),ncol = 1,align = 'v')
#
tN <- theme(legend.position = 'none')
pG1 <- ggplotGrob(plts[[2]]+tN)
pG2 <- ggplotGrob(plts[[3]]+tN)
pG3 <- ggplotGrob(plts[[4]]+tN)
pG4 <- ggplotGrob(plts[[5]]+tN)
pG5 <- ggplotGrob(plts[[7]]+tN)
pG6 <- ggplotGrob(plts[[6]]+tN)
pG7 <- ggplotGrob(plts[[8]]+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("img/nigricans-cairo.svg"))))
pueGrob <- gTree(children=gList(pictureGrob(readPicture("img/puella-cairo.svg"))))
uniGrob <- gTree(children=gList(pictureGrob(readPicture("img/unicolor-cairo.svg"))))
belGrob <- gTree(children=gList(pictureGrob(readPicture("img/belize-cairo.svg"))))
honGrob <- gTree(children=gList(pictureGrob(readPicture("img/honduras-cairo.svg"))))
panGrob <- gTree(children=gList(pictureGrob(readPicture("img/panama-cairo.svg"))))
globGrob <- gTree(children=gList(pictureGrob(readPicture("img/caribbean-cairo2.svg"))))
library(tidyr)
BW <- read.csv('subs/glob_between.interchrom.hap.ld',sep='\t')%>% select(R.2) %>% mutate(type='Global',run='Global')
for(j in 1:6){
flS <- c("boc","bel","hon","pue","nig","uni")
flL <- c("Panana","Belize","Honduras","H. puella","H. nigricans","H. unicolor")
flT <- c(rep('Geo',3),rep('Spec',3))
k <- flS[j]
q <- flL[j]
u <- flT[j]
print(j)
# WI <- WI %>%
# rbind(.,(read.csv(paste('glob_within.',k,'.hap.ld',sep=''),sep='\t') %>% select(R.2) %>% mutate(type='wi',run=k)))
BW <- BW %>%
rbind(.,(read.csv(paste('subs/glob_between.',k,'.interchrom.hap.ld',sep=''),sep='\t') %>% select(R.2) %>% mutate(type=u,run=q)))
}
#dt <- rbind(WI,BW)
BW$run <- factor(BW$run,levels=c('Global',"Panana","Belize","Honduras",
"H. nigricans","H. puella","H. unicolor"))
dt2 <- BW %>%
group_by(run) %>%
summarise(mn=mean(R.2,na.rm = T),md=median(R.2,na.rm = T)) %>%
gather(key = 'type',value = 'val',2:3)
BC <-rgb(.7,.7,.7)#viridis_pal(option='A')(7)[6]
clr <-colorRampPalette(colors = c('white',BC,'black'))(8)
pBOX <- ggplot(BW,aes(x=run,y=R.2))+#geom_violin(fill='#3bb33b')+
geom_boxplot(fill=BC,width=.7,outlier.size = .1)+
# coord_cartesian(ylim=c(0,.031))+
coord_fixed(ylim=c(0,.031),ratio = 250)+
geom_point(inherit.aes = F, data=dt2,aes(x=run,y=val,fill=type),shape=23,size=1)+
scale_x_discrete(labels = expression(Global,Panama,Belize,Honduas,italic("H. nigricans"),italic("H. puella"),italic("H. unicolor")))+
scale_y_continuous('r²')+
scale_fill_manual('',values = clr[c(1,6)],labels=c('median','mean'))+
# scale_color_manual('',values = clr[c(2,4)],labels=c('median','mean'))+
# scale_shape_manual(values=c(21,23))+
theme(legend.position = c(-.4,1.15),legend.direction = 'horizontal',
text = element_text(size=10),
axis.title.x = element_blank(),
axis.text.y = element_text(size=7),
axis.text.x = element_text(size=7,angle=45,hjust = 1),
panel.background = element_blank(),
plot.background = element_blank());pBOX
leg <- get_legend(plts[[2]])
sY <- .035;
sX <- -.03;
ggdraw(pGr1)+
draw_plot(pBOX,-.1,.55,.45,.21)+
draw_plot(pGRAD,.307,-.018,.16,.53)+
draw_grob(pGr2,0,0,1,1)+
draw_grob(leg,.15,-.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:4],x = c(rep(0.01,3),.52),
y=c(.99,.81,.57,.57), size = 14)
# draw_plot_label(letters[1:7],
# x = c(rep(0.01,4),rep(.52,3)),
# y=c(.99,rep(c(.57,.39,.2),2)),
# size = 14)+
#getwd()
#ggsave('ms_fig_4.pdf',width = 183,height = 210,units = 'mm',device = cairo_pdf)
ggsave('LDtriangle.pdf',width = 183,height = 230,units = 'mm',device = cairo_pdf)
ggsave('LDtriangle.png',width = 183,height = 230,units = 'mm',dpi = 150)
pG8 <- ggplotGrob(plts[[9]]+tN)
pG9 <- ggplotGrob(plts[[10]]+tN)
pG10 <- ggplotGrob(plts[[11]]+tN)
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("img/pue-nig-cairo.svg"))))
puGrob <- gTree(children=gList(pictureGrob(readPicture("img/uni-pue-cairo.svg"))))
nuGrob <- gTree(children=gList(pictureGrob(readPicture("img/nig-uni-cairo.svg"))))
leg2 <- get_legend(plts[[2]]+theme(legend.text = element_text(size = 5),
legend.title = element_text(size = 5),
legend.key.size = unit(7,'pt')))
ggdraw(pGr8)+
draw_grob(pGr9,0,0,1,1)+
draw_grob(leg2,-.65,.04,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)+
draw_grob(puGrob, 0.7, .04, 0.28, 0.1)
ggsave('LDtriangle_pw.pdf',width = 91.5,height = 145,units = 'mm',device = cairo_pdf)
ggsave('LDtriangle_pw.png',width = 91.5,height = 145,units = 'mm',dpi = 150)
\ No newline at end of file
setwd('~/reSEQ/LDspoiler/')
library(ggplot2)
library(dplyr)
library(scales)
library(reshape2)
library(cowplot)
library(grid)
library(gridSVG)
library(grImport2)
#library(grConvert)
#install.packages("grConvert")
#data <- read.csv('sox-hox-10000-bins.txt',sep='\t')
sel <- 1
trplot <- function(sel){
nm <- c('mac02','_maf01.boc','_maf01.bel','_maf01.hon','_maf01.pue','_maf01.nig','_maf01.uni')
#nm <- c('no_maf','_maf01_incl_perfect','_maf01.boc','_maf01.bel','_maf01.hon')
if(sel == 1){
check = 'global';
data <- read.csv('extData/extLD-10000-bins.txt',sep='\t')
#datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 2){
check = 'bocas';
data <- read.csv('extData/extLD.boc-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.boc.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 3){
check = 'belize';
data <- read.csv('extData/extLD.bel-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.bel.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 4){
check = 'honduras';
data <- read.csv('extData/extLD.hon-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 5){
check = 'puella';
data <- read.csv('extData/extLD.pue-10000-bins.txt',sep='\t')
#datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 6){
check = 'nigricans';
data <- read.csv('extData/extLD.nig-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 7){
check = 'unicolor';
data <- read.csv('extData/extLD.uni-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 8){
check = 'nigricans-puella';
data <- read.csv('extData/extLD.nig-pue-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.boc.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 9){
check = 'nigricans-unicolor';
data <- read.csv('extData/extLD.nig-uni-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.bel.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}else if(sel == 10){
check = 'puella-unicolor';
data <- read.csv('extData/extLD.pue-uni-10000-bins.txt',sep='\t')
# datadens <- read.csv(gzfile('altData/sox-hox-ops.maf01.hon.perfect_linked.txt.gz'),sep ='\t') %>% filter(!is.na(R.2))
}
print(check)
plim <- c(-20,250)
s1=5760000;s2=1884000;s3=12860000;s4=17745000;s5=20149000;s6=22085000;s7=22445000;s8=13730000
e2=2134000;e3=13110000;e5=20399000;e6=22335000
stp=20000;lALL=250000#;l2=250000;l3=250000;l4=250000;l5=250000;l6=250000;l7=250000;l8=250000
scling <- c(-s1,
-s2+(lALL+stp),
-s3+((lALL+stp)*2),
-s4+((lALL+stp)*3),
-s5+((lALL+stp)*4),
-s6+((lALL+stp)*5),
-s7+((lALL+stp)*6),
-s8+((lALL+stp)*7))
dt <- data %>% mutate(miX = floor(Mx/10000),miY=floor(My/10000))
genes <- data.frame(start=c(5884279,2009113,12992396,17871610,20186151,22228342,22553187,22556763,22561894,22573388,13862931),
end=c(5878614,2010303,12976894,17873915,20347811,22225149,22555052,22559742,22566321,22575503,13884003),
sclr=c(1,2,3,4,5,6,7,7,7,7,8),
LG=c("LG04","LG08","LG08","LG09","LG12","LG12","LG17","LG17","LG17","LG17","LG20"),
name=c("ednrb","foxd3","rorb","sox10","casz1","hoxc13a","sws2a\u03B2","sws2a\u03B1","sws2b","lws","invs"))
genes <- genes %>% mutate(Nx1 = (start+scling[sclr])/10000,Nx2 = (end+scling[sclr])/10000,
labPOS = (Nx1+Nx2)/2) #%>%
genes$labPOS[genes$name %in% c("sws2a\u03B1","sws2a\u03B2","sws2b","lws")] <- genes$labPOS[genes$name %in% c("sws2a\u03B1","sws2a\u03B2","sws2b","lws")]+c(-7,