Commit d7507191 authored by Kosmas Hench's avatar Kosmas Hench
Browse files

update fig 1 documentation

parent 25c7332e
Loading
Loading
Loading
Loading
+384 −1
Original line number Diff line number Diff line
---
output: html_document
editor_options: 
  chunk_output_type: console
---
# F1
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = './F1')
```

## Summary

This is the accessory documentation of Figure 1.
It explains the steps of script **2.2.0.1.maf2circos.sh** and documents progression from the whole genome alignment of the *de novo* *Hypoplectrus puella* genome assembly with the *Gasterosteus aculeatus* reference genome. 
 No newline at end of file
The Figure can be recreated by using runing the **R** script F1.R:
```sh
cd $WORK/3_figures/F1

Rscript --vanilla F1.R
rm Rplots.pdf

```

## Details of F1.R

In the following, the individual steps of the R script are documented. Is an executable R script that depends on a variety of image manipulation and geodata packages.
Aditionally, the supporting R script (**F1.functions.R**) in needed:


```{r head,results='hide',message=FALSE}
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")
```

### Fig. 1 a)

First, we load the data sheet containing the sampling location as well as the population specifics for each sample and subsample it to include only the samples included in the resequencing study.

The Coordinates need to be transformed from (deg min sec) to decimal degrees.
The sample inventory is summarised over populations and merged with the coordinates later used for the sample pointers.


```{r mapSampleData, results='hide', message=FALSE, warning=FALSE}
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'))
```

Then, we read in the coordinates of helper points that will be used to draw the bezier curves of the sample pointers.
We also set the spacial range of the map.

```{r mapSampleAnnotation, results='hide', message=FALSE, warning=FALSE}
crvs <- read.csv('../../0_data/0_resources/F1.curve.csv')
xlimW = c(-100,-53); ylimW = c(7,30.8)
```

The base map relies on the *world* data from the *maps* package.
This data is loaded and clipped to the extend of the plotting area.

```{r mapBaseData, results='hide', message=FALSE, warning=FALSE}
worldmap = map_data("world")
names(worldmap) <- c("X","Y","PID","POS","region","subregion")
worldmap = clipPolys(worldmap, xlim=xlimW,ylim=ylimW, keepExtra=TRUE)
```

External images are loaded for annotation.
This includes the hamlet images, the contry flags and legend.

```{r mapAnnotationImgs, results='hide', message=FALSE, warning=FALSE}
# compas and legend
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"))))
# hamlets
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"))))
# flags
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"))))
# globe and pairwise legend
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"))))
```

We create a summary data table, containing the sample size of each population and the position for the annotation within the map.

```{r mapSampleSizes, results='hide', message=FALSE, warning=FALSE}
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
```

Next, we read in the shappe files containg the distribution of the halet species.
One file contains the overlap of the distribution ranges of *H. nigricans*, *H. puella* and *H. unicolor*, while the other contains the distribution range of the whole genus (*Hypoplectrus*).
These distributions rely on the data provided by the [*Smithsonian Tropical Research Institute*](http://biogeodb.stri.si.edu/caribbean/en/research/index/range) that were processed using [`gimp`](https://www.gimp.org/) and [`qgis`](https://qgis.org/en/site/) (including the [`Georeferencer Plugin`](https://docs.qgis.org/2.8/en/docs/user_manual/plugins/plugins_georeferencer.html))

```{r mapDistributionData, results='hide', message=FALSE, warning=FALSE}
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)
```

We set some colors used for plotting:

```{r mapColors, results='hide', message=FALSE, warning=FALSE}
clrShps <- c("#ffc097","#b1d4f8") # hamlet distributions
cFILL <- rgb(.7,.7,.7) # landmass
ccc<-rgb(0,.4,.8) # sample pointers
```

Now, we got everything needed for the base map (Fig. 1 a):

```{r mapBase, results='hide', message=FALSE, warning=FALSE}
p1 <- ggplot()+
  # set coordinates
  coord_map(projection = 'mercator',xlim=xlimW,ylim=ylimW)+ 
  # axes labels
  labs(y="Latitude",x="Longitude") + 
  # preset theme
  theme_mapK +
  # landmass layer
  geom_polygon(data=worldmap,aes(X,Y,group=PID),fill = cFILL,
               col=rgb(0,0,0),lwd=.2)+
  # genus distribution layer
  geom_polygon(data=hS_f,aes(X,Y,group=PID),fill=clrShps[2],col=rgb(.3,.3,.3),lwd=.2)+
  # species distribution layer
  geom_polygon(data=uS_f,aes(X,Y,group=PID),fill=clrShps[1],col=clrShps[1],lwd=.2)+
  # point background
  geom_point(data=sampSize,aes(x=x,y=y),col='white',size=7)+
  # sample size annotation
  geom_text(data=sampSize,aes(x=x,y=y,label=n),size=3)+
  # sample pointer curves
  geom_bezier(data = crvs,aes(x= x, y = y, group = loc),col=ccc,lwd=.4) +
  # sample pointer points
  geom_point(data = crvs[c(3,6,9),],aes(x= x, y = y, group = loc),size=5,col=ccc,shape=21,fill='white')+
  # sample pointer labels
  geom_text(data = crvs[c(3,6,9),],aes(x= x, y = y, label=c('B','P','H')),size=3)+
  # sample point anchors
  geom_point(data=dataSum,aes(x=w,y=n),size=1.5,shape=21,fill='white')
