F3.figure3.Rmd

---
title: "F3.figure3"
author: "Daniele Filiault"
date: "7/28/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(RColorBrewer)
library(grid)
library(gridExtra)
library(ggpubr)
library(cowplot)
library(ggh4x)
library(viridis)
```


## Introduction

Generates figure 3 of manuscript. 
panel A is cumulative beta distribution
panel B is home allele frequency by beta bin - proportions
panel C is home allele frequency by beta bin - numbers

## 01.  Get data
Data and analysis is from script 51

```{r get data}
# for panel A
load(file="./data/51.data/n.af.sum.Rdat")
load(file="./data/51.data/s.af.sum.Rdat")

# for panels B and C
load(file="./data/51.data/sbeta.sum.Rdat")
load(file="./data/51.data/nbeta.sum.Rdat")
```

  
## 02. make cumulative beta plots
```{r cumul beta data prep}
af.sum <- rbind(n.af.sum, s.af.sum)
rm(n.af.sum)
rm(s.af.sum)
af.loc <- do.call(rbind, strsplit(af.sum$exp,"_"))
af.sum$exp.site <- af.loc[,1]
af.sum$year <- af.loc[,2]
af.sum$ahome.bins <- factor(af.sum$ahome.bins, levels=rev(levels(af.sum$ahome.bins)))

## also need borders of fixed for vertical lines
af.fixed <- aggregate(index~exp, data=af.sum[af.sum$ahome==1,], FUN=min)
af.loc <- do.call(rbind, strsplit(af.fixed$exp,"_"))
af.fixed$exp.site <- af.loc[,1]
af.fixed$year <- af.loc[,2]

```

```{r cumul beta plot}

### plot weighed cumulative beta
panel.a <- ggplot(data=af.sum, aes(x=index,y=csum.wb.norm)) +
  geom_point(aes(color=ahome.bins)) +
  scale_color_brewer(palette = "Paired", name="allele\nfrequency\nbin") +
  #scale_color_viridis(direction=-1,name="home\nallele\nfrequency") +
  facet_nested(exp.site + year~.)  +
  #facet_wrap(~exp, ncol=2) +
  geom_hline(yintercept=0,linetype="dashed", color = "red") +
  xlab("number of SNPs") +
  ylab("cumulative weighted beta (normalized)") + 
  geom_vline(data=af.fixed, aes(xintercept=index), linetype="dashed", color="purple") +
  theme(plot.margin=unit(c(1,0.5,0.5,0.5), "lines"), legend.position="None")

```

## 03. make barplots (panels B and C)
```{r AF comp by beta North}

af.beta.abs.n <- as.list(1:4)
af.beta.prop.n <- as.list(1:4)
up.exp <- levels(nbeta.sum$exp)

for(up in 1:4){
  up.dat <- nbeta.sum[nbeta.sum$exp==up.exp[up],]
  up.dat$beta.bin <- cut(up.dat$home.beta,10)
  up.dat$nb.bin <- as.numeric(up.dat$beta.bin)
  #get counts
  up.tab <- with(up.dat, table(ahome.bins, nb.bin))
  up.df <- as.data.frame(up.tab)
  up.df$exp <- up.exp[up]
  af.beta.abs.n[[up]] <- up.df
  #get as proportions
  up.ptab <- prop.table(up.tab, margin=2)
  up.pdf <- as.data.frame(up.ptab)
  up.pdf$exp <- up.exp[up]
  af.beta.prop.n[[up]] <- up.pdf
}
af.beta.abs.n <- do.call(rbind, af.beta.abs.n)
af.beta.prop.n <- do.call(rbind, af.beta.prop.n)
```

```{r AF comp by beta South}

af.beta.abs.s <- as.list(1:4)
af.beta.prop.s <- as.list(1:4)
up.exp <- levels(sbeta.sum$exp)

for(up in 1:4){
  up.dat <- sbeta.sum[sbeta.sum$exp==up.exp[up],]
  up.dat$beta.bin <- cut(up.dat$home.beta,10)
  up.dat$nb.bin <- as.numeric(up.dat$beta.bin)
  #get counts
  up.tab <- with(up.dat, table(ahome.bins, nb.bin))
  up.df <- as.data.frame(up.tab)
  up.df$exp <- up.exp[up]
  af.beta.abs.s[[up]] <- up.df
  #get as proportions
  up.ptab <- prop.table(up.tab, margin=2)
  up.pdf <- as.data.frame(up.ptab)
  up.pdf$exp <- up.exp[up]
  af.beta.prop.s[[up]] <- up.pdf
}
af.beta.abs.s <- do.call(rbind, af.beta.abs.s)
af.beta.prop.s <- do.call(rbind, af.beta.prop.s)
```

```{r combine n and s data}
af.beta.abs <- rbind(af.beta.abs.n, af.beta.abs.s)
af.beta.prop <- rbind(af.beta.prop.n, af.beta.prop.s)
af.loc <- do.call(rbind, strsplit(af.beta.abs$exp,"_"))
af.beta.abs$exp.site <- af.loc[,1]
af.beta.abs$year <- af.loc[,2]
af.beta.prop$exp.site <- af.loc[,1]
af.beta.prop$year <- af.loc[,2]

