Skip to contents

This function computes pairwise Fst. The following methods are implemented:

  • 'Hudson': Hudson's formulation, as derived in Bhatia et al (2013) for diploids. This is the only method that can also be used with pseudohaploid data.

  • 'Nei87' : Fst according to Nei (1987) - includes the correction for heterozygosity when computing Ht (it uses the same formulation as in hierfstat::pairwise.neifst()),

  • 'WC84' : Weir and Cockerham (1984), correcting for missing data (it uses the same formulation as in hierfstat::pairwise.WCfst()).

Usage

pairwise_pop_fst(.x, ...)

# S3 method for class 'gen_tbl'
pairwise_pop_fst(
  .x,
  type = c("tidy", "pairwise"),
  by_locus = FALSE,
  by_locus_type = c("tidy", "matrix", "list"),
  method = c("Hudson", "Nei87", "WC84"),
  return_num_dem = FALSE,
  n_cores = bigstatsr::nb_cores(),
  ...
)

Arguments

.x

a grouped gen_tibble (as obtained by using dplyr::group_by())

...

other arguments passed to specific methods, currently unused.

type

type of object to return One of "tidy" or "pairwise" for a pairwise matrix of populations. Default is "tidy". If data are grouped by more than one variable, the "tidy" output has one column per grouping variable for each population of a pair (e.g. population_1, sex_1, population_2, sex_2), and populations are otherwise labelled by the values of the grouping variables pasted with "." (e.g. pop1.F).

by_locus

boolean, determining whether Fst should be returned by locus(TRUE), or as a single genome wide value obtained by taking the ratio of the mean numerator and denominator (FALSE, the default).

by_locus_type

type of object to return. One of "tidy", "matrix" or "list". Default is "tidy".

method

one of 'Hudson', 'Nei87', and 'WC84'

return_num_dem

returns a list of numerators and denominators for each locus. This is useful for creating windowed estimates of Fst (as we need to compute the mean numerator and denominator within each window). Default is FALSE.

n_cores

number of cores to be used, it defaults to bigstatsr::nb_cores()

Value

if type=tidy, a tibble of genome-wide pairwise Fst values with each pairwise combination as a row if "by_locus=FALSE", else a list including the tibble of genome-wide values as well as a matrix with pairwise Fst by locus with loci as rows and and pairwise combinations as columns. If type=pairwise, a matrix of genome-wide pairwise Fst values is returned.

Details

For all formulae, the genome wide estimate is obtained by taking the ratio of the mean numerators and denominators over all relevant SNPs.

Loci with no genotypes in one of the two populations of a pair (e.g. loci that failed in all the individuals of a population) carry no information on the differentiation between them, so they are dropped for that pair with all methods: their Fst by locus is NaN, and they do not contribute to the mean numerator and denominator of the genome wide estimate.

Difference from hierfstat for 'WC84': hierfstat::pairwise.WCfst() (and hierfstat::wc()) keep these loci. Their among-population variance component (\(\sigma^2_a\), the numerator) is set to 0, while their within-population components (\(\sigma^2_b\) and \(\sigma^2_w\), estimated from the one population with data) still add to the denominator. So hierfstat gives these loci an Fst of 0 rather than NaN, and its genome wide estimate is pulled towards 0 (i.e. it is usually lower than the one computed here). The two estimates are identical when every locus has data in both populations, or if such loci are removed before using hierfstat. This applies to every function that computes Fst with 'WC84' (e.g. windows_pairwise_pop_fst(), jackknife_pairwise_pop_fst() and nwise_pop_pbs()). 'Nei87' gives the same results as hierfstat::pairwise.neifst(), which also drops these loci.

References

Bhatia G, Patterson N, Sankararaman S, Price AL. (2013) Estimating and Interpreting FST: The Impact of Rare Variants. Genome Research, 23(9):1514–1521.

Nei, M. (1987) Molecular Evolutionary Genetics. Columbia University Press

Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution, 38(6): 1358–1370.

Examples

example_gt <- load_example_gt("gen_tbl")

# For a basic global pairwise Fst calculation:
example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87")
#> # A tibble: 3 × 3
#>   population_1 population_2  value
#>   <chr>        <chr>         <dbl>
#> 1 pop1         pop2         0.0320
#> 2 pop1         pop3         0.143 
#> 3 pop2         pop3         0.0500

# With a pairwise matrix:
example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87", type = "pairwise")
#>           pop1      pop2      pop3
#> pop1        NA 0.0320197 0.1428571
#> pop2 0.0320197        NA 0.0500000
#> pop3 0.1428571 0.0500000        NA
#> attr(,"class")
#> [1] "pairwise_matrix" "matrix"          "array"          

# To calculate Fst by locus:
example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Hudson", by_locus = TRUE)
#> $Fst_by_locus
#> # A tibble: 18 × 3
#>    loci  stat_name       value
#>    <chr> <chr>           <dbl>
#>  1 rs1   fst_pop1.pop2  -0.24 
#>  2 rs1   fst_pop1.pop3   0.6  
#>  3 rs1   fst_pop2.pop3   0.667
#>  4 rs2   fst_pop1.pop2   0.6  
#>  5 rs2   fst_pop1.pop3   0.114
#>  6 rs2   fst_pop2.pop3   0    
#>  7 rs3   fst_pop1.pop2 NaN    
#>  8 rs3   fst_pop1.pop3   0    
#>  9 rs3   fst_pop2.pop3   0    
#> 10 rs4   fst_pop1.pop2  -0.167
#> 11 rs4   fst_pop1.pop3   0.333
#> 12 rs4   fst_pop2.pop3   0    
#> 13 rs5   fst_pop1.pop2  -0.1  
#> 14 rs5   fst_pop1.pop3  -0.6  
#> 15 rs5   fst_pop2.pop3  -0.5  
#> 16 rs6   fst_pop1.pop2  -0.25 
#> 17 rs6   fst_pop1.pop3  -0.333
#> 18 rs6   fst_pop2.pop3  -0.5  
#> 
#> $Fst
#> # A tibble: 3 × 3
#>   population_1 population_2  value
#>   <chr>        <chr>         <dbl>
#> 1 pop1         pop2         0.0345
#> 2 pop1         pop3         0.0556
#> 3 pop2         pop3         0     
#>