```

<center>
```{r mapBaseSHOW, echo=FALSE}
p1
```
</center>

### Fig. 1 b)

Next we create the subfigure containing the *F~ST~* information (Fig. 1 b).
We first read in the data and define the color palette.

```{r fstData, results='hide', message=FALSE, warning=FALSE}
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='-'))

# defining the pairwise color palette
clrP <- c('#fb8620','#1b519c','#d93327')
```

And than generate the bar plot:

```{r fstBase, results='hide', message=FALSE, warning=FALSE}
p5 <- ggplot(fstdata,aes(x=x,y=fst,fill=pairclass))+
  # bar layer
  geom_bar(stat='identity')+
  # adding annotations
  annotation_custom(legPGrob,xmin =0.15,xmax =5,ymin =.008,ymax = .035)+
  # set color palette
  scale_fill_manual(values = clrP)+
  # y label and axis
  scale_y_continuous(name=expression(italic("F"['ST'])),
                     expand = c(0,0))+
  # x labels and axis
  scale_x_continuous(expand = c(.01,0),breaks = 1:9,
                     labels = c('H','H','B','H','P','B','B','P','P'))+
  # preset theme
  theme_mapK+
  # theme adjustments
  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(size = 8));
```

<center>
```{r fstBaseSHOW, echo=FALSE}
p5
```
</center>

### Fig. 1 c)

The last subfigure contains the PCA output for each location.
For those, we define the colorpalette for hamlet species and read in the PCA data generated by `pcadapt`.

```{r pcaData, results='hide', message=FALSE, warning=FALSE}

clr<-c('#000000','#d45500','#000000')
fll<-c('#000000','#d45500','#ffffff')
eVclr <- c(rep('darkred',2),rep(cFILL,8))

# Belize data
belPCA <- readRDS('../../2_output/08_popGen/04_pca/bel/belpca.Rds')
# Honduras data
honPCA <- readRDS('../../2_output/08_popGen/04_pca/hon/honpca.Rds')
# Panama data
panPCA <- readRDS('../../2_output/08_popGen/04_pca/pan/panpca.Rds')
```

The data for each PCA is merged with metadata about the samples and than plotted.
Further, the explained variance of each PC is extracted for the axis labels.

```{r honPCA, results='hide', message=FALSE, warning=FALSE}
# Honduras PCA (data)
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)')

# Honduras PCA (plot)
p2 <- ggplot(dataHON,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,shape=21)+
  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),
        plot.margin = unit(rep(-10,4),'pt'))+#coord_equal()+
  scale_x_continuous(name=xlabHON,breaks = c(-.2,.5))+
  scale_y_continuous(name=ylabHON,breaks = c(-.4,.2))
