Commit 489dd18f authored by Kosmas Hench's avatar Kosmas Hench

revision_R2

parent a5aeb3c5
......@@ -2,11 +2,11 @@ create_K_plot <- function(searchLG,gfffile,xr,searchgene,secondary_genes,searchs
source('../../0_data/0_scripts/F3.functions.R');
source('../../0_data/0_scripts/F3.getFSTs.R')
source('../../0_data/0_scripts/F3.getDXY.R')
source('../../0_data/0_scripts/S15.getGxP.R')
source('../../0_data/0_scripts/S15.getPI.R')
source('../../0_data/0_scripts/S15.getPIpw.R')
source('../../0_data/0_scripts/S15.getIHH12.R')
source('../../0_data/0_scripts/S15.getXPEHH.R')
source('../../0_data/0_scripts/S17.getGxP.R')
source('../../0_data/0_scripts/S17.getPI.R')
source('../../0_data/0_scripts/S17.getPIpw.R')
source('../../0_data/0_scripts/S17.getIHH12.R')
source('../../0_data/0_scripts/S17.getXPEHH.R')
highclr <- '#3bb33b'
theme_set(theme_minimal(base_size = 6))
......
#PBS -l elapstim_req=12:00:00
#PBS -l memsz_job=80gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N pcabel
#PBS -q cllong
#PBS -o 2.2.4.pca_bel.stdout
#PBS -e 2.2.4.pca_bel.stderr
cd $WORK/2_output/08_popGen/04_pca
mkdir -p bel_ld
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_bel.pop \
--thin 25000 \
--recode --stdout | sed 's=|=/=g' > bel_ld/pca.vcf
cd bel_ld && Rscript --vanilla ../pca_run.R pca.vcf bel $PWD
rm tmp.pcadapt && rm .Rhistory && rm pca.vcf
cd ..
echo "--done--"
#PBS -l elapstim_req=12:00:00
#PBS -l memsz_job=80gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N pcabel
#PBS -q cllong
#PBS -o 2.2.4.pca_hon.stdout
#PBS -e 2.2.4.pca_hon.stderr
cd $WORK/2_output/08_popGen/04_pca
mkdir -p hon_ld
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_hon.pop \
--thin 25000 \
--recode --stdout | sed 's=|=/=g' > hon_ld/pca.vcf
cd hon_ld && Rscript --vanilla ../pca_run.R pca.vcf hon $PWD
rm tmp.pcadapt && rm .Rhistory && rm pca.vcf
cd ..
echo "--done--"
#PBS -l elapstim_req=12:00:00
#PBS -l memsz_job=80gb
#PBS -b 1
#PBS -l cpunum_job=1
#PBS -N pcabel
#PBS -q cllong
#PBS -o 2.2.4.pca_pan.stdout
#PBS -e 2.2.4.pca_pan.stderr
cd $WORK/2_output/08_popGen/04_pca
mkdir -p pan_ld
vcftools --gzvcf $WORK/2_output/07_phased_variants/6_phased_mac2.vcf.gz \
--keep $WORK/0_data/0_resources/vcfpops/vcftools_boc.pop \
--thin 25000 \
--recode --stdout | sed 's=|=/=g' > pan_ld/pca.vcf
cd pan_ld && Rscript --vanilla ../pca_run.R pca.vcf pan $PWD
rm tmp.pcadapt && rm .Rhistory && rm pca.vcf
cd ..
echo "--done--"
......@@ -209,3 +209,4 @@ F1 <- ggdraw()+
ggsave(plot = F1,filename = '../output/F1.pdf',width = 183,height = 145,units = 'mm',device = cairo_pdf)
#ggsave(plot = F1,filename = '../output/F1.png',width = 183,height = 145,units = 'mm',dpi = 200)
#ggsave(plot = F1,filename = '../output/F1.eps',width = 183,height = 145,units = 'mm',device = cairo_ps)
\ No newline at end of file
......@@ -161,4 +161,6 @@ F2 <- ggdraw(p1)+
draw_grob(honGrob, 0.954, boxes$y[1]+.78*yd, .045, .045)+
draw_grob(panGrob, 0.954, boxes$y[3]+.78*yd, .045, .045)
ggsave(plot = F2,filename = '../output/F2.png',width = 183,height = 183,dpi = 200,units = 'mm')
#ggsave(plot = F2,filename = '../output/F2.png',width = 183,height = 183,dpi = 200,units = 'mm')
#ggsave(plot = F2,filename = '../output/F2.eps',width = 183,height = 183,units = 'mm',device = cairo_ps)
ggsave(plot = F2,filename = '../output/F2.pdf',width = 183,height = 183,units = 'mm',device = cairo_pdf)
\ No newline at end of file
......@@ -42,7 +42,6 @@ F3 <- plot_grid(NULL,NULL,NULL,NULL,
'','','',''),label_size = 10)+
draw_grob(legGrob, 0.1, 0, .8, 0.04)
ggsave(plot = F3,filename = '../output/F3.pdf',
width = 183,height = 125,units = 'mm',device = cairo_pdf)
#ggsave(plot = F3,filename = '../output/F3.png',
# width = 183,height = 125,dpi = 200,units = 'mm')
ggsave(plot = F3,filename = '../output/F3.pdf', width = 183,height = 125,units = 'mm',device = cairo_pdf)
#ggsave(plot = F3,filename = '../output/F3.png', width = 183,height = 125,dpi = 200,units = 'mm')
#ggsave(plot = F3,filename = '../output/F3.eps', width = 183,height = 125,units = 'mm',device = cairo_ps)
\ No newline at end of file
......@@ -72,4 +72,6 @@ F4 <- ggdraw(pGr1)+
y=c(.99,.81,.81,.57,.57),
size = 14)
ggsave(plot = F4,filename = '../output/F4.png',width = 183,height = 235,units = 'mm',dpi = 150)
#ggsave(plot = F4,filename = '../output/F4.png',width = 183,height = 235,units = 'mm',dpi = 150)
#ggsave(plot = F4,filename = '../output/F4.eps',width = 183,height = 235,units = 'mm',device = cairo_ps)
ggsave(plot = F4,filename = '../output/F4.pdf',width = 183,height = 235,units = 'mm',device = cairo_pdf)
\ No newline at end of file
#!/usr/bin/env Rscript
library(grid)
library(gridSVG)
library(grImport2)
library(grConvert)
library(tidyverse)
library(rtracklayer)
library(hrbrthemes)
library(cowplot)
require(rtracklayer)
source('../../0_data/0_scripts/F3.functions.R');
source('../../0_data/0_scripts/S15.plot_fun.R')
p1 <- create_K_plot(searchLG = "LG09",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(17800000,18000000),
searchgene = c("SOX10_1"),secondary_genes = c("RNASEH2A"),searchsnp = c(17871737,17872597,17873443),
muskID = 'A')
p2 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(20080000,20500000),
searchgene = c("Casz1_2"),secondary_genes = c(),searchsnp = c(20316944,20317120,20323661,20323670,20333895,20347263),
muskID = 'B')
p3 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22100000,22350000),
searchgene = c("hoxc13a"),secondary_genes = c('hoxc10a',"hoxc11a","hoxc12a","calcoco1_1"),
searchsnp = c(),
muskID = 'C')
p4 <- create_K_plot(searchLG = "LG17",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22500000,22670000),
searchgene = c('LWS',"SWS2abeta","SWS2aalpha","SWS2b"),searchsnp = c(22553970,22561903,22566254),
secondary_genes = c("Hcfc1","HCFC1_2","HCFC1_1","GNL3L","TFE3_0","MDFIC2_1","CXXC1_3","CXXC1_1",'Mbd1','CCDC120'),
muskID = 'D')
legGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend-pw-single-cairo.svg"))))
S15 <- plot_grid(NULL,NULL,NULL,NULL,
p1,NULL,p2,NULL,
NULL,NULL,NULL,NULL,
p3,NULL,p4,NULL,
NULL,NULL,NULL,NULL,
ncol=4,rel_heights = c(.03,1,.03,1,.1),rel_widths = c(1,.025,1,.04),
labels = c('a','','b','',
'','','','',
'c','','d','',
'','','','',
'','','',''),label_size = 10)+
draw_grob(legGrob, 0.05, 0, .9, 0.04)
ggsave(plot = S15,filename = '../output/S15.pdf',width = 183,height = 210,units = 'mm',device = cairo_pdf)
#ggsave(plot = S15,filename = '../output/S15.png',width = 183,height = 210,units = 'mm',dpi = 200)
source("../../0_data/0_scripts/F1.functions.R")
dataAll <- read.csv('../../0_data/0_resources/F1.sample.txt',sep='\t') %>%
mutate(loc=substrRight(as.character(id),3))
cFILL <- rgb(.7,.7,.7)
ccc<-rgb(0,.4,.8)
clr<-c('#000000','#d45500','#000000')
fll<-c('#000000','#d45500','#ffffff')
eVclr <- c(rep('darkred',2),rep(cFILL,8))
belPCA <- readRDS('../../2_output/08_popGen/04_pca/bel_ld/belpca.Rds')
honPCA <- readRDS('../../2_output/08_popGen/04_pca/hon_ld/honpca.Rds')
panPCA <- readRDS('../../2_output/08_popGen/04_pca/pan_ld/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)+
ggtitle('Honduras')+
scale_fill_manual(values=fll,guide=F)+theme_mapK+
theme(legend.position='bottom',plot.title = element_text(size = 11,hjust = .45),
panel.border = element_rect(color=rgb(.9,.9,.9),fill=rgb(1,1,1,0)),
plot.margin = unit(c(3,rep(5,3)),'pt'))+#coord_equal()+
scale_x_continuous(name=xlabHON,breaks = c(-.3,.2))+
scale_y_continuous(name=ylabHON,breaks = c(-.4,.4))
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),'%)')
p1 <- 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)+
ggtitle('Belize')+
scale_fill_manual(values=fll,guide=F)+theme_mapK+
theme(legend.position='bottom',plot.title = element_text(size = 11,hjust = .45),
panel.border = element_rect(color=rgb(.9,.9,.9),fill=rgb(1,1,1,0)),
plot.margin = unit(c(3,rep(5,3)),'pt'))+
scale_x_continuous(name=xlabBEL,breaks = c(-.15,.25))+
scale_y_continuous(name=ylabBEL,breaks = c(-.35,.25))
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),'%)')
p3 <- 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)+
ggtitle('Panama')+
scale_fill_manual(values=fll,guide=F)+theme_mapK+
theme(legend.position='bottom',plot.title = element_text(size = 11,hjust = .45),
panel.border = element_rect(color=rgb(.9,.9,.9),fill=rgb(1,1,1,0)),
plot.margin = unit(c(3,rep(5,3)),'pt'))+
scale_x_continuous(name=xlabPAN,breaks = c(-.15,.2))+
scale_y_continuous(name=ylabPAN,breaks = c(-.2,.2))
S15 <- cowplot::plot_grid(p2, p1, p3, nrow = 1, labels = letters[1:3])
ggsave(plot = S15,filename = '../output/S15.pdf',width = 183,height = 63,units = 'mm',device = cairo_pdf)
#ggsave(plot = S15,filename = '../output/S15.png',width = 183,height = 63,units = 'mm',dpi = 200)
#ggsave(plot = S15,filename = '../output/S15.eps',width = 183,height = 63,units = 'mm',device = cairo_ps)
\ No newline at end of file
#!/usr/bin/env Rscript
library(tidyverse)
dat_dir <- '../../2_output/08_popGen/06_GxP/smoothed/50kb/'
files <- c("nig-uni-gemma.lm.50k.5k.smooth.gz", "nig-uni.lmm.50k.5k.smooth.gz",
"pue-nig-gemma.lm.50k.5k.smooth.gz", "pue-nig.lmm.50k.5k.smooth.gz",
"pue-uni-gemma.lm.50k.5k.smooth.gz", "pue-uni.lmm.50k.5k.smooth.gz")
runs <- files %>% str_remove(.,'-gemma') %>% str_remove(.,'.50k.5k.smooth.gz') %>% str_replace(.,'.lm','_lm')
karyo <- read_delim('../../0_data/0_resources/F2.karyo.txt', delim = '\t') %>%
mutate(GSTART=lag(cumsum(END),n = 1,default = 0),
GEND=GSTART+END,GROUP=rep(letters[1:2],12)) %>%
select(CHROM,GSTART,GEND,GROUP)
import <- function(file,run){
read_delim(file = gzfile(file), delim = "\t") %>%
dplyr::left_join(karyo, by = "CHROM") %>%
dplyr::mutate(POS = (BIN_START + BIN_END)/2,
GPOS = GSTART + POS) %>%
mutate(RUN = run)
}
data <- purrr::map2(.x = str_c(dat_dir,'/',files),runs,import) %>%
bind_rows() %>% separate(RUN,into = c('COMP','TYPE'), sep = '_')
clr <- c('#fb8620','#d93327','#1b519c')
yl <- expression(-log[10](italic(p)~value))
S16 <- ggplot(data,aes(x=GPOS,y=avgp_wald,col=COMP))+
geom_rect(inherit.aes = F,data=karyo,aes(xmin=GSTART,xmax=GEND,
ymin=-Inf,ymax=Inf,fill=GROUP))+
geom_point(size=.01)+
facet_grid(COMP+TYPE~.)+
scale_x_continuous(name = "", expand = c(0, 0), breaks = (karyo$GSTART + karyo$GEND)/2,
labels = 1:24, position = "top")+
ylab(yl)+
scale_fill_manual(values = c(NA,"#E6E6E6"), guide = F)+
scale_color_manual(name='comparison',values=clr[c(3,1,2)])+
theme_bw(base_size = 10, base_family = "Helvetica") +
theme(plot.background = element_blank(),
panel.background = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
axis.title.x = element_blank(),
axis.line = element_line(),
legend.position = 'none',
strip.background = element_blank())
#ggsave(plot = S16, filename = '../output/S16.pdf',width = 183,height = 60,units = 'mm',device = cairo_pdf)
ggsave(plot = S16, filename = '../output/S16.png',width = 183,height = 90,units = 'mm',dpi = 300)
#!/usr/bin/env Rscript
library(grid)
library(gridSVG)
library(grImport2)
library(grConvert)
library(tidyverse)
library(rtracklayer)
library(hrbrthemes)
library(cowplot)
require(rtracklayer)
source('../../0_data/0_scripts/F3.functions.R');
source('../../0_data/0_scripts/S17.plot_fun.R')
p1 <- create_K_plot(searchLG = "LG09",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(17800000,18000000),
searchgene = c("SOX10_1"),secondary_genes = c("RNASEH2A"),searchsnp = c(17871737,17872597,17873443),
muskID = 'A')
p2 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(20080000,20500000),
searchgene = c("Casz1_2"),secondary_genes = c(),searchsnp = c(20316944,20317120,20323661,20323670,20333895,20347263),
muskID = 'B')
p3 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22100000,22350000),
searchgene = c("hoxc13a"),secondary_genes = c('hoxc10a',"hoxc11a","hoxc12a","calcoco1_1"),
searchsnp = c(),
muskID = 'C')
p4 <- create_K_plot(searchLG = "LG17",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22500000,22670000),
searchgene = c('LWS',"SWS2abeta","SWS2aalpha","SWS2b"),searchsnp = c(22553970,22561903,22566254),
secondary_genes = c("Hcfc1","HCFC1_2","HCFC1_1","GNL3L","TFE3_0","MDFIC2_1","CXXC1_3","CXXC1_1",'Mbd1','CCDC120'),
muskID = 'D')
legGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend-pw-single-cairo.svg"))))
S17 <- plot_grid(NULL,NULL,NULL,NULL,
p1,NULL,p2,NULL,
NULL,NULL,NULL,NULL,
p3,NULL,p4,NULL,
NULL,NULL,NULL,NULL,
ncol=4,rel_heights = c(.03,1,.03,1,.1),rel_widths = c(1,.025,1,.04),
labels = c('a','','b','',
'','','','',
'c','','d','',
'','','','',
'','','',''),label_size = 10)+
draw_grob(legGrob, 0.05, 0, .9, 0.04)
ggsave(plot = S17,filename = '../output/S17.pdf',width = 183,height = 210,units = 'mm',device = cairo_pdf)
#ggsave(plot = S17,filename = '../output/S17.png',width = 183,height = 210,units = 'mm',dpi = 200)
This diff is collapsed.
---
output: html_document
editor_options:
chunk_output_type: console
---
# Supplementary Figure 16
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = './F_scripts')
```
## Summary
This is the accessory documentation of Supplementary Figure 16.
The Figure can be recreated by running the **R** script `S16.R`:
```sh
cd $WORK/3_figures/F_scripts
Rscript --vanilla S16.R
rm Rplots.pdf
```
## Details of `S16.R`
In S16.R script is the `F3.R` script run including additional data sets.
Is an executable R script that depends on the **tidyverse** package for data managing and plotting.
```{r import, results='hide', message=FALSE, warning=FALSE}
library(tidyverse)
```
### Setup
First, we define the directory storing the GxP results and create an inventory of data files to load.
We also create labels for the individual data sets.
```{r setup_data, results='hide', message=FALSE, warning=FALSE}
dat_dir <- '../../2_output/08_popGen/06_GxP/smoothed/50kb/'
files <- c("nig-uni-gemma.lm.50k.5k.smooth.gz", "nig-uni.lmm.50k.5k.smooth.gz",
"pue-nig-gemma.lm.50k.5k.smooth.gz", "pue-nig.lmm.50k.5k.smooth.gz",
"pue-uni-gemma.lm.50k.5k.smooth.gz", "pue-uni.lmm.50k.5k.smooth.gz")
runs <- files %>%
str_remove(.,'-gemma') %>%
str_remove(.,'.50k.5k.smooth.gz') %>%
str_replace(.,'.lm','_lm')
```
We load information about hamlet reference genome and create a import function to load the data.
```{r karyo, results='hide', message=FALSE, warning=FALSE}
karyo <- read_delim('../../0_data/0_resources/F2.karyo.txt', delim = '\t') %>%
mutate(GSTART=lag(cumsum(END),n = 1,default = 0),
GEND=GSTART+END,GROUP=rep(letters[1:2],12)) %>%
select(CHROM,GSTART,GEND,GROUP)
import <- function(file,run){
read_delim(file = gzfile(file), delim = "\t") %>%
dplyr::left_join(karyo, by = "CHROM") %>%
dplyr::mutate(POS = (BIN_START + BIN_END)/2,
GPOS = GSTART + POS) %>%
mutate(RUN = run)
}
```
We use the `map2()` function from the **purrr** package to load the data sets.
The we define some colors an the y-axis label.
```{r lead_data, results='hide', message=FALSE, warning=FALSE}
data <- purrr::map2(.x = str_c(dat_dir,'/',files),runs,import) %>%
bind_rows() %>% separate(RUN,into = c('COMP','TYPE'), sep = '_')
clr <- c('#fb8620','#d93327','#1b519c')
yl <- expression(-log[10](italic(p)~value))
```
Finally we plot the data.
```{r plot, results='hide', message=FALSE, warning=FALSE}
S16 <- ggplot(data,aes(x=GPOS,y=avgp_wald,col=COMP))+
geom_rect(inherit.aes = F,data=karyo,aes(xmin=GSTART,xmax=GEND,
ymin=-Inf,ymax=Inf,fill=GROUP))+
geom_point(size=.01)+
facet_grid(COMP+TYPE~.)+
scale_x_continuous(name = "", expand = c(0, 0), breaks = (karyo$GSTART + karyo$GEND)/2,
labels = 1:24, position = "top")+
ylab(yl)+
scale_fill_manual(values = c(NA,"#E6E6E6"), guide = F)+
scale_color_manual(name='comparison',values=clr[c(3,1,2)])+
theme_bw(base_size = 10, base_family = "Helvetica") +
theme(plot.background = element_blank(),
panel.background = element_blank(),
panel.grid = element_blank(),
panel.border = element_blank(),
axis.title.x = element_blank(),
axis.line = element_line(),
legend.position = 'none',
strip.background = element_blank())
```
The Figure is then exported using `ggsave()`:
```{r, eval=FALSE}
#ggsave(plot = S16, filename = '../output/S16.pdf',width = 183,height = 60,units = 'mm',device = cairo_pdf)
ggsave(plot = S16, filename = '../output/S16.png',width = 183,height = 90,units = 'mm',dpi = 300)
```
<center>
```{r s16SHOW, echo=FALSE, fig.width=183*0.03937008, fig.height=90*0.03937008, warning=FALSE, message=FALSE}
S16
```
</center>
---
\ No newline at end of file
---
output: html_document
editor_options:
chunk_output_type: console
---
# Supplementary Figure 17
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = './F_scripts')
```
## Summary
This is the accessory documentation of Supplementary Figure 17.
The Figure can be recreated by running the **R** script `S17.R`:
```sh
cd $WORK/3_figures/F_scripts
Rscript --vanilla S17.R
rm Rplots.pdf
```
## Details of `S17.R`
In S17.R script is the `F3.R` script run including additional data sets.
Is an executable R script that depends on a variety of image manipulation and data managing and genomic data packages.
It depends on the same helper script as `F3.R` (`F3.functions.R`). Additionally it depends on `S17.plot_fun.R` (modified `F3.plot_fun.R`) which itself depends on `F3.getDXY.R`, `F3.getFSTs.R`,`S17.getGxP.R`,`S17.getIHH12.R`, `S17.getPIpw.R`, `S17.getPI.R` and `S17.getXPEHH.R` (all located under `$WORK/0_data/0_scripts`).
For a detailed documentation please refer to [F3.R](figure-3.html). The only differences between the two scripts are located in `S17.plot_fun.R` (compared to `F3.plot_fun.R`).
### Additions in `S17.plot_fun.R`
In the following additions to `S17.plot_fun.R` are shown which are missing in the original.
In the beginning, four new helper scripts are loaded which include the data import for the additional parameter which are plotted in Supplementary Figure 15 ($\pi$, pairwise $\pi$, iHH~12~ and xpEHH).
```{r, eval=FALSE}
#....#
source('../../0_data/0_scripts/S17.getGxP.R')
source('../../0_data/0_scripts/S17.getPI.R')
source('../../0_data/0_scripts/S17.getPIpw.R')
source('../../0_data/0_scripts/S17.getIHH12.R')
source('../../0_data/0_scripts/S17.getXPEHH.R')
#....#
```
The imported data are stored in `data_*` vectors.
```{r, eval=FALSE}
#....#
# get pfst values
pfst_list <- getGxP(searchLG,xr)
# data_pfst_pw<-pfst_list$data_pfst_pw;
data_pfst<-pfst_list$data_pfst
# get pi values
piPW_list <- getPIpw(searchLG,xr)
data_piPW<-piPW_list$data_piPW
#....#
# get xpEHH values
xpEHH_list <- getXPEHH(searchLG,xr)
data_xpEHH<-xpEHH_list$data_xpEHH
#....#
```
For each new parameter, an additional plot is created.
```{r, eval=FALSE}
#....#
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)+
#....#
p15 <- ggplot()+coord_cartesian(xlim=xr)+
geom_line(data=(data_piPW %>% filter(POS > xr[1],POS<xr[2]))
,aes(x=POS,y=PI,col=run),lwd=LW)+
#....#
p18 <- ggplot()+coord_cartesian(xlim=xr)+
geom_line(data=(data_xpEHH %>% filter(POS > xr[1],POS<xr[2]))
,aes(x=POS,y=AVGxxEHH,col=run),lwd=LW)+
#....#
#....#
```
The plot composition returned by the function includes the additional new plots:
```{r, eval=FALSE}
#....#
p2 <- plot_grid(p11,p12,p13,p14,p15,p16,p17,p18,
ncol = 1,align = 'v',axis = 'r',rel_heights = c(.85,rep(1,6)))
return(p2)
```
Comparing to the original script (`F3.plot_fun.R`):
<div class="notthescript">
```{r oldScript2, eval=FALSE,}
#....#
p2 <- plot_grid(p11,p12,p13,p14,
ncol = 1,align = 'v',axis = 'r',rel_heights = c(1.3,rep(1,3)))
return(p2)
```
</div>
```{r head, include=FALSE}
library(grid)
library(gridSVG)
library(grImport2)
library(grConvert)
library(tidyverse)
library(rtracklayer)
library(hrbrthemes)
library(cowplot)
require(rtracklayer)
source('../../0_data/0_scripts/F3.functions.R');
source('../../0_data/0_scripts/S17.plot_fun.R')
p1 <- create_K_plot(searchLG = "LG09",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(17800000,18000000),
searchgene = c("SOX10_1"),secondary_genes = c("RNASEH2A"),searchsnp = c(17871737,17872597,17873443),
muskID = 'A')
p2 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(20080000,20500000),
searchgene = c("Casz1_2"),secondary_genes = c(),searchsnp = c(20316944,20317120,20323661,20323670,20333895,20347263),
muskID = 'B')
p3 <- create_K_plot(searchLG = "LG12",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22100000,22350000),
searchgene = c("hoxc13a"),secondary_genes = c('hoxc10a',"hoxc11a","hoxc12a","calcoco1_1"),
searchsnp = c(),
muskID = 'C')
p4 <- create_K_plot(searchLG = "LG17",gfffile = '../../1-output/09_gff_from_IKMB/HP.annotation.named.gff',xr = c(22500000,22670000),
searchgene = c('LWS',"SWS2abeta","SWS2aalpha","SWS2b"),searchsnp = c(22553970,22561903,22566254),
secondary_genes = c("Hcfc1","HCFC1_2","HCFC1_1","GNL3L","TFE3_0","MDFIC2_1","CXXC1_3","CXXC1_1",'Mbd1','CCDC120'),
muskID = 'D')
legGrob <- gTree(children=gList(pictureGrob(readPicture("../../0_data/0_img/legend-pw-cairo.svg"))))
S17 <- plot_grid(NULL,NULL,NULL,NULL,
p1,NULL,p2,NULL,
NULL,NULL,NULL,NULL,
p3,NULL,p4,NULL,
NULL,NULL,NULL,NULL,
ncol=4,rel_heights = c(.03,1,.03,1,.1),rel_widths = c(1,.025,1,.02),
labels = c('a','','b','',
'','','','',
'c','','d','',
'','','','',
'','','',''),label_size = 10)+
draw_grob(legGrob, 0.1, 0, .8, 0.04)
```
<center>
```{r s09SHOW, echo=FALSE, fig.width=183*0.03937008, fig.height=210*0.03937008, warning=FALSE, message=FALSE}
S17
```
</center>
---
\ No newline at end of file
......@@ -62,9 +62,9 @@ edgsB <- rbind(edgR(ptsB1,edgB1),edgR(ptsB2,edgB2),edgR(ptsB3,edgB3))
plotList <- data.frame(x=c(0,.2,.66,.84,1,1.9,2.1,1.5,2.5,
2.64,2.64,2.64,2.64,0,3.5,3.5,.45),
y=c(rep(1,9),-1,-.8,-.6,-.4,-.05,0,.75,1),
label=c('S01','F1','F2','S07','S12','F3','S09, S08',
'S03','S15','S14','S11',
'F4','S14','S13','S06','S15','S02'),
label=c('S01','F1, S15','F2','S07','S12','F3','S09, S08',
'S03, S16','S15','S14','S11',
'F4','S14','S13','S06','S17','S02'),
status=c(rep(steps[5],14),rep(steps[6],2),steps[5]))
```
<center>
......@@ -117,7 +117,7 @@ p2 <- ggplot()+
arrow = arrow(length = unit(4,'pt'),type = 'closed'))+
geom_text(data=edgsB,aes(x=(x+xend)*.5,y=(y+yend)*.5,label=label,angle=angle),
vjust=-1)+
#geom_text(data=ptsB%>% filter(!(status==steps[5] & nr == 20)),aes(x=x,y=y,label=nr),color='white')+
# geom_text(data=ptsB%>% filter(!(status==steps[5] & nr == 20)),aes(x=x,y=y,label=nr),color='white')+
geom_text(data=plotList,
aes(x=x,y=y,label=label),color=clr2[4])+
scale_x_continuous(expand=c(.03,.03))+
......
......@@ -2,7 +2,7 @@ book_filename: "Script repository Hench et al. 2018"
chapter_name: ""
repo: https://github.com/k-hench/bookdown
output_dir: ../docs
rmd_files: ["index.Rmd","Workflow.Rmd","F1.Rmd","F2.Rmd","F3.Rmd","F4.Rmd","S01.Rmd","S02.Rmd","S03.Rmd","S05.Rmd","S06.Rmd","S07.Rmd","S08.Rmd","S09.Rmd","S10.Rmd","S11.Rmd","S12.Rmd","S13.Rmd","S14.Rmd","S15.Rmd"]
rmd_files: ["index.Rmd","Workflow.Rmd","F1.Rmd","F2.Rmd","F3.Rmd","F4.Rmd","S01.Rmd","S02.Rmd","S03.Rmd","S05.Rmd","S06.Rmd","S07.Rmd","S08.Rmd","S09.Rmd","S10.Rmd","S11.Rmd","S12.Rmd","S13.Rmd","S14.Rmd","S15.Rmd","S16.Rmd","S17.Rmd"]
clean: [packages.bib, bookdown.bbl]
new_session: yes
delete_merged_file: True
......@@ -34,7 +34,8 @@ A more detailed documentation exists for all the figures of the manuscript:
as well as for all the supplementary figures:
[S01](supplementary-figure-01.html), [S02](supplementary-figure-02.html), [S03](supplementary-figure-03.html), [S05](supplementary-figure-05.html), [S06](supplementary-figure-06.html), [S07](supplementary-figure-07.html), [S08](supplementary-figure-08.html), [S09](supplementary-figure-09.html), [S10](supplementary-figure-10.html), [S11](supplementary-figure-11.html), [S12](supplementary-figure-12.html), [S13](supplementary-figure-13.html), [S14](supplementary-figure-14.html) & [S15](supplementary-figure-15.html)
[S01](supplementary-figure-01.html), [S02](supplementary-figure-02.html), [S03](supplementary-figure-03.html), [S05](supplementary-figure-05.html), [S06](supplementary-figure-06.html), [S07](supplementary-figure-07.html), [S08](supplementary-figure-08.html), [S09](supplementary-figure-09.html), [S10](supplementary-figure-10.html), [S11](supplementary-figure-11.html), [S12](supplementary-figure-12.html), [S13](supplementary-figure-13.html), [S14](supplementary-figure-14.html), [S15](supplementary-figure-15.html),
[S16](supplementary-figure-16.html) & [S17](supplementary-figure-17.html)
The only exception to this is the supplementary figure S04. This figure is a byproduct of the anchoring step during the assembly and was produced by the **Allmaps** software. Afterwards, **Inkscape** was used to adjust the coloration and labels of the linkage maps.
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.