include fasteprr

0_data/0_scripts/S14.functions.R → 0_data/0_scripts/E1.functions.R

 get_anno_df_single_line <- function(searchLG,gfffile,
                                     xrange,
-                                    anno_rown=4){
+                                    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))
@@ -13,13 +15,16 @@ get_anno_df_single_line <- function(searchLG,gfffile,
            pe=ifelse(strand=='-',checkStart,checkEnd),
            labelx=mean(c(sort(c(xrange[1],ps))[2],
                          sort(c(xrange[2],pe))[1])),
-           window='bold(Gene)') %>% 
+           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),by='Parent',all.x=T) %>% 
+    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))}
+

0_data/0_scripts/F3.plot_fun.R → 0_data/0_scripts/E1.plot_fun.R

 create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchsnp,muskID){
-  source('../../0_data/0_scripts/F3.functions.R');
-  source('../../0_data/0_scripts/F3.getFSTs.R')
-  source('../../0_data/0_scripts/F3.getGxP.R')
-  source('../../0_data/0_scripts/F3.getDXY.R')
-  
+  source('../../0_data/0_scripts/E1.functions.R');
+  source('../../0_data/0_scripts/E1.getFSTs.R')
+  source('../../0_data/0_scripts/E1.getGxP.R')
+  source('../../0_data/0_scripts/E1.getDXY.R')
+
   highclr <- '#3bb33b'
   theme_set(theme_minimal(base_size = 6))
   if(xr[1]>xr[2]){print('error: ill conditioned x-range')}
   wdth <- .3
-  
-  
+
+
   # get gene models
   df_list <- get_anno_df_single_line(searchLG=searchLG,gfffile=gfffile,xrange=xr,
                                      genes_of_interest = searchgene,genes_of_sec_interest=secondary_genes,anno_rown=3)
-  
+
   # get fst values
   fst_list <- getFSTS(searchLG,xr,searchsnp,highclr)
-  # data_fst_pw<-fst_list$data_fst_pw; 
-  global_fst<-fst_list$global_fst; 
+  # data_fst_pw<-fst_list$data_fst_pw;
+  global_fst<-fst_list$global_fst;
   data_fst<-fst_list$data_fst
-  
+
   # get pfst values
   pfst_list <- getGxP(searchLG,xr)
-  #  data_pfst_pw<-pfst_list$data_pfst_pw; 
+  #  data_pfst_pw<-pfst_list$data_pfst_pw;
   data_pfst<-pfst_list$data_pfst
-  
+
   # get dxy values
   dxy_list <- getDXY(searchLG,xr)
   #data_dxy_pw<-dxy_list$data_dxy_pw;
   data_dxy<-dxy_list$data_dxy
-  
+
   clr <- c('#fb8620','#1b519c','#d93327')
   annoclr <- c('lightgray',highclr,rgb(.3,.3,.3))[1:3]
-  df_list[[1]] <- df_list[[1]] %>% 
+  df_list[[1]] <- df_list[[1]] %>%
     mutate(label=paste("italic(",tolower(Parentgenename),")",sep='') )
-  
+
   LW <- .3;lS <- 9;tS <- 6
   plotSET <-  theme(rect = element_blank(),
                     text=element_text(size=tS,color='black'),
@@ -52,7 +52,7 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
                     strip.placement = "outside",
                     panel.border = element_blank()
   )
-  
+
   p11 <- ggplot()+coord_cartesian(xlim=xr/1000)+
     geom_rect(data=df_list[[2]],aes(xmin=ps/1000,xmax=pe/1000,ymin=yl-(wdth/2),ymax=yl+(wdth/2),
                                     fill=as.factor(clr),col=as.factor(clr),group=Parent),lwd=.1)+
@@ -66,7 +66,7 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
     facet_grid(window~.,scales='free_y',
                switch = 'y',labeller = label_parsed,as.table = T)+
     # color settings
-    scale_color_manual(values=annoclr,breaks=c("x","y","z"),guide=F)+  
+    scale_color_manual(values=annoclr,breaks=c("x","y","z"),guide=F)+
     scale_fill_manual(values=annoclr,guide=F)+
     scale_x_continuous(name=paste0(searchLG,' (',muskID,', kb)'),expand=c(0,0),position = 'top')+
     scale_y_continuous(breaks = seq(0,.75,length.out = 4))+
@@ -86,11 +86,11 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
           strip.placement = "outside",
           panel.border = element_blank()
     );
