In nalcala/PopGenBounds: Bounds on population-genetic statistics

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This vignette shows an example analysis of yellow toad microsatellite data from Pohele et al. 2021, used to produce figure 6 from Alcala et al., In prep.

library(PopGenBounds)
library(tidyverse)
library(patchwork)

We load the dataset

Toad_microsat_freq = PopGenBounds::Toad_microsat_freq

Plotting the allele frequency tables

ggfreqtable(Toad_microsat_freq[[1]])

Computing the statistics

For all subpopulations together (K=47)

Diff_toad = lapply(Toad_microsat_freq,Diff)

For pairs of subpopulations (K=2)

Diff_toad_pairs = lapply(Toad_microsat_freq,function(f){sapply(1:46, function(i){ sapply((i+1):47, function(j) Diff(f[c(i,j),])$value)})} )

Plotting the data

For K=47

toad.tib = tibble(M=unlist( sapply(Diff_toad,function(x){x[1,2]})),
                  FST=unlist( sapply(Diff_toad,function(x){x[2,2]})),
                  GpST=unlist( sapply(Diff_toad,function(x){x[3,2]})),
                  D=unlist( sapply(Diff_toad,function(x){x[4,2]})))

ggtoad = ggbounds(M=toad.tib$M,FST=toad.tib$FST,GpST=toad.tib$GpST,D=toad.tib$D,K=47)

ggtoad[[1]]+ ggtoad[[2]]+ ggtoad[[3]]

For K=2

toad.tib_pairs = tibble(M=as.numeric(sapply(Diff_toad_pairs, function(x) unlist(sapply(x,function(y){y[1,]})))),
                  FST=as.numeric(sapply(Diff_toad_pairs, function(x) unlist(sapply(x,function(y){y[2,]})))),
                  GpST=as.numeric(sapply(Diff_toad_pairs, function(x) unlist(sapply(x,function(y){y[3,]})))),
                  D=as.numeric(sapply(Diff_toad_pairs, function(x) unlist(sapply(x,function(y){y[4,]})))) )

# when M=1, no polymorphism, so we exclude the locus
toad.tib_pairs = toad.tib_pairs %>% filter(M<1)

ggtoad_pairs = ggbounds(M=toad.tib_pairs$M,FST=toad.tib_pairs$FST,GpST=toad.tib_pairs$GpST,D=toad.tib_pairs$D,K=2)

ggtoad_pairs[[1]]+ ggtoad_pairs[[2]]+ ggtoad_pairs[[3]]

#ggsave(filename = "Fig_toad.svg",((ggtoad_pairs[[1]]+ scale_color_viridis(limits=c(0,5000)) +guides(color="none")) + 
#  (ggtoad_pairs[[2]]+ scale_color_viridis(limits=c(0,5000)) +guides(color="none") ) + 
#  (ggtoad_pairs[[3]] + scale_color_viridis(limits=c(0,5000)))) / ((ggtoad[[1]]+ scale_color_viridis(limits=c(0,7)) +guides(color="none")) + 
#  (ggtoad[[2]]+ scale_color_viridis(limits=c(0,7)) +guides(color="none") ) + 
#  (ggtoad[[3]] + scale_color_viridis(limits=c(0,7)))),height = 2.3*2,width=3*3)