Commit 06738b06 authored by Kosmas Hench's avatar Kosmas Hench

relabel -> new format

parent 79ce4b8c
get_anno_df_single_line <- function(searchLG,gfffile,
xrange,
genes_of_interest,
genes_of_sec_interest,
anno_rown=3){
gff_filter <- list(seqid=searchLG)
data <- as.data.frame(readGFF(gfffile,filter=gff_filter)) %>% mutate(Parent=as.character(Parent))#%>%
#rowwise() %>% mutate(Parent=ifelse(length(Parent)==0,ID,Parent))
mRNAs <- data %>% filter(type=='mRNA',end>xrange[1],start<xrange[2]) %>%
ungroup() %>% mutate(yl=row_number()%%anno_rown+2) %>% rowwise()%>%
mutate(checkStart =ifelse(start<xrange[1],-Inf,start),
checkEnd =ifelse(end>xrange[2],Inf,end),
ps=ifelse(strand=='-',checkEnd,checkStart),
pe=ifelse(strand=='-',checkStart,checkEnd),
labelx=mean(c(sort(c(xrange[1],ps))[2],
sort(c(xrange[2],pe))[1])),
window='bold(Gene)',
clr=ifelse(Parentgenename %in% genes_of_interest,"y",
ifelse(Parentgenename %in% genes_of_sec_interest,"z","x"))) %>%
select(-Parent);
names(mRNAs)[names(mRNAs)=='ID'] <- 'Parent'
exons <- data %>% filter(type=='exon',end>xrange[1],start<xrange[2]) %>%
merge(.,mRNAs %>% select(Parent,yl,clr),by='Parent',all.x=T) %>%
mutate(ps=ifelse(strand=='-',end,start),
pe=ifelse(strand=='-',start,end),
window='bold(Gene)')
return(list(mRNAs,exons))}
This diff is collapsed.
create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchsnp,muskID){ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchsnp,muskID){
source('../../0_data/0_scripts/E1.functions.R'); source('../../0_data/0_scripts/F3.functions.R');
source('../../0_data/0_scripts/E1.getFSTs.R') source('../../0_data/0_scripts/F3.getFSTs.R')
source('../../0_data/0_scripts/E1.getGxP.R') source('../../0_data/0_scripts/F3.getGxP.R')
source('../../0_data/0_scripts/E1.getDXY.R') source('../../0_data/0_scripts/F3.getDXY.R')
highclr <- '#3bb33b' highclr <- '#3bb33b'
theme_set(theme_minimal(base_size = 6)) theme_set(theme_minimal(base_size = 6))
...@@ -87,24 +87,24 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs ...@@ -87,24 +87,24 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
panel.border = element_blank() panel.border = element_blank()
); );
if(muskID=='A'){ if(muskID=='a'){
p11$layers[[4]]$data$label <- c("italic(sox10)","italic(rnaseh2a)") p11$layers[[4]]$data$label <- c("italic(sox10)","italic(rnaseh2a)")
} else if(muskID=='B'){ } else if(muskID=='b'){
p11$layers[[4]]$data$label <- c("italic(casz1.3)","italic(casz1.2)","italic(casz1.1)") p11$layers[[4]]$data$label <- c("italic(casz1.3)","italic(casz1.2)","italic(casz1.1)")
} else if(muskID=='C'){ } else if(muskID=='c'){
p11$layers[[4]]$data$label <- c('italic(hoxc10a)',"italic(hoxc11a)","italic(hoxc12a)","italic(hoxc13a)","italic(calcoco1)") p11$layers[[4]]$data$label <- c('italic(hoxc10a)',"italic(hoxc11a)","italic(hoxc12a)","italic(hoxc13a)","italic(calcoco1)")
}else if(muskID=='E'){ }else if(muskID=='e'){
p11$layers[[4]]$data$label <- c("italic(polr1d)","italic(ednrb)") p11$layers[[4]]$data$label <- c("italic(polr1d)","italic(ednrb)")
} else if(muskID=='F'){ } else if(muskID=='f'){
p11$layers[[4]]$data$label <- c("italic(foxd3)","italic(alg6)","italic(efcab7)") p11$layers[[4]]$data$label <- c("italic(foxd3)","italic(alg6)","italic(efcab7)")
} else if(muskID=='D'){ } else if(muskID=='d'){
p11$layers[[4]]$data$label <- c("italic(hcfc1)","italic(hcfc1[2])","italic(hcfc1[1])",'italic(paste(sws2a,"\u03B2"))', p11$layers[[4]]$data$label <- c("italic(hcfc1)","italic(hcfc1[2])","italic(hcfc1[1])",'italic(paste(sws2a,"\u03B2"))',
'italic(paste(sws2a,"\u03B1"))',"italic(sw2b)","italic(lws)", 'italic(paste(sws2a,"\u03B1"))',"italic(sw2b)","italic(lws)",
"italic(gnl3l)","italic(tfe3)","italic(mdfic2)","italic(cxxc1[3])", "italic(gnl3l)","italic(tfe3)","italic(mdfic2)","italic(cxxc1[3])",
"italic(cxxc1[1])",'italic(mbd1)','italic(ccdc120)') "italic(cxxc1[1])",'italic(mbd1)','italic(ccdc120)')
}else if(muskID=='G'){ }else if(muskID=='g'){
p11$layers[[4]]$data$label <- c("italic(lgals3bpb)","italic(mpnd)","italic(sh3gl1)","italic(rorb)") p11$layers[[4]]$data$label <- c("italic(lgals3bpb)","italic(mpnd)","italic(sh3gl1)","italic(rorb)")
}else if(muskID=='H'){ }else if(muskID=='h'){
p11$layers[[4]]$data$label <- c("italic(alg2)","italic(sec61b)","italic(nr4a3)", p11$layers[[4]]$data$label <- c("italic(alg2)","italic(sec61b)","italic(nr4a3)",
"italic(stx17)","italic(erp44)","italic(invs)") "italic(stx17)","italic(erp44)","italic(invs)")
} }
......
se <- function(x){
x2 <- x[!is.na(x)]
return(sd(x2)/sqrt(length(x2)))}
trplot <- function(sel){
files <- c("global","boc","bel","hon","pue","nig","uni","nig-pue","nig-uni","pue-uni")
data <- read.csv(paste('../../2_output/08_popGen/07_LD/',files[sel],'-10000-bins.txt',sep=''),sep='\t')
message(files[sel])
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)
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)])
Gcol <- '#3bb33b'
Zcol = rgb(.94,.94,.94)
DG <- rgb(.4,.4,.4)
LGoffset <- 15
GLABoffset <- 8
if(sel %in% c(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))
textSCALE1 <- c(1.8,rep(NA,6),.7)
textSCALE2 <- c(2,rep(NA,6),1)
textSCALE3 <- c(3.5,rep(NA,6),1.75)
p1 <- ggplot(dt %>% filter(!is.na(Mval)),aes(fill=Mval))+
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=textSCALE1[sel])+
geom_text(inherit.aes = F,data=genes,
aes(x=labPOS+GLABoffset,y=labPOS-GLABoffset,label=name),
angle=-45,fontface='italic',size=textSCALE2[sel])+
geom_text(inherit.aes = F,data=zmLab,
aes(x=Nx+LGoffset-.8,y=Nx-LGoffset+.8,label=label),
angle=-45,fontface='bold',size=textSCALE3[sel],col='white')+
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))+
theme_void()+
theme(legend.position = c(.9,.75),
legend.direction = 'vertical')
} else {
p1 <- ggplot(dt %>% filter(!is.na(Mval)),aes(fill=Mval))+
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')+
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)
}
\ No newline at end of file
create_K_plot <- function(searchLG,gfffile,xr,xr2,muskID){ create_K_plot <- function(searchLG,gfffile,xr,xr2,muskID){
source('../../0_data/0_scripts/S07.functions.R'); source('../../0_data/0_scripts/S07.functions.R');
source('../../0_data/0_scripts/E1.getFSTs.R') source('../../0_data/0_scripts/F3.getFSTs.R')
highclr <- '#3bb33b' highclr <- '#3bb33b'
theme_set(theme_minimal(base_size = 5)) theme_set(theme_minimal(base_size = 5))
......
...@@ -2,11 +2,11 @@ trplot <- function(sel){ ...@@ -2,11 +2,11 @@ trplot <- function(sel){
files <- c("global","boc","bel","hon","pue","nig","uni","nig-pue","nig-uni","pue-uni") files <- c("global","boc","bel","hon","pue","nig","uni","nig-pue","nig-uni","pue-uni")
data <- read.csv(paste('../../2_output/08_popGen/08_LD_extended/',files[sel],'.ext-10000-bins.txt',sep=''),sep='\t') data <- read.csv(paste('../../2_output/08_popGen/08_LD_extended/',files[sel],'.ext-10000-bins.txt',sep=''),sep='\t')
message(files[sel]) message(files[sel])
s1=5760000;s2=1884000;s3=12860000;s4=17745000;s5=20149000;s6=22085000;s7=22445000;s8=13730000 s1=5760000;s2=1884000;s3=12860000;s4=17745000;s5=20149000;s6=22085000;s7=22445000;s8=13730000
e2=2134000;e3=13110000;e5=20399000;e6=22335000 e2=2134000;e3=13110000;e5=20399000;e6=22335000
stp=20000;lALL=250000 stp=20000;lALL=250000
scling <- c(-s1, scling <- c(-s1,
-s2+(lALL+stp), -s2+(lALL+stp),
-s3+((lALL+stp)*2), -s3+((lALL+stp)*2),
...@@ -15,7 +15,7 @@ trplot <- function(sel){ ...@@ -15,7 +15,7 @@ trplot <- function(sel){
-s6+((lALL+stp)*5), -s6+((lALL+stp)*5),
-s7+((lALL+stp)*6), -s7+((lALL+stp)*6),
-s8+((lALL+stp)*7)) -s8+((lALL+stp)*7))
dt <- data %>% mutate(miX = floor(Mx/10000),miY=floor(My/10000)) 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), 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), end=c(5878614,2010303,12976894,17873915,20347811,22225149,22555052,22559742,22566321,22575503,13884003),
...@@ -25,7 +25,7 @@ trplot <- function(sel){ ...@@ -25,7 +25,7 @@ trplot <- function(sel){
genes <- genes %>% mutate(Nx1 = (start+scling[sclr])/10000,Nx2 = (end+scling[sclr])/10000, genes <- genes %>% mutate(Nx1 = (start+scling[sclr])/10000,Nx2 = (end+scling[sclr])/10000,
labPOS = (Nx1+Nx2)/2) 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,-3,5,9) 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,-3,5,9)
GZrS <- 40000;GZrS2 <- 20000; #gene zoom offset GZrS <- 40000;GZrS2 <- 20000; #gene zoom offset
BS <- c(-2,2, # ednrb BS <- c(-2,2, # ednrb
2,-2, # foxd3 2,-2, # foxd3
...@@ -61,11 +61,11 @@ trplot <- function(sel){ ...@@ -61,11 +61,11 @@ trplot <- function(sel){
genes$start[9],genes$end[9],genes$end[9]-GZrS+BS[17],genes$start[9]-GZrS+BS[18],genes$start[9], genes$start[9],genes$end[9],genes$end[9]-GZrS+BS[17],genes$start[9]-GZrS+BS[18],genes$start[9],
genes$start[10],genes$end[10],genes$end[10]-GZrS+BS[19]+OPSoffset,genes$start[10]-GZrS+BS[20]+OPSoffset,genes$start[10], genes$start[10],genes$end[10],genes$end[10]-GZrS+BS[19]+OPSoffset,genes$start[10]-GZrS+BS[20]+OPSoffset,genes$start[10],
genes$start[11],genes$end[11],genes$end[11]-GZrS+BS[21],genes$start[11]-GZrS+BS[22],genes$start[11])-GZrS2, genes$start[11],genes$end[11],genes$end[11]-GZrS+BS[21],genes$start[11]-GZrS+BS[22],genes$start[11])-GZrS2,
grp=rep(letters[1:11],each=5)) %>% grp=rep(letters[1:11],each=5)) %>%
mutate(sclr=rep(c(1,2,3,4,5,6,7,7,7,7,8),each=5), mutate(sclr=rep(c(1,2,3,4,5,6,7,7,7,7,8),each=5),
Nx1 = (x+scling[sclr])/10000, Nx1 = (x+scling[sclr])/10000,
Nx2 = (y+scling[sclr])/10000) Nx2 = (y+scling[sclr])/10000)
clr = c(viridis::inferno(5)[c(1,1:5)]) clr = c(viridis::inferno(5)[c(1,1:5)])
Gcol <- '#3bb33b' Gcol <- '#3bb33b'
Zcol = rgb(.94,.94,.94) Zcol = rgb(.94,.94,.94)
...@@ -89,7 +89,7 @@ trplot <- function(sel){ ...@@ -89,7 +89,7 @@ trplot <- function(sel){
s6,s6+lALL,s6+lALL-rS,s6-rS,s6, s6,s6+lALL,s6+lALL-rS,s6-rS,s6,
s7,s7+lALL,s7+lALL-rS,s7-rS,s7, s7,s7+lALL,s7+lALL-rS,s7-rS,s7,
s8,s8+lALL,s8+lALL-rS,s8-rS,s8), s8,s8+lALL,s8+lALL-rS,s8-rS,s8),
grp=rep(letters[1:8],each=5)) %>% grp=rep(letters[1:8],each=5)) %>%
mutate(sclr=rep(1:8,each=5), mutate(sclr=rep(1:8,each=5),
Nx1 = (x+scling[sclr])/10000-1, Nx1 = (x+scling[sclr])/10000-1,
Nx2 = (y+scling[sclr])/10000) Nx2 = (y+scling[sclr])/10000)
...@@ -110,11 +110,11 @@ trplot <- function(sel){ ...@@ -110,11 +110,11 @@ trplot <- function(sel){
s6,s6+lALL,s6+lALL-rS2,s6-rS2,s6, s6,s6+lALL,s6+lALL-rS2,s6-rS2,s6,
s7,s7+lALL,s7+lALL-rS2,s7-rS2,s7, s7,s7+lALL,s7+lALL-rS2,s7-rS2,s7,
s8,s8+lALL,s8+lALL-rS2,s8-rS2,s8), s8,s8+lALL,s8+lALL-rS2,s8-rS2,s8),
grp=rep(letters[1:8],each=5)) %>% grp=rep(letters[1:8],each=5)) %>%
mutate(sclr=rep(1:8,each=5), mutate(sclr=rep(1:8,each=5),
Nx1 = (x+scling[sclr])/10000-1, Nx1 = (x+scling[sclr])/10000-1,
Nx2 = (y+scling[sclr])/10000) Nx2 = (y+scling[sclr])/10000)
rS3 <- .07*rS rS3 <- .07*rS
zmRange3 <- data.frame(x=c(s1,s1+lALL,s1+lALL+rS3,s1+rS3,s1, zmRange3 <- data.frame(x=c(s1,s1+lALL,s1+lALL+rS3,s1+rS3,s1,
s2,s2+lALL,s2+lALL+rS3,s2+rS3,s2, s2,s2+lALL,s2+lALL+rS3,s2+rS3,s2,
...@@ -132,16 +132,16 @@ trplot <- function(sel){ ...@@ -132,16 +132,16 @@ trplot <- function(sel){
s6,s6+lALL,s6+lALL-rS3,s6-rS3,s6, s6,s6+lALL,s6+lALL-rS3,s6-rS3,s6,
s7,s7+lALL,s7+lALL-rS3,s7-rS3,s7, s7,s7+lALL,s7+lALL-rS3,s7-rS3,s7,
s8,s8+lALL,s8+lALL-rS3,s8-rS3,s8), s8,s8+lALL,s8+lALL-rS3,s8-rS3,s8),
grp=rep(letters[1:8],each=5)) %>% grp=rep(letters[1:8],each=5)) %>%
mutate(sclr=rep(1:8,each=5), mutate(sclr=rep(1:8,each=5),
Nx1 = (x+scling[sclr])/10000-1, Nx1 = (x+scling[sclr])/10000-1,
Nx2 = (y+scling[sclr])/10000) Nx2 = (y+scling[sclr])/10000)
zmLab <- data.frame(x=c(s1+.5*lALL,s2+.5*lALL,s3+.5*lALL,s4+.5*lALL, zmLab <- data.frame(x=c(s1+.5*lALL,s2+.5*lALL,s3+.5*lALL,s4+.5*lALL,
s5+.5*lALL,s6+.5*lALL,s7+.5*lALL,s8+.5*lALL), s5+.5*lALL,s6+.5*lALL,s7+.5*lALL,s8+.5*lALL),
label=c("LG4 (E)","LG08 (F)","LG08 (G)",'LG09 (A)','LG12 (B)','LG12 (C)','LG17 (D)',"LG20 (H)")) %>% label=c("LG4 (e)","LG08 (f)","LG08 (g)",'LG09 (a)','LG12 (b)','LG12 (c)','LG17 (d)',"LG20 (h)")) %>%
mutate(sclr=1:8, Nx= (x+scling[sclr])/10000) mutate(sclr=1:8, Nx= (x+scling[sclr])/10000)
zmEND <- data.frame(x=c(s1,s1+lALL, zmEND <- data.frame(x=c(s1,s1+lALL,
s2,s2+lALL, s2,s2+lALL,
s3,s3+lALL, s3,s3+lALL,
...@@ -149,17 +149,17 @@ trplot <- function(sel){ ...@@ -149,17 +149,17 @@ trplot <- function(sel){
s5,s5+lALL, s5,s5+lALL,
s6,s6+lALL, s6,s6+lALL,
s7,s7+lALL, s7,s7+lALL,
s8,s8+lALL)) %>% s8,s8+lALL)) %>%
mutate(sclr=rep(1:8,each=2), mutate(sclr=rep(1:8,each=2),
Nx = ((x+scling[sclr])/10000)+rep(c(2.5,-4),8), Nx = ((x+scling[sclr])/10000)+rep(c(2.5,-4),8),
label=round((x/1000000),2)) label=round((x/1000000),2))
DG <- rgb(.4,.4,.4) DG <- rgb(.4,.4,.4)
gnShift <- c(rep(0,6),3.5,0,0,3.5,0) gnShift <- c(rep(0,6),3.5,0,0,3.5,0)
textSCALE1 <- c(1.8,rep(NA,6),.9) textSCALE1 <- c(1.8,rep(NA,6),.9)
textSCALE2 <- c(2,rep(NA,6),1) textSCALE2 <- c(2,rep(NA,6),1)
textSCALE3 <- c(2.5,rep(NA,6),1.25) textSCALE3 <- c(2.5,rep(NA,6),1.25)
p1 <- ggplot(dt %>% filter(!is.na(Mval)),aes(fill=Mval))+ p1 <- ggplot(dt %>% filter(!is.na(Mval)),aes(fill=Mval))+
coord_equal()+ coord_equal()+
geom_polygon(inherit.aes = F,data=zmRange, geom_polygon(inherit.aes = F,data=zmRange,
......
create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchsnp,muskID){ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchsnp,muskID){
source('../../0_data/0_scripts/E1.functions.R'); source('../../0_data/0_scripts/F3.functions.R');
source('../../0_data/0_scripts/E1.getFSTs.R') source('../../0_data/0_scripts/F3.getFSTs.R')
source('../../0_data/0_scripts/E1.getGxP.R') source('../../0_data/0_scripts/F3.getGxP.R')
source('../../0_data/0_scripts/E1.getDXY.R') source('../../0_data/0_scripts/F3.getDXY.R')
source('../../0_data/0_scripts/S15.getPI.R') source('../../0_data/0_scripts/S15.getPI.R')
source('../../0_data/0_scripts/S15.getPIpw.R') source('../../0_data/0_scripts/S15.getPIpw.R')
source('../../0_data/0_scripts/S15.getIHH12.R') source('../../0_data/0_scripts/S15.getIHH12.R')
......
This diff is collapsed.
--- ---
output: html_document output: html_document
editor_options: editor_options:
chunk_output_type: console chunk_output_type: console
--- ---
# Figure 1 # Figure 1
```{r setup, include=FALSE} ```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = './F_scripts') knitr::opts_knit$set(root.dir = './F_scripts')
...@@ -49,7 +49,7 @@ library(marmap) ...@@ -49,7 +49,7 @@ library(marmap)
source("../../0_data/0_scripts/F1.functions.R") source("../../0_data/0_scripts/F1.functions.R")
``` ```
### Fig. 1 a) ### Fig. 1 A)
First, we load the data sheet containing the sampling location as well as the population specifics for each sample and sub sample it to include only the samples included in the re-sequencing study. First, we load the data sheet containing the sampling location as well as the population specifics for each sample and sub sample it to include only the samples included in the re-sequencing study.
...@@ -62,7 +62,7 @@ dataAll <- read.csv('../../0_data/0_resources/F1.sample.txt',sep='\t') %>% mutat ...@@ -62,7 +62,7 @@ dataAll <- read.csv('../../0_data/0_resources/F1.sample.txt',sep='\t') %>% mutat
data <- dataAll %>% filter(sample=="sample") 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)), 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))) %>% 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)) %>% 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')) as.data.frame() %>% cbind(.,read.csv('../../0_data/0_resources/F1.pointer.csv'))
``` ```
...@@ -146,9 +146,9 @@ Now, we got everything needed for the base map (Fig. 1 a): ...@@ -146,9 +146,9 @@ Now, we got everything needed for the base map (Fig. 1 a):
```{r mapBase, results='hide', message=FALSE, warning=FALSE} ```{r mapBase, results='hide', message=FALSE, warning=FALSE}
p1 <- ggplot()+ p1 <- ggplot()+
# set coordinates # set coordinates
coord_map(projection = 'mercator',xlim=xlimW,ylim=ylimW)+ coord_map(projection = 'mercator',xlim=xlimW,ylim=ylimW)+
# axes labels # axes labels
labs(y="Latitude",x="Longitude") + labs(y="Latitude",x="Longitude") +
# preset theme # preset theme
theme_mapK + theme_mapK +
# landmass layer # landmass layer
...@@ -178,17 +178,17 @@ p1 ...@@ -178,17 +178,17 @@ p1
``` ```
</center> </center>
### Fig. 1 b) ### Fig. 1 B)
Next we create the sub figure containing the *F~ST~* information (Fig. 1 b). Next we create the sub figure containing the *F~ST~* information (Fig. 1 B).
We first read in the data and define the color palette. We first read in the data and define the color palette.
```{r fstData, results='hide', message=FALSE, warning=FALSE} ```{r fstData, results='hide', message=FALSE, warning=FALSE}
fstdata <- read.csv('../../2_output/08_popGen/05_fst/genome_wide_weighted_mean_fst.txt', 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')) %>% sep='\t',skip = 1,head=F,col.names = c('pair','fst')) %>%
filter(!pair %in% c('NU','PN','PU')) %>% filter(!pair %in% c('NU','PN','PU')) %>%
arrange(fst) %>% arrange(fst) %>%
mutate(x=row_number()) %>% mutate(x=row_number()) %>%
rowwise() %>% rowwise() %>%
mutate(loc = substr(as.character(pair),2,2), mutate(loc = substr(as.character(pair),2,2),
pop1 = substr(as.character(pair),1,1), pop1 = substr(as.character(pair),1,1),
...@@ -230,7 +230,7 @@ p5 ...@@ -230,7 +230,7 @@ p5
``` ```
</center> </center>
### Fig. 1 c) ### Fig. 1 C)
The last sub figure contains the PCA output for each location. The last sub figure contains the PCA output for each location.
For those, we define the color palette for hamlet species and read in the PCA data generated by `pcadapt`. For those, we define the color palette for hamlet species and read in the PCA data generated by `pcadapt`.
...@@ -257,7 +257,7 @@ Further, the explained variance of each PC is extracted for the axis labels. ...@@ -257,7 +257,7 @@ Further, the explained variance of each PC is extracted for the axis labels.
dataHON <- cbind(dataAll %>% filter(sample=='sample',loc=='hon')%>% select(id,spec), dataHON <- cbind(dataAll %>% filter(sample=='sample',loc=='hon')%>% select(id,spec),
honPCA$scores); names(dataHON)[3:12]<- paste('PC',1:10,sep='') honPCA$scores); names(dataHON)[3:12]<- paste('PC',1:10,sep='')
exp_varHON <- (honPCA$singular.values[1:10])^2/length(honPCA$maf) exp_varHON <- (honPCA$singular.values[1:10])^2/length(honPCA$maf)
xlabHON <- paste('PC1 (',sprintf("%.1f",exp_varHON[1]*100),'%)')# explained varinace)') xlabHON <- paste('PC1 (',sprintf("%.1f",exp_varHON[1]*100),'%)')# explained varinace)')
ylabHON <- paste('PC2 (',sprintf("%.1f",exp_varHON[2]*100),'%)')# explained varinace)') ylabHON <- paste('PC2 (',sprintf("%.1f",exp_varHON[2]*100),'%)')# explained varinace)')
# Honduras PCA (plot) # Honduras PCA (plot)
...@@ -355,11 +355,11 @@ F1 <- ggdraw()+ ...@@ -355,11 +355,11 @@ F1 <- ggdraw()+
draw_grob(nigGrob, xf[1]+xS[1], yf[1]+yS[1], hSC[1], hSC[2])+ # Honduras draw_grob(nigGrob, xf[1]+xS[1], yf[1]+yS[1], hSC[1], hSC[2])+ # Honduras
draw_grob(pueGrob, xf[2]+xS[1], yf[2]+yS[1], hSC[1], hSC[2])+ # Honduras draw_grob(pueGrob, xf[2]+xS[1], yf[2]+yS[1], hSC[1], hSC[2])+ # Honduras
draw_grob(uniGrob, xf[3]+xS[1], yf[3]+yS[1], hSC[1], hSC[2])+ # Honduras draw_grob(uniGrob, xf[3]+xS[1], yf[3]+yS[1], hSC[1], hSC[2])+ # Honduras
draw_grob(nigGrob, xf[1]+xS[2], yf[1]+yS[2], hSC[1], hSC[2])+ # Panama draw_grob(nigGrob, xf[1]+xS[2], yf[1]+yS[2], hSC[1], hSC[2])+ # Panama
draw_grob(pueGrob, xf[2]+xS[2], yf[2]+yS[2], hSC[1], hSC[2])+ # Panama draw_grob(pueGrob, xf[2]+xS[2], yf[2]+yS[2], hSC[1], hSC[2])+ # Panama
draw_grob(uniGrob, xf[3]+xS[2], yf[3]+yS[2], hSC[1], hSC[2])+ # Panama draw_grob(uniGrob, xf[3]+xS[2], yf[3]+yS[2], hSC[1], hSC[2])+ # Panama
draw_grob(nigGrob, xf[1]+xS[3], yf[1]+yS[3], hSC[1], hSC[2])+ # Belize draw_grob(nigGrob, xf[1]+xS[3], yf[1]+yS[3], hSC[1], hSC[2])+ # Belize
draw_grob(pueGrob, xf[2]+xS[3], yf[2]+yS[3], hSC[1], hSC[2])+ # Belize draw_grob(pueGrob, xf[2]+xS[3], yf[2]+yS[3], hSC[1], hSC[2])+ # Belize
draw_grob(uniGrob, xf[3]+xS[3], yf[3]+yS[3], hSC[1], hSC[2])+ # Belize draw_grob(uniGrob, xf[3]+xS[3], yf[3]+yS[3], hSC[1], hSC[2])+ # Belize
...@@ -367,21 +367,21 @@ F1 <- ggdraw()+ ...@@ -367,21 +367,21 @@ F1 <- ggdraw()+
draw_grob(panGrob, 0.3, 0.485, 0.042, 0.042)+ draw_grob(panGrob, 0.3, 0.485, 0.042, 0.042)+
draw_grob(belGrob, 0.2, 0.67, 0.042, 0.042)+ draw_grob(belGrob, 0.2, 0.67, 0.042, 0.042)+
draw_grob(honGrob, 0.28, 0.61, 0.042, 0.042)+ draw_grob(honGrob, 0.28, 0.61, 0.042, 0.042)+
draw_grob(honGrob, 0.285+.09, 0.37, 0.042, 0.042)+ draw_grob(honGrob, 0.285+.09, 0.37, 0.042, 0.042)+
draw_grob(belGrob, 0.535+.09, 0.37, 0.042, 0.042)+ draw_grob(belGrob, 0.535+.09, 0.37, 0.042, 0.042)+
draw_grob(panGrob, 0.785+.09, 0.37, 0.042, 0.042)+ draw_grob(panGrob, 0.785+.09, 0.37, 0.042, 0.042)+
# annotations (subfigure labels) # annotations (subfigure labels)
draw_plot_label(x = c(.005,.005,.25), draw_plot_label(x = c(.005,.005,.25),
y = c(.995,.42,.42), y = c(.995,.42,.42),
label =letters[1:3]) label = LETTERS[1:3])
``` ```
It is then exported using `ggsave()`. It is then exported using `ggsave()`.
```{r exportF1, eval=FALSE} ```{r exportF1, eval=FALSE}
ggsave(plot = F1,filename = '../output/F1.pdf',width = 183,height = 145,units = 'mm',device = cairo_pdf) ggsave(plot = F1,filename = '../output/F1.pdf',width = 183,height = 145,units = 'mm',device = cairo_pdf)
``` ```