```
---
```{r belPCA, results='hide', message=FALSE, warning=FALSE}
# Belize PCA (data)
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),'%)')

# Belize PCA (plot)
p3 <- ggplot(dataBEL,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,shape=21)+
  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),
        plot.margin = unit(rep(-10,4),'pt'))+
  scale_x_continuous(name=xlabBEL,breaks = c(-.15,.25))+
  scale_y_continuous(name=ylabBEL,breaks = c(-.3,.2))
```
---
```{r panPCA, results='hide', message=FALSE, warning=FALSE}
# Panama PCA (data)
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),'%)')

# Panama PCA (plot)
p4 <- ggplot(dataPAN,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,shape=21)+
  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),
        plot.margin = unit(rep(-10,4),'pt'))+
  scale_x_continuous(name=xlabPAN,breaks = c(-.15,.2))+
  scale_y_continuous(name=ylabPAN,breaks = c(-.2,.2))
```

<center>
```{r pcaSHOW, echo=FALSE}
plot_grid(p2,p3,p4,ncol = 3,nrow = 1)
```
</center>

### Fig 1 (complete)

Finally, we have all the pieces and can put together Fig. 1.

We need a few helper variables for the placing and scaling of the annotations:

```{r fig1Anno, results='hide', message=FALSE, warning=FALSE}
# scaling (x,y)
hSC <- c(.07,.07)
# hamlet placing (within group)
yf <- c(0.07,0.035,0.00);xf<-c(0.07,0,0.07)
# hamlet group placing
yS <- c(.485,.577,.805);xS<-c(.09,.417,.156)
```

Using `ggdraw()` form the `cowplot` package, the figure is put together.

```{r fig1, results='hide', message=FALSE, warning=FALSE}
F1 <- ggdraw()+
  # include subplots
  draw_plot(plot = p1,0,.38,.81,.64)+
  draw_plot(plot = p5,0,.01,.24,.35)+
  draw_plot(plot = p2,.25,.01,.23,.33)+
  draw_plot(plot = p3,.5,.01,.23,.33)+
  draw_plot(plot = p4,.75,.01,.23,.33)+
  # annotations (legend)
  draw_grob(legGrob, .375,.395,.65,.65)+
  # annotations (compas)
  draw_grob(compGrob, 0.605, 0.482, 0.06, 0.06)+
  # annotations (hamlets)
  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(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(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(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(uniGrob, xf[3]+xS[3], yf[3]+yS[3], hSC[1], hSC[2])+ # Belize
  # annotations (flags)
  draw_grob(panGrob, 0.3, 0.485, 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.285+.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)+
  # annotations (subfigure labels)  
  draw_plot_label(x = c(.005,.005,.25),
                  y = c(.995,.42,.42),
                  label =letters[1:3])
```
It is then exported using `ggsave()`.
```{r exportF1, eval=FALSE}
ggsave(plot = F1,filename = '../output/F1.pdf',width = 183,height = 145,units = 'mm',device = cairo_pdf)
  
```

<center>
```{r f1SHOW, echo=FALSE, fig.height=145*0.03937008, fig.width=183*0.03937008}
F1
```
</center>
---
 No newline at end of file
+49 −49
Original line number Diff line number Diff line
@@ -55,7 +55,6 @@ sampSize <- data %>% group_by(spec,loc) %>% summarise(n=length(id)) %>% filter(s
        x=c(rep(-94.7,3)+xpos, # belize
            rep(-78,3)+xpos, # panama
            rep(-98.75,3)+xpos)) # honduras
ccc<-rgb(0,.4,.8)

uS <- readShapeSpatial('../../0_data/0_resources/F1_shps/F1.npu_distribution.shp');
uS_f <- fortify(uS,region="DN");
@@ -69,6 +68,7 @@ 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") + 
@@ -84,6 +84,31 @@ p1 <- ggplot()+coord_map(projection = 'mercator',xlim=xlimW,ylim=ylimW)+
  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))
@@ -92,6 +117,22 @@ 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=1.1,shape=21)+
  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),
        plot.margin = unit(rep(-10,4),'pt'))+#coord_equal()+
  scale_x_continuous(name=xlabHON,breaks = c(-.2,.5))+
  scale_y_continuous(name=ylabHON,breaks = c(-.4,.2))

dataBEL <- cbind(dataAll %>% filter(sample=='sample',loc=='bel') %>% select(id,spec),
                 belPCA$scores); names(dataBEL)[3:12]<- paste('PC',1:10,sep='')

@@ -99,7 +140,7 @@ 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),'%)')

p2 <- ggplot(dataBEL,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,shape=21)+
p3 <- ggplot(dataBEL,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,shape=21)+
  scale_color_manual(values=clr,guide=F)+
  scale_fill_manual(values=fll,guide=F)+theme_mapK+
  theme(legend.position='bottom',
@@ -108,23 +149,8 @@ p2 <- ggplot(dataBEL,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,sh
        axis.title.x = element_text(vjust = 6),
        plot.margin = unit(rep(-10,4),'pt'))+
  scale_x_continuous(name=xlabBEL,breaks = c(-.15,.25))+
  scale_y_continuous(name=ylabBEL,breaks = c(-.3,.2));p2
  scale_y_continuous(name=ylabBEL,breaks = c(-.3,.2))

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)')
p3 <- ggplot(dataHON,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,shape=21)+
  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),
        plot.margin = unit(rep(-10,4),'pt'))+#coord_equal()+
  scale_x_continuous(name=xlabHON,breaks = c(-.2,.5))+
  scale_y_continuous(name=ylabHON,breaks = c(-.4,.2));p3

dataPAN <- cbind(dataAll %>% filter(sample=='sample',loc=='boc',spec!='gum')%>% select(id,spec),
                 panPCA$scores); names(dataPAN)[3:12]<- paste('PC',1:10,sep='')
@@ -140,39 +166,13 @@ p4 <- ggplot(dataPAN,aes(x=PC1,y=PC2,col=spec,fill=spec))+geom_point(size=1.1,sh
        axis.title.x = element_text(vjust = 6),
        plot.margin = unit(rep(-10,4),'pt'))+
  scale_x_continuous(name=xlabPAN,breaks = c(-.15,.2))+
  scale_y_continuous(name=ylabPAN,breaks = c(-.2,.2));p4


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));
  scale_y_continuous(name=ylabPAN,breaks = c(-.2,.2))

hSC <- c(.07,.07)
yf <- c(0.07,0.035,0.00);xf<-c(0.07,0,0.07)
yS <- c(.485,.577,.805);xS<-c(.09,.417,.156)

ggdraw()+
F1 <- ggdraw()+
  draw_plot(plot = p1,0,.38,.81,.64)+
  draw_plot(plot = p5,0,.01,.24,.35)+
  draw_plot(plot = p2,.25,.01,.23,.33)+
@@ -197,12 +197,12 @@ ggdraw()+
  draw_grob(belGrob, 0.2, 0.67, 0.042, 0.042)+
  draw_grob(honGrob, 0.28, 0.61, 0.042, 0.042)+
  
  draw_grob(belGrob, 0.285+.09, 0.37, 0.042, 0.042)+
  draw_grob(honGrob, 0.535+.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(panGrob, 0.785+.09, 0.37, 0.042, 0.042)+
  draw_plot_label(x = c(.005,.005,.25),
                  y = c(.995,.42,.42),
                  label =letters[1:3])

ggsave('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)
  
 No newline at end of file
+1 −1
Original line number Diff line number Diff line
@@ -3,7 +3,7 @@ bookdown::gitbook:
  split_by: chapter
  config:
    toc:
      collapse: subsection
      collapse: subsubsection
      before: |
        <li><a href="./">Script repository</a></li>
      after: |
+390 −1

File changed.

Preview size limit exceeded, changes collapsed.

+5.44 KiB

File added.

No diff preview for this file type.

Loading