Function to calculate Fst from Relate trees for pairs of populations specified in poplabels.

This function will calculate Fst in blocks of prespecified size for all pairs of populations specified in the poplabels file. Please refer to the Relate documentation for input file formats (https://myersgroup.github.io/relate/). The output is in the same format as for f2_blocks_from_Relate.

Fst_blocks_from_Relate(
  file_anc,
  file_mut,
  poplabels,
  file_map = NULL,
  chrs = NULL,
  blgsize = NULL,
  mu = NULL,
  tmin = NULL,
  t = NULL,
  transitions = NULL,
  use_muts = NULL,
  minMAF = NULL,
  Fst = NULL,
  dump_blockpos = NULL,
  apply_corr = NULL
)

Arguments

file_anc: Filename of anc file. If chrs is specified, this should only be the prefix, resulting in filenames of ${file_anc}_chr${chr}.anc(.gz).
file_mut: Filename of mut file. If chrs is specified, this should only be the prefix, resulting in filenames of ${file_anc}_chr${chr}.anc(.gz).
poplabels: Filename of poplabels file
file_map: (Optional) File prefix of recombination map. Not needed if blgsize is given in base-pairs, i.e. blgsize > 100
chrs: (Optional) Vector of chromosome IDs
blgsize: (Optional) SNP block size. Default is 10000 bases. If blgsize is <100 then it will be interpreted in Morgan. If blgsize is 100 or greater, if will be interpreted as base pair distance. If blgsize is negative, every tree is its own block.
mu: (Optional) Per base per generation mutation rate to scale f2 values. Default: 1.25e-8
tmin: (Optional) Minimum time cutof in generations. Any lineages younger than tmin will be excluded from the analysis. Default: t = 0.
t: (Optional) Time cutoff in generations. Default: Inf
transitions: (Optional) Set this to FALSE to exclude transition SNPs. Only meaningful with use_muts
use_muts: (Optional) Calculate traditional f2 statistics by only using mutations mapped to Relate trees. Default: false.
minMAF: (Optional) Minimum frequency cutoff. Default: 1 (i.e. excl singletons)
dump_blockpos: (Optional) Filename of blockpos file.
apply_corr: (Optional) Use small sample size correction. Default: true.

Value

3d array of dimension #groups x #groups x #blocks. Analogous to output of f2_from_geno in admixtools.

Examples

file_anc  <- system.file("sim/msprime_ad0.8_split250_1_chr1.anc.gz", package = "twigstats")
file_mut  <- system.file("sim/msprime_ad0.8_split250_1_chr1.mut.gz", package = "twigstats")
poplabels <- system.file("sim/msprime_ad0.8_split250_1.poplabels", package = "twigstats")
file_map  <- system.file("sim/genetic_map_combined_b37_chr1.txt.gz", package = "twigstats")

#Calculate f2s between all pairs of populations
Fst_blocks <- Fst_blocks_from_Relate(file_anc, file_mut, poplabels, file_map)

#Use a cutoff of 500 generations
Fst_blocks <- Fst_blocks_from_Relate(file_anc, file_mut, poplabels, file_map, dump_blockpos = "test.pos", t = 500)