GppFst: GppFst: function to simulate posterior predictive...

In radamsRHA/GppFst: Posterior Predictive Simulations for generating theoretical distributions of Fst and Dxy under the neutral coalescent model for two species

Description

This function returns a vector of Fst values simulated from a posterior distribution

Usage

GppFst(posterior.samples, loci.per.step, SNP.per.locus, sample.vec.0,
  sample.vec.1, sequence.length.vec)

Arguments

posterior.samples	List of posterior samples with columns for single column for each parameter and a row for each step of the MCMC chain [see read.posterior]
loci.per.step	Number of genealogies to simulate for each joint parameter combination of the posterior MCMC chain
SNP.per.locus	Total number of fst values taken. This is arbitrary, and can be set to the locus length to get all simulated SNPs
sample.vec.0	List of number of sampled individuals for population 0 for each locus in the empirical distribution
sample.vec.1	List of number of sampled individuals for population 1 for each locus in the empirical distribution
sequence.length.vec	List of the number of sites for each locus in the empirical dataset (can be set to c(190,190) to simulate all loci with length of 190 nucleotides)

Examples

## Not run: 
library(GppFst)
library(Geneland)
library(phybase)

experimental_params <- read.table(file = '~/Desktop/GppFST_Tutorial/ExperimentalParameters.txt', header = T) # Read tab-delimited file with experimental parameters
pop0.samples <- experimental_params$pop0.samples # Extract pop0 samples per empirical locus
pop1.samples <- experimental_params$pop1.samples # Extract pop1 samples per empirical locus
locus.lengths <- experimental_params$locus.length # Extract locus lengths per empirical locus

MCMC.samples <- read.posterior(posterior.file = '~/Desktop/GppFST_Tutorial/atrox_snap_gamma2.log', format = "tab", burnin = .95)

Gppfst.results <- GppFst(posterior.samples = MCMC.samples, loci.per.step = 10, SNP.per.locus = 1, sample.vec.0 = pop0.samples, sample.vec.1 = pop1.samples, sequence.length.vec = locus.lengths)
mean(as.numeric(Gppfst.results) # get mean
hist(as.numeric(Gppfst.results) # plot distribution

Compare empirical Fst vs. simulated Fst
emp.fst <- experimental_params$WEIR_AND_COCKERHAM_FST

par(mfrow=c(2,1))

hist(as.numeric(Gppfst.results) # plot distribution
hist(as.numeric(emp.fst) # plot distribution

## End(Not run)

radamsRHA/GppFst documentation built on Nov. 9, 2019, 7:08 p.m.