af.beta.abs$ahome.bins <- factor(af.beta.abs$ahome.bins, levels=rev(levels(af.beta.abs$ahome.bins)))
af.beta.prop$ahome.bins <- factor(af.beta.prop$ahome.bins, levels=rev(levels(af.beta.prop$ahome.bins)))
```


```{r make proportions bin plot}

panel.b <- ggplot(data=af.beta.prop, aes(x=nb.bin,y=Freq,fill=ahome.bins)) +
  geom_bar(position="stack",stat="identity") + 
  coord_flip()+
  ylab("proportion of SNPs") +
  xlab("home beta bin") +
  scale_fill_brewer(palette = "Paired", name="allele\nfrequency\nbin", guide=guide_legend(reverse=T)) +
  #facet_nested(exp.site + year~.)  +
  facet_grid(exp ~.) +
  theme(strip.background=element_blank(), strip.text=element_blank(), plot.margin=unit(c(1,0.5,0.5,1), "lines"), legend.position = "none")
```

```{r make absolute numbers bin plot}

panel.c <- ggplot(data=af.beta.abs, aes(x=nb.bin,y=Freq,fill=ahome.bins)) +
  geom_bar(position="stack",stat="identity") + 
  coord_flip()+
  ylab("number of SNPs") +
  xlab("home beta bin") +
  scale_fill_brewer(palette = "Paired", name="allele\nfrequency\nbin", guide=guide_legend(reverse=T)) +
  #facet_grid(exp.site + year~.)
  facet_nested(exp.site + year~.)+
  #facet_grid(exp ~.) +
  theme(axis.text.y=element_blank(), plot.margin=unit(c(1,1,0.5,0), "lines")) 

```


## 04. combine panels to final figure

```{r generate final figure}

fig3 <- ggarrange(panel.a, panel.b, panel.c + rremove("ylab"),
                    ncol = 3,
                    labels=c("A","B","C"),
                    font.label = list(size = 10, color = "black", face = "bold", family = NULL, position = "top"),
                    widths=c(1.5,1.2,1.5))

ggsave(fig3, file="./manuscript.figures/Figure3.jpg", width=10, height=9)

```

## 05.  subsets of figure for talks

```{r subsets for talks}
#cumulative betas for NM_2011 and SU_2011
af.subset <- af.sum[af.sum$exp%in%c("NM_2011", "SU_2011"),]

csum.subset.plot <- ggplot(data=af.subset, aes(x=index,y=csum.wb.norm)) +
  geom_point(aes(color=ahome.bins)) +
  scale_color_brewer(palette = "Paired", name="allele\nfrequency\nbin", guide=guide_legend(reverse=T)) +
  facet_wrap(~exp, ncol=1) +
  geom_hline(yintercept=0,linetype="dashed", color = "red") +
  xlab("number of SNPs") +
  ylab("cumulative weighted beta (normalized)") + 
  geom_vline(data=af.fixed[af.fixed$exp%in%c("NM_2011", "SU_2011"),], aes(xintercept=index), linetype="dashed", color="purple") +
  theme(plot.margin=unit(c(1,0.5,0.5,0.5), "lines"))
ggsave(csum.subset.plot, file="./figures/99.talk.figures/subset.cumsum.jpg", width=6, height=5)


#proportions for NM_2011 and SU_2011
prop.subset.plot <- ggplot(data=af.beta.prop[af.beta.prop$exp%in%c("NM_2011", "SU_2011"),], aes(x=nb.bin,y=Freq,fill=ahome.bins)) +
  geom_bar(position="stack",stat="identity") + 
  coord_flip()+
  ylab("proportion of SNPs") +
  xlab("home beta bin") +
  scale_fill_brewer(palette = "Paired", name="allele\nfrequency\nbin", guide=guide_legend(reverse=T)) +
  #facet_nested(exp.site + year~.)  +
  facet_wrap(~exp, ncol=1) +
  theme(plot.margin=unit(c(1,0.5,0.5,1), "lines"))
ggsave(prop.subset.plot, file="./figures/99.talk.figures/subset.prop.jpg", width=4, height=5)

#absolute numbers for NM_2011 and SU_2011
absolute.subset.plot <- ggplot(data=af.beta.abs[af.beta.abs$exp%in%c("NM_2011", "SU_2011"), ], aes(x=nb.bin,y=Freq,fill=ahome.bins)) +
  geom_bar(position="stack",stat="identity") + 
  coord_flip()+
  ylab("number of SNPs") +
  xlab("home beta bin") +
  scale_fill_brewer(palette = "Paired", name="allele\nfrequency\nbin", guide=guide_legend(reverse=T)) +
  facet_wrap(~exp, ncol=1) +
  theme(plot.margin=unit(c(1,0.5,0.5,1), "lines")) 
ggsave(absolute.subset.plot, file="./figures/99.talk.figures/subset.absolute.jpg", width=4, height=5)


subset.homeAF.plot <- ggarrange(csum.subset.plot, prop.subset.plot, absolute.subset.plot+ rremove("ylab"), ncol=3, common.legend = TRUE, legend="right", widths=c(2,1,0.85))
ggsave(subset.homeAF.plot, file="./figures/99.talk.figures/subset.homeAF.plot.jpg", width=9, height=5)

```