-  
+
  if(muskID=='A'){
    p11$layers[[4]]$data$label <- c("italic(sox10)","italic(rnaseh2a)")
  } 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'){
      p11$layers[[4]]$data$label <- c('italic(hoxc10a)',"italic(hoxc11a)","italic(hoxc12a)","italic(hoxc13a)","italic(calcoco1)")
    }else if(muskID=='E'){
@@ -108,7 +108,7 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
      p11$layers[[4]]$data$label <- c("italic(alg2)","italic(sec61b)","italic(nr4a3)",
                                         "italic(stx17)","italic(erp44)","italic(invs)")
    }
-  
+
   p12 <-  ggplot()+coord_cartesian(xlim=xr)+
     geom_point(data=global_fst,aes(x=POS,y=WEIR_AND_COCKERHAM_FST),col=global_fst$clr,size=.2)+
     geom_point(data=global_fst,aes(x=POS,y=WEIR_AND_COCKERHAM_FST),col=global_fst$clr2,size=.3)+
@@ -117,7 +117,7 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
               ,aes(x=POS,y=WEIGHTED_FST,col=run),lwd=LW)+
     #  geom_line(data=(data_fst_pw %>% filter(POS > xr[1],POS<xr[2]))
     #            ,aes(x=POS,y=WEIGHTED_FST,col=run,linetype=group),lwd=1)+
-    scale_color_manual(values=c(clr,annoclr),breaks=c("nig-pue","nig-uni","pue-uni"),guide=F)+  
+    scale_color_manual(values=c(clr,annoclr),breaks=c("nig-pue","nig-uni","pue-uni"),guide=F)+
     scale_fill_manual(values=annoclr,guide=F)+
     facet_grid(window~.,scales='free_y',
                switch = 'y',labeller = label_parsed,as.table = T)+
@@ -125,34 +125,34 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
     scale_y_continuous(breaks = seq(0,.75,length.out = 4))+
     scale_linetype(name='location',label=c('Belize','Honduras','Panama'))+
     guides(linetype= guide_legend(override.aes = list(color = 'black')))+plotSET
-  
+
   p13 <-  ggplot()+coord_cartesian(xlim=xr)+
     geom_line(data=(data_pfst %>% filter(POS > xr[1],POS<xr[2]))
               ,aes(x=POS,y=avgp_wald,col=run),lwd=LW)+
     # geom_line(data=(data_pfst_pw %>% filter(POS > xr[1],POS<xr[2]))
     #            ,aes(x=POS,y=avgp_wald,col=run,linetype=group),lwd=1)+
     scale_color_manual(values=c(clr,annoclr),
-                       breaks=c("nig-pue","nig-uni","pue-uni","x","y","z"),guide=F)+  
+                       breaks=c("nig-pue","nig-uni","pue-uni","x","y","z"),guide=F)+
     facet_grid(window~.,scales='free_y',
                switch = 'y',labeller = label_parsed,as.table = T)+
     scale_x_continuous(name=searchLG,expand=c(0,0),position = 'top')+
     scale_y_continuous(breaks = seq(0,9,length.out = 4))+
     scale_linetype(name='location',label=c('Belize','Honduras','Panama'))+
     guides(linetype= guide_legend(override.aes = list(color = 'black')))+plotSET
-  
+
   p14 <-  ggplot()+coord_cartesian(xlim=xr)+
     geom_line(data=(data_dxy %>% filter(POS > xr[1],POS<xr[2]))
               ,aes(x=POS,y=dxy,col=run),lwd=LW)+
     # geom_line(data=(data_dxy_pw %>% filter(POS > xr[1],POS<xr[2]))
     #            ,aes(x=POS,y=dxy,col=run,linetype=group),lwd=1)+
-    scale_color_manual(values=c(clr,annoclr),breaks=c("nig-pue","nig-uni","pue-uni","x","y","z"),guide=F)+  
+    scale_color_manual(values=c(clr,annoclr),breaks=c("nig-pue","nig-uni","pue-uni","x","y","z"),guide=F)+
     facet_grid(window~.,scales='free_y',
                switch = 'y',labeller = label_parsed,as.table = T)+
     scale_x_continuous(name=searchLG,expand=c(0,0),position = 'top')+
     scale_y_continuous(breaks = seq(0,.01,length.out = 5))+
     scale_linetype(name='location',label=c('Belize','Honduras','Panama'))+
     guides(linetype= guide_legend(override.aes = list(color = 'black')))+plotSET
-  
+
   p2 <- plot_grid(p11,p12,p13,p14,
                   ncol = 1,align = 'v',axis = 'r',
                   rel_heights = c(1.3,rep(1,3)))

0_data/0_scripts/F4.functions.R

-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

0_data/0_scripts/S02.functions.R

+substrAllRight <- function(x, n){substr(x, n+1, nchar(x))}
+substrRight <- function(x, n){substr(x, nchar(x)-n+1, nchar(x))}
+
+theme_mapK <- theme_minimal()+
+  theme(panel.grid = element_blank(),
+        text=element_text(size=10),
+        axis.text.y = element_text(size = 7),
+        axis.text.x = element_text(size = 7))
+
+GeomCustom <- ggproto(
+  "GeomCustom",
+  Geom,
+  setup_data = function(self, data, params) {
+    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,
+                         data = NULL,
+                         stat = "identity",
+                         position = "identity",
+                         na.rm = FALSE,
+                         show.legend = NA,
+                         inherit.aes = FALSE,
+                         ...) {
+  layer(
+    geom = GeomCustom,
+    mapping = mapping,
+    data = data,
+    stat = stat,
+    position = position,
+    show.legend = show.legend,
+    inherit.aes = inherit.aes,
+    params = list(na.rm = na.rm, ...)
+  )
+}
+
+
+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 %>%
+    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',
+                            '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"))
+
+  groblist <- tibble(RUN=gsub('XX',loc,c('nigXX-pueXX','nigXX-uniXX','pueXX-uniXX')),
+                     grob = gl )
+
+  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~.)+
+    geom_hline(yintercept=c(0,.5,1),col=rgb(.9,.9,.9))+
+    geom_custom(data=groblist, aes(grob=grob), inherit.aes = FALSE,
+                x = 1, y=0.535)+
+    geom_point(inherit.aes = F, data=dotdat,
+               aes(x=x,y=-.09,col=spec),size=2)+
+    geom_point(inherit.aes = F, data=dotdat,
+               aes(x=x,y=-.09),col='black',size=2,shape=1)+
+    geom_bar(aes(y=post_prob),stat = 'identity',position = "stack")+
+    scale_fill_manual('class',values=clr)+
+    scale_color_manual(values = c('#000000','#d45500','#ffffff'))+
+    scale_y_continuous('',breaks = c(0,.5,1))+
+    scale_x_continuous(breaks=1:(Grobx-1)-.5,labels=XL[[sel]])+
+    theme_minimal()+
+    theme(panel.grid = element_blank(),
+          axis.text.y = element_text(size = 7),
+          axis.text.x = element_text(angle = 90,size = 7),
+          axis.title.x = element_blank(),
+          legend.position = 'none',
+          strip.background = element_blank(),
+          strip.text = element_blank())
+
+  return(p)
+}
+
+plotPCA <- function(pca,dataAll,locIN,introgression_candidates,clr,fll){
+  data <- cbind(dataAll %>% filter(sample=='sample',loc==locIN) %>% 
+                  select(id,spec),
+                  pca$scores);
+  names(data)[3:12]<- paste('PC',1:10,sep='')
+  exp_var <- (pca$singular.values[1:10])^2/length(pca$maf)
+  xlab <- paste('PC1 (',sprintf("%.1f",exp_var[1]*100),'%)')# explained varinace)') 
+  ylab <- paste('PC2 (',sprintf("%.1f",exp_var[2]*100),'%)')# explained varinace)')
+  p <- ggplot(data,aes(x=PC1,y=PC2,col=spec,fill=spec))+
+    geom_point(size=1.1,shape=21)+
+    geom_text_repel(inherit.aes = F,
+                    data=(data %>% filter(gsub('[a-z,A-Z]','',id) %in% introgression_candidates)),
+                    aes(x=PC1,y=PC2,col=spec,label=substr(id,1,5)),angle=30,size=2.5,nudge_y = .01)+
+    geom_point(inherit.aes = F,data=(data %>% filter(gsub('[a-z,A-Z]','',id) %in% introgression_candidates)),
+               aes(x=PC1,y=PC2),col=rgb(.5,.5,.5),size=2.5,shape=1)+
+    scale_color_manual(values=clr,guide=F)+
+    scale_fill_manual(values=fll,guide=F)+
+    theme_mapK+
+    theme(legend.position='bottom',
+          panel.border = element_rect(color=rgb(.9,.9,.9),fill=rgb(1,1,1,0)),
+          axis.title.y = element_text(vjust = -7),
+          axis.title.x = element_text(vjust = 6),
+          axis.text = element_blank(),
+          plot.margin = unit(c(10,5,10,30),'pt'))+#coord_equal()+
+    scale_x_continuous(name=xlab,breaks = c(-.2,.5),position = "top")+
+    scale_y_continuous(name=ylab,breaks = c(-.4,.2),position = "right");
+  return(p)
+}
\ No newline at end of file

0_data/0_scripts/S14.plot_fun.R → 0_data/0_scripts/S07.plot_fun.R

 create_K_plot <- function(searchLG,gfffile,xr,xr2,muskID){
-  source('../../0_data/0_scripts/S14.functions.R');
-  source('../../0_data/0_scripts/F3.getFSTs.R')
-  
+  source('../../0_data/0_scripts/S07.functions.R');
+  source('../../0_data/0_scripts/E1.getFSTs.R')
+
   highclr <- '#3bb33b'
-  theme_set(theme_minimal(base_size = 6))
+  theme_set(theme_minimal(base_size = 5))
   if(xr[1]>xr[2]){print('error: ill conditioned x-range')}
   wdth <- .3
-  
-  
+
+
   # get gene models
   df_list <- get_anno_df_single_line(searchLG=searchLG,gfffile=gfffile,xrange=xr,anno_rown=4)
-  
+
   # get fst values
   fst_list <- getFSTS(searchLG,xr,c(),highclr)
-  # data_fst_pw<-fst_list$data_fst_pw; 
-  global_fst<-fst_list$global_fst; 
+  # data_fst_pw<-fst_list$data_fst_pw;
+  global_fst<-fst_list$global_fst;
   data_fst<-fst_list$data_fst
-  
+
   clr <- c('#fb8620','#1b519c','#d93327')
   annoclr <- c('lightgray',highclr,rgb(.3,.3,.3))[1:3]
-  df_list[[1]] <- df_list[[1]] %>% 
+  df_list[[1]] <- df_list[[1]] %>%
     rowwise() %>%
     mutate(label=unlist(strsplit(tolower(Parentgenename), "_"))[1])
-  
-  LW <- .3;lS <- 9;tS <- 6
+
+  LW <- .3;lS <- 7;tS <- 5
   plotSET <-  theme(rect = element_blank(),
                     text=element_text(size=tS,color='black'),
                     panel.grid.major = element_line(colour = rgb(.92,.92,.92)),
@@ -43,7 +43,7 @@ create_K_plot <- function(searchLG,gfffile,xr,xr2,muskID){
   Gclr=rgb(.3,.3,.3)
   range_df <- data.frame(x=xr2[1],x2=xr2[2])
   range_clr <- rgb(0,0,.4,.1)
-  
+
   p11 <- ggplot()+coord_cartesian(xlim=xr/1000)+
     geom_rect(data=range_df,aes(xmin = x/1000,xmax=x2/1000,ymin=-Inf,ymax=Inf),
               col=range_clr,fill=range_clr)+
@@ -60,7 +60,7 @@ create_K_plot <- function(searchLG,gfffile,xr,xr2,muskID){
     facet_grid(window~.,scales='free_y',
                switch = 'y',labeller = label_parsed,as.table = T)+
     # color settings
-    #scale_color_manual(values=annoclr,breaks=c("x","y","z"),guide=F)+  
+    #scale_color_manual(values=annoclr,breaks=c("x","y","z"),guide=F)+
     #scale_fill_manual(values=annoclr,guide=F)+
     scale_x_continuous(name=paste0(searchLG,' (',muskID,', kb)'),expand=c(0,0),position = 'top')+
     scale_y_continuous(breaks = seq(0,.75,length.out = 4))+
@@ -91,7 +91,7 @@ create_K_plot <- function(searchLG,gfffile,xr,xr2,muskID){
               ,aes(x=POS,y=WEIGHTED_FST,col=run),lwd=LW)+
     #  geom_line(data=(data_fst_pw %>% filter(POS > xr[1],POS<xr[2]))
     #            ,aes(x=POS,y=WEIGHTED_FST,col=run,linetype=group),lwd=1)+
-    scale_color_manual(values=c(clr,annoclr),breaks=c("nig-pue","nig-uni","pue-uni"),guide=F)+  
+    scale_color_manual(values=c(clr,annoclr),breaks=c("nig-pue","nig-uni"